anchor.diff:        f2c(:,3,i)=rhou(3,i)*u(:,i)+[0d0,0d0,p(i)]
         f3c(:,i)=(rhoE(i)+p(i))*u(:,i)
       end forall
  +    !$omp end workshare
       ! Roe flux difference splitting
  +    !$omp parallel do default(shared)&
  +    !$omp& private(m,n,rhoAvg,uAvg,Havg,cAvg,uNormAvg,rhoJump,pJump,uJump,uNormJump,&
  +    !$omp&         dFEntropy,dFAcoustic1,dFAcoustic2,flow)
       do i=1,grid%nP
         m=grid%iEP(1,i)
         n=grid%iEP(2,i)
  
positive.code:         &         dot_product(0.5d0*(f2c(:,3,m)+f2c(:,3,n)),grid%normP(:,i)),&
          &         dot_product(0.5d0*(f3c(:,m)+f3c(:,n)),grid%normP(:,i))]&
          &        -0.5d0*(dFEntropy(:)+dFAcoustic1(:)+dFAcoustic2(:)))
          if(n<=grid%nC)then
            dRho(m)=dRho(m)-flow(1)
            dRho(n)=dRho(n)+flow(1)
  
negative.code:     allocate(f1c(DIMS,k))
      allocate(f2c(DIMS,DIMS,k))
      allocate(f3c(DIMS,k))
      !$omp workshare
      forall(i=1:k)
        u(:,i)=rhou(:,i)/rho(i)
        H(i)=(rhoE(i)+p(i))/rho(i)
  