anchor.diff:      call associateVector(newVector,y)
       
       rhoE(1:grid%nC)=x(1:grid%nC)
  -    forall(i=1:grid%nE)
  +    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()
  
positive.code:     call energyEq%getNIt(n)
      nItEnergy=max(nItEnergy,n)
      rhoE(1:grid%nC)=tmpRhoE(:)
    end subroutine
    
    !> RHS of the momentum equation in the form of a fix point problem
  
negative.code:   
    !> set the boundary conditions
    subroutine setBC()
      double precision::pGst,TGst,uGst(DIMS),rhoGst,rhouGst(DIMS),rhoEGst,Mach,pP(DIMS),Tt,pt
      
      do i=1,grid%nP
        m=grid%iEP(1,i)
  