anchor.diff:          dudt_bt = 0d0
           dvdt_bt = 0d0
   
  -        do i = 2,nx-1
  +        do i = 1,nx-1
               do j = 1,ny-1
  -                dudt_bt(i,j) = -g_vec(1)*(etanew(i,j) - etanew(i-1,j))/(dx)
  -                unew(i,j,:) = unew(i,j,:) + 23d0*dt*dudt_bt(i,j)/12d0
  +                do k = 1,layers
  +                    dudt_bt(i,j) = -g_vec(1)*(etanew(i,j) - etanew(i-1,j))/(dx)
  +                    unew(i,j,k) = unew(i,j,k) + 23d0*dt*dudt_bt(i,j)/12d0
  +                end do
               end do
           end do
   
           do i = 1,nx-1
  -            do j = 2,ny-1
  -                dvdt_bt(i,j) = -g_vec(1)*(etanew(i,j) - etanew(i,j-1))/(dy)
  -                vnew(i,j,:) = vnew(i,j,:) + 23d0*dt*dvdt_bt(i,j)/12d0
  +            do j = 1,ny-1
  +                do k = 1,layers
  +                    dvdt_bt(i,j) = -g_vec(1)*(etanew(i,j) - etanew(i,j-1))/(dy)
  +                    vnew(i,j,k) = vnew(i,j,k) + 23d0*dt*dvdt_bt(i,j)/12d0
  +                end do
               end do
           end do
   
           ! Add the barotropic tendency to the baroclinic tendency so that it gets used in subsequent time steps.
  -        do i = 2,nx-1
  +        do i = 1,nx-1
               do j = 1,ny-1
                   do k = 1,layers
                       dudt(i,j,k) = dudt(i,j,k) + dudt_bt(i,j)
  
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:     b = 0d0
  !   no baroclinic pressure contribution to the first layer Bernoulli potential 
  !   (the barotropic pressure contributes, but that's not done here).
      ! do j = 1,ny-1
      !     do i = 1,nx-1
      !         b(i,j,1) = (u(i,j,1)**2+u(i+1,j,1)**2+v(i,j,1)**2+v(i,j+1,1)**2)/4.0d0
      !     end do
      ! end do
  
  !   for the rest of the layers we get a baroclinic pressure contribution
      do k = 1,layers !move through the different layers of the model
        do j=1,ny-1 !move through longitude
          do i=1,nx-1 ! move through latitude
            b(i,j,k)= M(i,j,k) + (u(i,j,k)**2+u(i+1,j,k)**2+v(i,j,k)**2+v(i,j+1,k)**2)/4.0d0
  