anchor.diff:      allocate(H(grid%nC))
       allocate(massFlowSensU(DIMS,2,grid%nP))
       allocate(massFlowSensP(2,grid%nP))
  -    allocate(pbcDiagPrec((DIMS+1)*grid%nC))
  +    allocate(pbcPrecI((size(grid%neib,1)+1)*(DIMS+1)*grid%nC))
  +    allocate(pbcPrecJ((size(grid%neib,1)+1)*(DIMS+1)*grid%nC))
  +    allocate(pbcPrecVal((size(grid%neib,1)+1)*(DIMS+1)*grid%nC))
       allocate(pbcXFact((DIMS+1)*grid%nC))
       allocate(pbcRFact((DIMS+1)*grid%nC))
       allocate(fullXFact((DIMS+2)*grid%nC))
  
positive.code:   
    ! data for algebraic solver
    type(NewtonKrylov)::pbcEq !< momentum and pressure equations (isothermal) solved by Newton-GMRES
    type(fixPt)::energyEq !< energy equation as a fix point problem
    type(NewtonKrylov)::fullEq !< full NS equations solved by Newton-GMRES
    integer::nItPBC,nItEnergy,nItOuter,nItFull !< number of iterations
  
negative.code:     call associateVector(c_fScale,fScale)
      call associateVector(c_v,v)
      
      if(.not.allocated(pbcRst))then
        allocate(pbcRst(grid%nC*(DIMS+1)))
        allocate(pbcRHS(grid%nC*(DIMS+1)))
      end if
      
      ! find x increment from f increment v and save back to v
      pbcRHS(1:grid%nC*(DIMS+1))=v(1:grid%nC*(DIMS+1))
      call pbcPrec%solve(pbcRHS,pbcRst)
      v(1:grid%nC*(DIMS+1))=pbcRst(1:grid%nC*(DIMS+1))
      pbcPsol=0
      foo=dat
    end function
  