anchor.diff:      end do
     end if
   
  -  ! Now add this to the thickness tendency due to the flow field and
  -  ! sponge regions
  -  dhdt = 0d0
  +  return
  +end subroutine dhdt_vert_diff
  +
  +! ---------------------------------------------------------------------------
  +!> Use a first-order centered advection scheme to calculate the advctive
  +!! thickness tendency
  +
  +subroutine h_advec_1st_centered(dhdt_advec, h, u, v, dx, dy, nx, ny, layers)
  +  implicit none
  +
  +  ! dhdt is evaluated at the centre of the grid box
  +  double precision, intent(out) :: dhdt_advec(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in)  :: h(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in)  :: u(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in)  :: v(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in)  :: dx, dy
  +  integer, intent(in) :: nx, ny, layers
  +
  +  integer i, j, k
  +
  +  dhdt_advec = 0d0
   
     do k = 1, layers
       do j = 1, ny
         do i = 1, nx
  -        dhdt(i,j,k) = &
  -            dhdt_GM(i,j,k) & ! horizontal thickness diffusion
  -            + dhdt_vert_diff(i,j,k) & ! vetical thickness diffusion 
  +        dhdt_advec(i,j,k) = & 
               - ((h(i,j,k)+h(i+1,j,k))*u(i+1,j,k) &
  -               - (h(i-1,j,k)+h(i,j,k))*u(i,j,k))/(dx*2d0) & ! d(hu)/dx
  +             - (h(i-1,j,k)+h(i,j,k))*u(i,j,k))/(dx*2d0) & ! d(hu)/dx
               - ((h(i,j,k)+h(i,j+1,k))*v(i,j+1,k) &
  -              - (h(i,j-1,k)+h(i,j,k))*v(i,j,k))/(dy*2d0)  & ! d(hv)/dy
  -            + spongeTimeScale(i,j,k)*(spongeH(i,j,k)-h(i,j,k)) ! forced relaxtion in the sponge regions.
  +             - (h(i,j-1,k)+h(i,j,k))*v(i,j,k))/(dy*2d0)   ! d(hv)/dy
         end do
       end do
     end do
   
  -  ! Make sure the dynamics are only happening in the wet grid points.
  -  do k = 1, layers
  -    dhdt(:, :, k) = dhdt(:, :, k) * wetmask
  -  end do
  -
  -  call wrap_fields_3D(dhdt, nx, ny, layers)
  -
  -  return
  -end subroutine evaluate_dhdt
  +  return
  +end subroutine h_advec_1st_centered
   
   ! ---------------------------------------------------------------------------
   !> Calculate the tendency of zonal velocity for each of the active layers
  
positive.code:     dhdt_GM(:,:,layers) = -sum(dhdt_GM(:,:,:layers-1), 3)
    end if
  
    ! calculate vertical diffusive mass fluxes
    ! only evaluate vertical mass diff flux if more than 1 layer, or reduced gravity
    if (layers .eq. 1) then
      if (RedGrav) then
        do j = 1, ny
          do i = 1, nx
            dhdt_vert_diff(i,j,1) = kv/h(i,j,1)
          end do
        end do
      end if
  
negative.code:       do j = 1, ny
          do i = 1, nx
            if (k .eq. 1) then ! in top layer
              dhdt_kv(i,j,k) = kv/h(i,j,k) - kv/h(i,j,k+1)
            else if (k .eq. layers) then ! bottom layer
              dhdt_kv(i,j,k) = kv/h(i,j,k) - kv/h(i,j,k-1)
            else ! mid layer/s
              dhdt_kv(i,j,k) = 2d0*kv/h(i,j,k) -  &
                  kv/h(i,j,k-1) - kv/h(i,j,k+1)
            end if
          end do
  