anchor.diff:                          - cA(7,k,j,i)*p(k  ,j,i-1) - cA(7,k  ,j,i+1)*p(k  ,j,i+1) &
                           - cA(8,k,j,i)*p(k-1,j,i-1) 
                      d(k)   = cA(1,k,j,i)
  - 
  +
                      call tridiag(nz,d,ud,rhs,p(:,j,i)) !solve for vertical_coeff_matrix.p1d=rhs
   
                      ! December 10th, dev below is to try to by-pass the computation of the residual
  
positive.code: 
      do i = 1,nx
         do j = 1,ny
            
               z = b(k,j,i) - cA(1,k,j,i)*p(k,j,i)                           &
                    - cA(2,k,j,i)*p(k,j-1,i  ) - cA(2,k,j+1,i  )*p(k,j+1,  i)&
                    - cA(3,k,j,i)*p(k,j  ,i-1) - cA(3,k,j  ,i+1)*p(k,j  ,i+1)&
                    - cA(4,k,j,i)*p(k,j-1,i-1) - cA(4,k,j+1,i+1)*p(k,j+1,i+1)&
                    - cA(5,k,j,i)*p(k,j+1,i-1) - cA(5,k,j-1,i+1)*p(k,j-1,i+1)
  
               r(k,j,i) = z
               !          res = max(res,abs(r(k,j,i)))
               res = res+z*z
  
         enddo
      enddo
  
negative.code: 
                     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) &
  