anchor.diff:                      dvdt(i,j,k) = dvdt(i,j,k) - 1.0d0*ar*(v(i,j,k) - 1.0d0*v(i,j,k+1))
                   else if (k .eq. layers) then ! bottom layer
                       dvdt(i,j,k) = dvdt(i,j,k) - 1.0d0*ar*(v(i,j,k) - 1.0d0*v(i,j,k-1))
  +                    ! add bottom drag here in isopycnal version
                   else ! mid layer/s
                       dvdt(i,j,k) = dvdt(i,j,k) - 1.0d0*ar*(2.0d0*v(i,j,k) - 1.0d0*v(i,j,k-1) - 1.0d0*v(i,j,k+1))
                   endif
  
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:         end do
  
          do i = 1,nx-1
              do j = 1,ny-1
                  do k = 1,layers
                      dvdt(i,j,k) = dvdt(i,j,k) + dvdt_bt(i,j)
                  end do
  