anchor.diff:      allocate(YInit(1))
       allocate(visc(grid%nE))
       allocate(cond(grid%nE))
  -    allocate(fluxRho(DIMS,grid%nE))
       allocate(fluxRhou(DIMS,DIMS,grid%nE))
       allocate(fluxRhoH(DIMS,grid%nE))
       allocate(gradP(DIMS,grid%nC))
  
positive.code:     call findPresForce(grid,p,presF)
      forall(i=1:grid%nC)
        rhou(:,i)=rhou(:,i)+dt/grid%v*(presF(:,i)-presF1(:,i))
        !u(:,i)=rhou(:,i)/rho(i)
      end forall
    end subroutine
    
  
negative.code:       allocate(tmpP(grid%nC))
      end if
      tmpP(:)=p(1:grid%nC)
      call findAdv(grid,rho,rhou,flowRho) ! vary only the Rhie-Chow part of flowRho in RHS
      call pressureEq%solve(tmpP)
      p(1:grid%nC)=tmpP(:)
    end subroutine
  