anchor.diff:      allocate(visc(grid%nE))
       allocate(cond(grid%nE))
       allocate(fluxRhou(DIMS,DIMS,grid%nE))
  -    allocate(fluxRhoH(DIMS,grid%nE))
       allocate(gradP(DIMS,grid%nC))
       allocate(laP(grid%nC))
       allocate(laP1(grid%nC))
  
positive.code:       temp(i)=(rhoE(i)/rho(i)-0.5d0*dot_product(u(:,i),u(:,i)))/(1d0/(gamm-1d0))/Rgas
      end forall
      call setBC()
      forall(i=1:grid%nE)
        fluxRhoH(:,i)=(rhoE(i)+p(i))*u(:,i)
      end forall
      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)*(advRhoH(i)+condQ(i)) ! TODO add viscous heating
  
negative.code:     allocate(flowRho1(grid%nP))
      allocate(advRho(grid%nC))
      allocate(advRhou(DIMS,grid%nC))
      allocate(flowRhoH(grid%nP))
      allocate(advRhoH(grid%nC))
      ! load initial condition
      do i=1,grid%nE
  