anchor.diff:      
       p(1:grid%nC)=x(1:grid%nC)
       call setBC()
  +    call findGrad(grid,p,gradP)
       if(.not.allocated(tmpFlowRho))then
         allocate(tmpFlowRho(grid%nC))
       else if(size(tmpFlowRho)/=grid%nC)then
  
positive.code:   !> RHS of the momentum equation in the form of a fix point problem
    function momentumRHS(oldVector,newVector,dat)
      use iso_c_binding
      use modNumerics
      use modAdvection
      use modNewtonian
      use modRhieChow
      type(C_PTR),value::oldVector !< old N_Vector
      type(C_PTR),value::newVector !< new N_Vector
      type(C_PTR),value::dat !< optional user data object
  
negative.code:     end forall
      call setBC()
      call findViscForce(grid,u,visc,viscF)
      viscF(:,:)=0d0 ! FIXME: enable viscous force
      forall(i=1:grid%nE,j=1:DIMS)
        fluxRhou(:,j,i)=rhou(j,i)*u(:,i)
      end forall
      call findAdv(grid,rhou,fluxRhou,flowRhou)
      forall(i=1:grid%nC)
        rhou(:,i)=rhou0(:,i)+dt/grid%v(i)*(flowRhou(:,i)+presF(:,i)+viscF(:,i))
      end forall
      y=reshape(rhou(:,1:grid%nC),[grid%nC*DIMS])
      momentumRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
    end function
    
    !> residual function of the pressure equation
    function pressureRHS(oldVector,newVector,dat)
      use iso_c_binding
      use ieee_arithmetic
      use modNumerics
      use modGradient
      use modDiffusion
      use modRhieChow
      type(C_PTR),value::oldVector !< old N_Vector
  