anchor.diff:        call order3d_mpi(rho,rho4th,im_local,jm_local,kb)
   
   ! compute terms correct to 4th order
  -      do i=1,im
  -        do j=1,jm
  -          ddx(i,j)=0.
  -          d4(i,j)=0.
  -        end do
  -      end do
  -      do k=1,kb
  -        do j=1,jm
  -          do i=1,im
  -            rhou(i,j,k)=0.
  -            drho(i,j,k)=0.
  -          end do
  -        end do
  -      end do
  +      ddx = 0.
  +      d4  = 0.
  +      rhou = 0.
  +      drho = 0.
   
   ! compute DRHO, RHOU, DDX and D4
         do j=1,jm
  
positive.code:       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
        
  ! calculate x-component of baroclinic pressure gradient
        do j=2,jmm1
  
negative.code:       double precision xflux(im,jm,kb),yflux(im,jm,kb)
        integer i,j,k
  
        xflux = 0.
        yflux = 0.
  
  ! reestablish boundary conditions
        do k=1,kbm1
  