anchor.diff:      end if
     end function
     
  +  !> 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 energy equation in the form of a fix point problem
  +  function energyRHS(oldVector,newVector,dat)
  +    use iso_c_binding
  +    use ieee_arithmetic
  +    use modNumerics
  +    use modAdvection
  +    use modDiffusion
  +    type(C_PTR),value::oldVector !< old N_Vector
  +    type(C_PTR),value::newVector !< new N_Vector
  +    type(C_PTR),value::dat !< optional user data object
  +    integer(C_INT)::energyRHS !< error code
  +    double precision,pointer::x(:) !< fortran pointer associated with oldVector
  +    double precision,pointer::y(:) !< fortran pointer associated with newVector
  +    
  +    call associateVector(oldVector,x)
  +    call associateVector(newVector,y)
  +    
  +    rhoE(1:grid%nC)=x(1:grid%nC)
  +    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
  +    call setBC()
  +    call findVarFlow(grid,(rhoE+p)/rho,flowRho,flowRhoH)
  +    call findAdv(grid,flowRhoH,advRhoH)
  +    call findDiff(grid,temp,cond,condQ)
  +    forall(i=1:grid%nC)
  +      y(i)=rhoE0(i)+dt/grid%v(i)*(advRhoH(i)+condQ(i)) ! TODO add viscous heating
  +    end forall
  +    energyRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
  +    if(c_associated(dat))then
  +    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
  
positive.code:     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
  
negative.code:     if(c_associated(dat))then
      end if
    end function
    
    !> 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
    
    !> residual function of the pressure equation
    function pressureRHS(oldVector,newVector,dat)
  