anchor.diff:      call associateVector(newVector,y)
       
       rhou(:,1:grid%nC)=reshape(x,[DIMS,grid%nC])
  +    call setBC()
       call findViscForce(grid,u,visc,viscF)
       forall(i=1:grid%nE,j=1:DIMS)
         fluxRhou(:,j,i)=rhou(j,i)*u(:,i)
  
positive.code:   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
    
  contains
    
  
negative.code:       Y(:,i)=YInit(:)
      end do
      call recoverState(p,u,temp,Y,rho,rhou,rhoE)
      t=0d0
      tNext=tInt
      iOut=0
  