anchor.diff:    double precision,allocatable::flowRho(:) !< flow rate of rho through pairs [kg/s]
     double precision,allocatable::fluxRhou(:,:,:) !< flux of rhou [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 through pairs [J/s]
     double precision,allocatable::advRhoH(:) !< advection rate of rhoH into cell [J/s]
     
     ! state at the beginning of a time step
  
positive.code:     use modNumerics
      use modAdvection
      use modDiffusion
      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:       temp(i)=(rhoE(i)/rho(i)-0.5d0*dot_product(u(:,i),u(:,i)))/(1d0/(gamm-1d0))/Rgas
      end forall
      call setBC()
      call findVarFlow(grid,(rhoE+p)/rho,flowRho,flowRhoH)
      call findAdv(grid,flowRhoH,advRhoH)
      call findDiff(grid,temp,cond,condQ)
      forall(i=1:grid%nC)
        y(i)=rhoE0(i)+dt/grid%v(i)*(advRhoH(i)+condQ(i)) ! TODO add viscous heating
  