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 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)-R_AIR*temp(i)*dt**2/grid%v(i)*(laP(i)-laP1(i))-R_AIR*temp(i)/R_AIR/temp1(i)*p1(i)& ! FIXME derive the pressure correction eq!
        &    +R_AIR*temp(i)*(rho1(i)-rho0(i))-R_AIR*temp(i)*dt/grid%v(i)*tmpFlowRho(i)!tmpFlowRho(i)
      end forall
      write(*,*)'pressure',maxval(abs(y(1:grid%nC)))
      pressureRHS=0
  
negative.code:     call findPresForce(grid,p,presF)
      forall(i=1:grid%nC)
        rhou(:,i)=rhou(:,i)+dt/grid%v*(presF(:,i)-presF1(:,i))
      end forall
    end subroutine
    
  