anchor.diff:        Y(:,i)=YInit(:)
       end do
       call recoverState(p,u,temp,Y,rho,rhou,rhoE)
  -    call momentumEq%init(grid%nC*DIMS,momentumRHS,maa=10)
       t=0d0
       tNext=tInt
       iOut=0
  
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))
  