anchor.diff:      end if
     end function
     
  -  !> 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
  +  
     
   end module
  
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:     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
  