anchor.diff:      dhdt_GM(:,:,layers) = -sum(dhdt_GM(:,:,:layers-1), 3)
     end if
   
  +  return
  +end subroutine dhdt_hor_diff
  +
  +! ---------------------------------------------------------------------------
  +!> Calculate the tendency of layer thickness for each of the active layers
  +!! due to vertical thickness diffusion
  +subroutine dhdt_vert_diff(dhdt_kv, h, kv, nx, ny, layers, RedGrav)
  +  implicit none
  +
  +  ! dhdt is evaluated at the centre of the grid box
  +  double precision, intent(out) :: dhdt_kv(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in)  :: h(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in)  :: kv
  +  integer, intent(in) :: nx, ny, layers
  +  logical, intent(in) :: RedGrav
  +
  +  integer i, j, k
  +
  +  dhdt_kv = 0d0
  +
     ! 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)
  +          dhdt_kv(i,j,1) = kv/h(i,j,1)
           end do
         end do
       end if
  
positive.code:   integer i, j, k
    ! Thickness tendency due to thickness diffusion (equivalent to Gent
    ! McWilliams in a z coordinate model)
    double precision dhdt_GM(0:nx+1, 0:ny+1, layers)
  
    ! Thickness tendency due to vertical diffusion of mass
    double precision dhdt_vert_diff(0:nx+1, 0:ny+1, layers)
  
    ! Calculate tendency due to thickness diffusion (equivalent
    ! to GM in z coordinate model with the same diffusivity).
    dhdt_GM = 0d0
    dhdt_vert_diff = 0d0
  
    ! Loop through all layers except lowest and calculate
    ! thickness tendency due to horizontal diffusive mass fluxes
  
negative.code:     end do
    end if
  
    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_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,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
        end do
      end do
    end do
  
    return
  end subroutine h_advec_1st_centered
  
  ! ---------------------------------------------------------------------------
  !> Calculate the tendency of zonal velocity for each of the active layers
  