anchor.diff:      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))))
  +    if(c_associated(dat))then
  +    end if
     end function
     
     !> residual function of the pressure equation
  
positive.code:       TInit=298d0
        YInit=[1d0]
        ! FIXME remove above
        ! FIXME calculation of rho, rhoE based on cantera
        p(i)=pInit
        u(:,i)=uInit(:)
        temp(i)=TInit
  
negative.code:       y(i)=rhoE0(i)+dt/grid%v(i)*(flowRhoH(i)+condQ(i)) ! TODO add viscous heating
      end forall
      energyRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
      if(c_associated(dat))then
      end if
    end function
    
  end module
  