anchor.diff:        deallocate(tmpRhou)
         allocate(tmpRhou(DIMS*grid%nC))
       end if
  +    ! solve momentum equation with current pressure
       tmpRhou=reshape(rhou(:,1:grid%nC),[DIMS*grid%nC])
       call setBC()
       call findGrad(grid,p,gradP)
  
positive.code:     call momentumEq%getNIt(n)
      nItMomentum=max(nItMomentum,n)
      rhou(:,1:grid%nC)=reshape(tmpRhou,[DIMS,grid%nC])
    end subroutine
    
    !> solve the pressure, such that mass is conserved
    subroutine solvePressure()
      use modAdvection
      use modNewtonian
      double precision,allocatable,save::tmpP(:)
      integer::info
      
  
negative.code:     end do
    end subroutine
    
    !> predict the rhou1 and flowRho without acceleration by pressure
    subroutine predictMomentum()
      use modGradient
      use modPressure
  