anchor.diff:    double precision,allocatable::uScale !< velocity scale [m/s]
     double precision,allocatable::tempScale !< temperature scale [K]
     
  -  type(fixPtEq)::momentumEq !< momentum equation as a fix point problem
  -  
  -  !double precision,allocatable::x(:) !< solution vector of the nonlinear system of equations
  -  !double precision,allocatable::xscale(:) !< scaling factors for the solution
  -  !double precision,allocatable::rscale(:) !< scaling factors for the residual
  -  !integer(C_LONG)::nEq !< number of equations
  -  !integer(C_LONG)::ioutFKIN(100) !< integer output of FKINSOL
  -  !double precision::routFKIN(100) !< real output of FKINSOL
  +  ! data for algebraic solver
  +  type(fixPtEq)::momentumEq !< momentum equation as a fix point problem  
     
   contains
     
  
positive.code:     call associateVector(newVector,y)
      
      rhou(:,1:grid%nC)=reshape(x,[DIMS,grid%nC])
      call findViscForce(grid,u,visc,viscF)
      forall(i=1:grid%nE,j=1:DIMS)
        fluxRhou(:,j,i)=rhou(j,i)*u(:,i)
  
negative.code:     !    end if
      !  end do
      !end do
      ! initialize algebraic solver
      call momentumEq%init(grid%nC*DIMS,momentumRHS,maa=10)
      ! work space and initial state
      allocate(rho(grid%nE))
      allocate(rho0(grid%nE))
  