anchor.diff:          p(n)=pGst
           u(:,n)=uGst(:)
           temp(n)=TGst
  -        rho(n)=pGst/Rgas/TGst
  -        rhou(:,n)=rho(n)*uGst(:)
  -        rhoE(n)=rho(n)*(1d0/(gamm-1d0)*Rgas*TGst+0.5d0*dot_product(uGst,uGst))
  +        rho(n)=rho(m)!pGst/Rgas/TGst
  +        rhou(:,n)=-rhou(:,m)!rho(n)*uGst(:)
  +        rhoE(n)=rhoE(m)!rho(n)*(1d0/(gamm-1d0)*Rgas*TGst+0.5d0*dot_product(uGst,uGst))
         end if
       end do
     end subroutine
  
positive.code:     end if
      
      ! 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)/Rgas)
      rhoScale=maxval(rho)
      rhouScale=rhoScale*uScale
      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)
      call pressureEq%setTol(pScale*RTOL_PRESSURE)
      call densityEq%setTol(rhoScale*RTOL_DENSITY)
      call energyEq%setTol(rhoEScale*RTOL_ENERGY)
  
negative.code:   double precision,parameter::RTOL_MOMENTUM=1d-8 !< momentum prediction relative tolerance
    double precision,parameter::RTOL_PRESSURE=1d-8 !< pressure correction relative tolerance
    double precision,parameter::RTOL_DENSITY=1d-7 !< density equation relative tolerance
    double precision,parameter::RTOL_ENERGY=1d-7 !< energy equation relative tolerance
    
    integer,parameter::BC_WALL_TEMP=0 !< wall boundary with prescribed temperature
    integer,parameter::BC_WALL_TEMP_UDF=5 !< wall boundary with prescribed temperature by UDF
  