anchor.diff:        end if
   
         ! Apply the boundary conditions
  -      ! - Enforce no normal flow boundary condition
  -      !   and no flow in dry cells.
  -      ! - no/free-slip is done inside the dudt and dvdt subroutines.
  -      ! - hfacW and hfacS are zero where the transition between
  -      !   wet and dry cells occurs.
  -      ! - wetmask is 1 in wet cells, and zero in dry cells.
  -      do k = 1, layers
  -        unew(:, :, k) = unew(:, :, k) * hfacW * wetmask(:, :)
  -        vnew(:, :, k) = vnew(:, :, k) * hfacS * wetmask(:, :)
  -      end do
  +      call apply_boundary_conditions(unew, hfacW, wetmask, nx, ny, layers)
  +      call apply_boundary_conditions(vnew, hfacS, wetmask, nx, ny, layers)
       end if
   
       ! Stop layers from getting too thin
  
positive.code:     hnew = h+dt*(23d0*dhdt - 16d0*dhdtold + 5d0*dhdtveryold)/12d0
  
      ! Apply the boundary conditions
      ! - Enforce no normal flow boundary condition
      !   and no flow in dry cells.
      ! - no/free-slip is done inside the dudt and dvdt subroutines.
      ! - hfacW and hfacS are zero where the transition between
      !   wet and dry cells occurs.
      ! - wetmask is 1 in wet cells, and zero in dry cells.
      do k = 1, layers
        unew(:, :, k) = unew(:, :, k) * hfacW * wetmask(:, :)
        vnew(:, :, k) = vnew(:, :, k) * hfacS * wetmask(:, :)
      end do
  
      ! Do the isopycnal layer physics
      if (.not. RedGrav) then
  
negative.code:   v = v + dt*dvdtold
  
    ! Apply the boundary conditions
    call apply_boundary_conditions(u, hfacW, wetmask, nx, ny, layers)
    call apply_boundary_conditions(v, hfacS, wetmask, nx, ny, layers)
  
    ! Wrap fields around for periodic simulations
    call wrap_fields_3D(u, nx, ny, layers)
  