anchor.diff:        allocate(tmpP(grid%nC))
       end if
       tmpP(:)=p(1:grid%nC)
  -    call setBC()
  -    call findViscForce(grid,u,visc,viscF)
  -    forall(i=1:grid%nE,j=1:DIMS)
  -      fluxRhou(:,j,i)=rhou(j,i)*u(:,i)
  -    end forall
  -    call findAdv(grid,rhou,fluxRhou,flowRhou)
  -    forall(i=1:grid%nC)
  -      rhou1(:,i)=rhou0(:,i)+dt/grid%v(i)*(flowRhou(:,i)+viscF(:,i)) ! pressure force excluded
  -    end forall
       call pressureEq%solve(tmpP,info=info)
       needRetry=info<0
       call pressureEq%getNIt(n)
  
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:       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)
  