anchor.diff:          end do
   
           do i = 1,nx-1
  -            do j = 2,ny-1
  +            do j = 1,ny-1
                   do k = 1,layers
                       dvdt(i,j,k) = dvdt(i,j,k) + dvdt_bt(i,j)
                   end do
  
positive.code:     double precision eta(0:nx,0:ny)
      double precision freesurfFac
      double precision v(0:nx,ny,layers),vb(0:nx,ny)
      double precision proto_vb
  
      vb = 0d0
  
      ! add free surface elevation to the upper layer
      h_temp(:,:,1) = h(:,:,1) + eta*freesurfFac
      
      do i = 1,nx-1
          do j = 2,ny-1
              proto_vb = 0d0
              do k = 1,layers
                  proto_vb = proto_vb + v(i,j,k)*(h_temp(i,j,k)+h_temp(i,j-1,k))/2d0
              end do
              vb(i,j) = proto_vb
          end do
      end do
  
  
negative.code:                     dudt(i,j,k) = dudt(i,j,k) - 1.0d0*ar*(u(i,j,k) - 1.0d0*u(i,j,k+1))
                  else if (k .eq. layers) then ! bottom layer
                      dudt(i,j,k) = dudt(i,j,k) - 1.0d0*ar*(u(i,j,k) - 1.0d0*u(i,j,k-1))
                      ! add bottom drag here in isopycnal version
                  else ! mid layer/s
                      dudt(i,j,k) = dudt(i,j,k) - 1.0d0*ar*(2.0d0*u(i,j,k) - 1.0d0*u(i,j,k-1) - 1.0d0*u(i,j,k+1))
                  endif
  