anchor.diff:      
       p(1:grid%nC)=x(1:grid%nC)
       call setBC()
  +    call findGrad(grid,p,gradP)
       if(.not.allocated(tmpFlowRho))then
         allocate(tmpFlowRho(grid%nC))
       else if(size(tmpFlowRho)/=grid%nC)then
  
positive.code:       y(i)=p(i)-R_AIR*temp(i)*dt**2/grid%v(i)*(laP(i)-laP1(i))-R_AIR*temp(i)/R_AIR/temp1(i)*p1(i)&
        &    +R_AIR*temp(i)*(rho1(i)-rho0(i))-R_AIR*temp(i)*dt/grid%v(i)*tmpFlowRho(i)
      end forall
      pressureRHS=0
    end function
    
    !> RHS of the density equation in the form of a fix point problem
    function densityRHS(oldVector,newVector,dat)
      use iso_c_binding
      use modNumerics
      use modAdvection
      use modRhieChow
  
negative.code:   !> RHS of the momentum equation in the form of a fix point problem
    function momentumRHS(oldVector,newVector,dat)
      use iso_c_binding
      use ieee_arithmetic
      use modNumerics
      use modAdvection
      use modNewtonian
      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
  