anchor.diff:      ! solution and residual scales
       pScale=max(maxval(p)-minval(p),maxval(0.5d0*rho*norm2(u,1)**2),1d0)
       uScale=max(maxval(norm2(u,1)),sqrt(2d0*pScale/minval(rho)))
  -    tempScale=max(maxval(temp)-minval(temp),pScale/minval(rho)/287.058d0)
  +    tempScale=max(maxval(temp)-minval(temp),pScale/minval(rho)/Rgas)
       rhoScale=maxval(rho)
       rhouScale=rhoScale*uScale
  -    rhoEScale=max(rhoScale*287.058*tempScale,0.5d0*rhoScale*uScale**2)
  +    rhoEScale=max(rhoScale*Rgas*tempScale,0.5d0*rhoScale*uScale**2)
       
       ! set tolerance and maximum number of iterations, and etc.
       call momentumEq%setTol(rhoScale*RTOL_MOMENTUM)
  
positive.code:     integer(C_INT)::energyRHS !< error code
      double precision,pointer::x(:) !< fortran pointer associated with oldVector
      double precision,pointer::y(:) !< fortran pointer associated with newVector
      double precision,parameter::R_AIR=287.058d0 ! TODO fluid data
      
      call associateVector(oldVector,x)
      call associateVector(newVector,y)
      
      rhoE(1:grid%nC)=x(1:grid%nC)
      forall(i=1:grid%nE)
        temp(i)=(rhoE(i)/rho(i)-0.5d0*dot_product(u(:,i),u(:,i)))/(1d0/(1.4d0-1d0))/R_AIR
      end forall
      call setBC()
      forall(i=1:grid%nE)
  
negative.code:   public
    
    integer,parameter::DIMS=3 !< three dimensions
    
    integer,parameter::MAXIT_MOMENTUM=40 !< max number of momentum prediction equation iteration
    integer,parameter::MAXIT_PRESSURE=20 !< max number of pressure correction equation iteration
    integer,parameter::MAXIT_DENSITY=20 !< max number of density equation iteration
  