anchor.diff:      end forall
     end subroutine
     
  -  !> solve the pressure, such that mass is conserved
  -  subroutine solvePressure()
  -    use modAdvection
  -    double precision,allocatable,save::tmpP(:)
  -    integer::info
  -    
  -    if(.not.allocated(tmpP))then
  -      allocate(tmpP(grid%nC))
  -    else if(size(tmpP)/=grid%nC)then
  -      deallocate(tmpP)
  -      allocate(tmpP(grid%nC))
  -    end if
  -    tmpP(:)=p(1:grid%nC)
  -    call setBC()
  -    call findMassFlow(grid,rhou,flowRho1)
  -    call pressureEq%solve(tmpP,info=info)
  -    needRetry=info<0
  -    call pressureEq%getNIt(n)
  -    nItPressure=max(nItPressure,n)
  -    p(1:grid%nC)=tmpP(:)
  -  end subroutine
  -  
  -  !> correct the momentum, mass flow, mass advection with updated pressure
  -  subroutine correctMomentum()
  -    use modPressure
  -    use modGradient
  -    use modAdvection
  -    use modRhieChow
  -    integer,parameter::N_RHOU_CORRECTION=3 !< number of enabled momentum corrections
  -    
  -    if(nItPISO<=N_RHOU_CORRECTION)then
  -      call setBC()
  -      call findPresForce(grid,p,gradP,presF)
  -      forall(i=1:grid%nC)
  -        rhou(:,i)=rhou(:,i)+dt/grid%v(i)*(presF(:,i)-presF1(:,i))
  -      end forall
  -    end if
  -    call setBC()
  -    call findMassFlow(grid,rhou,flowRho)
  -    call findGrad(grid,p,gradP)
  -    call addRhieChow(grid,p,gradP,dt,flowRho)
  -    call findAdv(grid,flowRho,advRho)
  -  end subroutine
  -  
  -  !> solve the density
  -  subroutine solveDensity()
  -    double precision,allocatable,save::tmpRho(:)
  -    integer::info
  -    
  -    if(.not.allocated(tmpRho))then
  -      allocate(tmpRho(grid%nC))
  -    else if(size(tmpRho)/=grid%nC)then
  -      deallocate(tmpRho)
  -      allocate(tmpRho(grid%nC))
  -    end if
  -    tmpRho(:)=rho(1:grid%nC)
  -    call densityEq%solve(tmpRho,info)
  -    needRetry=info<0
  -    call densityEq%getNIt(n)
  -    nItDensity=max(nItDensity,n)
  -    rho(1:grid%nC)=tmpRho(:)
  -  end subroutine
  -  
  -  !> solve the (total) energy, during which T is updated
  -  subroutine solveEnergy()
  -    use modAdvection
  -    use modDiffusion
  -    double precision,allocatable,save::tmpRhoE(:)
  -    integer::info
  -    
  -    if(.not.allocated(tmpRhoE))then
  -      allocate(tmpRhoE(grid%nC))
  -    else if(size(tmpRhoE)/=grid%nC)then
  -      deallocate(tmpRhoE)
  -      allocate(tmpRhoE(grid%nC))
  -    end if
  -    tmpRhoE(:)=rhoE(1:grid%nC)
  -    call energyEq%solve(tmpRhoE,info)
  -    needRetry=info<0
  -    call energyEq%getNIt(n)
  -    nItEnergy=max(nItEnergy,n)
  -    rhoE(1:grid%nC)=tmpRhoE(:)
  -    forall(i=1:grid%nC)
  -      temp(i)=(rhoE(i)/rho(i)-0.5d0*dot_product(u(:,i),u(:,i)))/(1d0/(gamm-1d0))/Rgas
  -    end forall
  -  end subroutine
  -  
     !> RHS of the momentum equation in the form of a fix point problem
     function momentumRHS(oldVector,newVector,dat)
       use iso_c_binding
  
positive.code:     end if
    end function
    
    !> RHS of the density equation in the form of a fix point problem
    function densityRHS(oldVector,newVector,dat)
      use iso_c_binding
  
negative.code:     end if
    end function
    
    
    
  end module
  