anchor.diff:    double precision,allocatable::presF(:,:) !< pressure force applied on cell [N]
     double precision,allocatable::condQ(:) !< heat conduction into cell [W]
     double precision,allocatable::carrier(:,:) !< {rho,rhou,rhoH} of each cell [*]
  -  double precision,allocatable::flowRho(:) !< flow rate of rho [kg/s]
  +  double precision,allocatable::advRho(:) !< advection rate of rho [kg/s]
     double precision,allocatable::fluxRhou(:,:,:) !< flux of rhou [kg/m/s^2]
  -  double precision,allocatable::flowRhou(:,:) !< flow rate of rhou into cell [kg*m/s^2]
  +  double precision,allocatable::advRhou(:,:) !< advection rate of rhou into cell [kg*m/s^2]
     double precision,allocatable::fluxRhoH(:,:) !< flux of rhoH [J/m^2/s]
  -  double precision,allocatable::flowRhoH(:) !< flow rate of rhoH into cell [J/s]
  +  double precision,allocatable::advRhoH(:) !< advection rate of rhoH into cell [J/s]
     
     ! state at the beginning of a time step
     double precision,allocatable::rho0(:),rhou0(:,:),rhoE0(:),p0(:),u0(:,:),temp0(:),Y0(:,:)
     
     ! auxiliary variables for the PISO iteration
  -  double precision,allocatable::rho1(:),p1(:),laP1(:),presF1(:,:),temp1(:),rhou1(:,:),flowRho1(:)
  +  double precision,allocatable::rho1(:),p1(:),laP1(:),presF1(:,:),temp1(:),rhou1(:,:),advRho1(:)
     
     ! scales
     double precision,allocatable::rhoScale !< density scale [kg/m^3]
  
positive.code:     forall(i=1:grid%nE)
        fluxRhoH(:,i)=(rhoE(i)+p(i))*u(:,i)
      end forall
      call findAdv(grid,rhoE+p,fluxRhoH,flowRhoH)
      call addRhieChow(grid,rhoE+p,p,gradP,rho,dt,flowRhoH)
      call findDiff(grid,temp,cond,condQ)
      forall(i=1:grid%nC)
        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
  
negative.code:     integer(C_INT)::pressureRHS !< error code
      double precision,pointer::x(:) !< fortran pointer associated with oldVector
      double precision,pointer::y(:) !< fortran pointer associated with newVector
      double precision,allocatable,save::tmpFlowRho(:),tmpAdvRho(:)
      
      call associateVector(oldVector,x)
      call associateVector(newVector,y)
  