anchor.diff:         endif
   
          do rb = rbb,rbe ! Red black loop
  - 
  +
             do i = ib,ie
                do j = jb+mod(i+rb,rbi),je,rbi
  -              
  +
  +!!$       ib = 1
  +!!$       ie = nx
  +!!$       jb = 1
  +!!$       je = ny
  +!!$
  +!!$       fcb = 1
  +!!$       fce = 2
  +!!$       fci = 2
  +!!$
  +!!$       do fc1 = fcb, fce ! 
  +!!$
  +!!$          do fc2 = fcb, fce ! 
  +!!$
  +!!$             do i = ib+mod(fc1-1,fci),ie,fci
  +!!$
  +!!$                do j = jb+mod(fc2-1,fci),je,fci
  +
                      k=1 !lower level
                      rhs(k) = b(k,j,i)                                              &
                           - cA(3,k,j,i)*p(k+1,j-1,i)                                &
  
positive.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) &
  
negative.code:                         - 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
  