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: !! 0 means barrier
  !! 1 mean open
  
  subroutine calc_boundary_masks(wetmask,hfacW,hfacE,hfacS,hfacN,nx,ny)
    implicit none
  
    integer nx !< number of grid points in x direction
    integer ny !< number of grid points in y direction
    double precision wetmask(0:nx+1,0:ny+1)
    double precision temp(0:nx+1,0:ny+1)
    double precision hfacW(0:nx+1,0:ny+1)
    double precision hfacE(0:nx+1,0:ny+1)
    double precision hfacN(0:nx+1,0:ny+1)
    double precision hfacS(0:nx+1,0:ny+1)
    integer i,j
  
    hfacW = 1d0
  
    temp = 0.0
    do j = 0,ny+1
      do i = 1,nx+1
        temp(i,j) = wetmask(i-1,j)- wetmask(i,j)
      end do
    end do
  
    do j = 0,ny+1
      do i = 1,nx+1
        if (temp(i,j) .ne. 0.0) then
          hfacW(i,j) = 0d0
        end if
      end do
    end do
  
    ! and now for all  western cells
    hfacW(0,:) = hfacW(nx,:)
  
    hfacE = 1d0
  
    temp = 0.0
    do j = 0,ny+1
      do i = 0,nx
        temp(i,j) = wetmask(i,j)- wetmask(i+1,j)
      end do
    end do
  
    do j = 0,ny+1
      do i = 0,nx
        if (temp(i,j) .ne. 0.0) then
          hfacE(i,j) = 0d0
        end if
      end do
    end do
  
    ! and now for all  eastern cells
    hfacE(nx+1,:) = hfacE(1,:)
  
    hfacS = 1
  
    temp = 0.0
    do j = 1,ny+1
      do i = 0,nx+1
        temp(i,j) = wetmask(i,j-1)- wetmask(i,j)
      end do
    end do
  
    do j = 1,ny+1
      do i = 0,nx+1
        if (temp(i,j) .ne. 0.0) then
          hfacS(i,j) = 0d0
        end if
      end do
    end do
  
    ! all southern cells
    hfacS(:,0) = hfacS(:,ny)
  
    hfacN = 1
    temp = 0.0
    do j = 0,ny
      do i = 0,nx+1
        temp(i,j) = wetmask(i,j)- wetmask(i,j+1)
      end do
    end do
  
    do j = 0,ny
      do i = 0,nx+1
        if (temp(i,j) .ne. 0.0) then
          hfacN(i,j) = 0d0
        end if
      end do
    end do
    ! all northern cells
    hfacN(:,ny+1) = hfacN(:,1)
  
    return
  end subroutine calc_boundary_masks
  
  ! -----------------------------------------------------------------------------
  
  subroutine read_input_fileH(name,array,default,nx,ny,layers)
    implicit none
  
    character(30) name
    integer nx, ny, layers, k
    double precision array(0:nx+1,0:ny+1,layers), default(layers)
    double precision array_small(nx,ny,layers)
  
    if (name.ne.'') then
      open(unit = 10, form='unformatted', file=name)
      read(10) array_small
      close(10)
  
      array(1:nx,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
      do k = 1,layers
        array(:,:,k) = default(k)
      end do
    end if
  
  
negative.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
  