anchor.diff:      call findGrad(grid,u,gradU)
       call findPresForce(grid,p,gradP,presF)
       ! non-linear functional iterations on disturbed u and p (analogous to Jacobi iterations)
  -    do l=1,10
  +    do l=1,20
         call setBC()
         !call findGrad(grid,u,gradU)
         !call findViscForce(grid,u,gradU,visc,viscF)
  
positive.code:       close(FID)
      end if
      ! initialize algebraic solver
      call pbcEq%init(grid%nC*(DIMS+1),pbcRHS,maxl=MAXIT_PBC,pset=pbcPSet,psol=pbcPSol)
      call pbcEq%setMaxIt(MAXIT_PBC)
      call energyEq%init(grid%nC,energyRHS,maa=MAXIT_ENERGY)
      call energyEq%setMaxIt(MAXIT_ENERGY)
  
negative.code:   
    integer,parameter::DIMS=3 !< three dimensions
    
    integer,parameter::MAXIT_PBC=50 !< max number of momentum and pressure equation iterations
    integer,parameter::MAXL_PBC=20 !< max dimensions of PBC Krylov subspace
    integer,parameter::MAXIT_ENERGY=50 !< max number of energy equation iterations
    integer,parameter::MAXIT_OUTER=20 !< max number of outer iterations
    integer,parameter::MAXIT_FULL=50 !< max number of full NS equation iterations
  