anchor.diff:        y(i)=p(i)-R_AIR*temp(i)*dt**2/grid%v(i)*(laP(i)-laP1(i))-R_AIR*temp(i)/R_AIR/temp1(i)*p1(i)&
         &    +R_AIR*temp(i)*(rho1(i)-rho0(i))-R_AIR*temp(i)*dt/grid%v(i)*tmpFlowRho(i)
       end forall
  -    pressureRHS=0
  +    pressureRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
     end function
     
     !> RHS of the density equation in the form of a fix point problem
     function densityRHS(oldVector,newVector,dat)
       use iso_c_binding
  +    use ieee_arithmetic
       use modNumerics
       use modAdvection
       use modRhieChow
  
positive.code:     
      p(1:grid%nC)=x(1:grid%nC)
      call setBC()
      if(.not.allocated(tmpFlowRho))then
        allocate(tmpFlowRho(grid%nC))
      else if(size(tmpFlowRho)/=grid%nC)then
  
negative.code:   !> RHS of the momentum equation in the form of a fix point problem
    function momentumRHS(oldVector,newVector,dat)
      use iso_c_binding
      use ieee_arithmetic
      use modNumerics
      use modAdvection
      use modNewtonian
      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
  