anchor.diff:      forall(i=1:grid%nC)
         y(i)=rhoE0(i)+dt/grid%v(i)*(flowRhoH(i)+condQ(i)) ! TODO add viscous heating
       end forall
  -    energyRHS=0
  +    energyRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
     end function
     
   end module
  
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:     
      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
  