anchor.diff:      end do
     end do
   
  -  call wrap_fields_3D(dvdt,nx,ny,layers)
  +  call wrap_fields_3D(dvdt, nx, ny, layers)
   
     return
   end subroutine evaluate_dvdt
   
   ! -----------------------------------------------------------------------------
   !> Calculate the barotropic u velocity
  -subroutine calc_baro_u(ub,u,h,eta,freesurfFac,nx,ny,layers)
  -  implicit none
  -
  -  integer nx,ny,layers
  -  integer i,j,k
  -  double precision h(0:nx+1,0:ny+1,layers)
  -  double precision h_temp(0:nx+1,0:ny+1,layers)
  -  double precision eta(0:nx+1,0:ny+1)
  +subroutine calc_baro_u(ub, u, h, eta, freesurfFac, nx, ny, layers)
  +  implicit none
  +
  +  integer nx, ny, layers
  +  integer i, j, k
  +  double precision h(0:nx+1, 0:ny+1, layers)
  +  double precision h_temp(0:nx+1, 0:ny+1, layers)
  +  double precision eta(0:nx+1, 0:ny+1)
     double precision freesurfFac
  -  double precision u(0:nx+1,0:ny+1,layers)
  -  double precision ub(0:nx+1,0:ny+1)
  +  double precision u(0:nx+1, 0:ny+1, layers)
  +  double precision ub(0:nx+1, 0:ny+1)
   
     ub = 0d0
   
     h_temp = h
     ! add free surface elevation to the upper layer
  -  h_temp(:,:,1) = h(:,:,1) + eta*freesurfFac
  -
  -  do i = 1,nx
  -    do j = 1,ny
  -      do k = 1,layers
  +  h_temp(:, :, 1) = h(:, :, 1) + eta*freesurfFac
  +
  +  do i = 1, nx
  +    do j = 1, ny
  +      do k = 1, layers
           ub(i,j) = ub(i,j) + u(i,j,k)*(h_temp(i,j,k)+h_temp(i-1,j,k))/2d0
         end do
       end do
     end do
   
  -  call wrap_fields_2D(ub,nx,ny)
  +  call wrap_fields_2D(ub, nx, ny)
   
     return
   end subroutine calc_baro_u
  
positive.code: !> Calculate the free surface anomaly using the velocities
  !! timestepped with the tendencies excluding the free surface
  !! pressure gradient.
  subroutine calc_eta_star(ub,vb,eta,etastar,freesurfFac,nx,ny,dx,dy,dt)
    implicit none
  
    integer nx,ny,layers
    integer i,j,k
    double precision eta(0:nx+1,0:ny+1)
    double precision etastar(0:nx+1,0:ny+1)
    double precision ub(0:nx+1,0:ny+1)
    double precision vb(0:nx+1,0:ny+1)
    double precision freesurfFac, dx,dy,dt
  
    etastar = 0d0
  
    do i = 1,nx
      do j = 1,ny
        etastar(i,j) = freesurfFac*eta(i,j) - &
            dt*((ub(i+1,j) - ub(i,j))/dx + (vb(i,j+1) - vb(i,j))/dy)
      end do
    end do
  
    call wrap_fields_2D(etastar,nx,ny)
  
    return
  end subroutine calc_eta_star
  
negative.code: 
  ! -----------------------------------------------------------------------------
  
  subroutine read_input_fileU(name, array, default, nx, ny, layers)
    implicit none
  
    character(30) name
    integer nx, ny, layers
    double precision array(0:nx+1, 0:ny+1, layers), default
    double precision array_small(nx+1, ny, layers)
  
    if (name.ne.'') then
      open(unit = 10, form='unformatted', file=name)
      read(10) array_small
      close(10)
  
      array(1:nx+1, 1:ny, :) = array_small
      ! wrap array around for periodicity
      array(0, :, :) = array(nx, :, :)
      array(nx+1, :, :) = array(1, :, :)
      array(:, 0, :) = array(:, ny, :)
      array(:, ny+1, :) = array(:, 1, :)
    else
      array = default
    end if
  