anchor.diff:      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 findAdv(grid,rhoE+p,fluxRhoH,advRhoH)
  +    call addRhieChow(grid,rhoE+p,p,gradP,rho,dt,advRhoH)
       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
  +      y(i)=rhoE0(i)+dt/grid%v(i)*(advRhoH(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
  
positive.code:     call setBC()
      call findGrad(grid,p,gradP)
      if(.not.allocated(tmpFlowRho))then
        allocate(tmpFlowRho(grid%nC))
      else if(size(tmpFlowRho)/=grid%nC)then
        deallocate(tmpFlowRho)
        allocate(tmpFlowRho(grid%nC))
      end if
      tmpFlowRho(:)=flowRho(:)
      call addRhieChow(grid,rho,p,gradP,rho,dt,tmpFlowRho)
      call findDiff(grid,p,[(1d0,i=1,grid%nC)],laP)
      forall(i=1:grid%nC)
        y(i)=p(i)-Rgas*temp(i)*(rho0(i)+dt/grid%v(i)*(tmpFlowRho(i)+dt*(laP(i)-laP1(i))))
      end forall
      pressureRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
      if(c_associated(dat))then
  
negative.code:     forall(i=1:grid%nE,j=1:DIMS)
        fluxRhou(:,j,i)=rhou(j,i)*u(:,i)
      end forall
      call findAdv(grid,rhou,fluxRhou,advRhou)
      forall(i=1:grid%nC)
        rhou(:,i)=rhou0(:,i)+dt/grid%v(i)*(advRhou(:,i)+presF(:,i)+viscF(:,i))
      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))))
  