anchor.diff:  
                      d(k)   = cA(1,k,j,i)
                      ud(k)  = cA(2,k+1,j,i)
  - 
  +
                      do k = 2,nz-1 !interior levels
                         rhs(k) = b(k,j,i) &
                              - cA(3,k,j,i)*p(k+1,j-1,i) - cA(3,k-1,j+1,i)*p(k-1,j+1,i) &
  
positive.code:     !
      !     LOCAL 
      integer(kind=ip)            :: i,j,k,it,rb
      integer(kind=ip)            :: ib,ie,jb,je,rbb,rbe,rbi
      real(kind=rp),dimension(nz) :: rhs,d,ud
  
      real(kind=rp)    :: g1,g2,gamma
  
negative.code:                    ! seems to be the same
                     ! obviously, this residual is only an estimate, it is not the correct one
  
                  enddo ! j
               enddo    ! i
  
               ! don't call mpi at every pass if nh>1
               !if ((mod(it,nh) == 0).or.(it==nsweeps)) then
               call fill_halo(lev,p)
               !endif
  
  !!$          enddo  ! fc2
  !!$
  !!$       enddo  ! fc1
  
         enddo ! Red Black
  
      enddo  !it
  
      call toc(lev,'relax_line')
  
  