anchor.diff:        include 'pom.h'
         integer i,j,k
   
  -      do k=1,kb
  -        do j=1,jm
  -          do i=1,im
  -            rho(i,j,k)=rho(i,j,k)-rmean(i,j,k)
  -          end do
  -        end do
  -      end do
  +      rho = rho-rmean
         
   ! calculate x-component of baroclinic pressure gradient
         do j=2,jmm1
  
positive.code: 
  ! calculate y-component of velocity advection
  
        do k=1,kb
          do j=1,jm
            do i=1,im
              advy(i,j,k)=0.d0
              xflux(i,j,k)=0.d0
              yflux(i,j,k)=0.d0
            end do
          end do
        end do
  
  ! calculate horizontal advective fluxes
        do k=1,kbm1
  
negative.code:       integer i,j,k
  
  ! do vertical advection
        uf = 0.
  
        do k=2,kbm1
          do j=1,jm
  