anchor.diff:          uInit(3)=udf%eval(iUdf(5),pE,0d0)
         end if
         ! FIXME remove below
  -      YInit=[1d0]
  +      massFracInit=[1d0]
         ! FIXME remove above
         p(i)=pInit
         u(:,i)=uInit(:)
         temp(i)=TInit
  -      Y(:,i)=YInit(:)
  +      massFrac(:,i)=massFracInit(:)
       end do
  -    call recoverState(p,u,temp,Y,rho,rhou,rhoE)
  +    call recoverState(p,u,temp,massFrac,rho,rhou,rhoE)
       t=0d0
       tNext=tInt
       iOut=0
  -    deallocate(YInit)
  +    deallocate(massFracInit)
     end subroutine
     
     !> clear the simulation environment
  
positive.code:     p0(:)=p(:)
      u0(:,:)=u(:,:)
      temp0(:)=temp(:)
      Y0(:,:)=Y(:,:)
    end subroutine
    
    !> load {rho,rhoU,rhoE,p,u,temp,Y} from {rho0,rhoU0,rhoE0,p0,u0,temp0,Y0}
    subroutine loadState0()
      rho(:)=rho0(:)
      rhou(:,:)=rhou0(:,:)
  
negative.code:     p(:)=p0(:)
      u(:,:)=u0(:,:)
      temp(:)=temp0(:)
      massFrac(:,:)=massFrac0(:,:)
    end subroutine
    
    !> derive primitive state {p,u,T} from conserved state {rho,rhou,rhoE}
    subroutine deriveState(rhoi,rhoui,rhoEi,massFraci,po,uo,tempo)
      double precision,intent(in)::rhoi(:) !< density
      double precision,intent(in)::rhoui(:,:) !< momentum
      double precision,intent(in)::rhoEi(:) !< total energy
      double precision,intent(in)::massFraci(:,:) !< mass fraction of species
      double precision,intent(inout)::po(:) !< pressure
      double precision,intent(inout)::uo(:,:) !< velocity
      double precision,intent(inout)::tempo(:) !< temperature
  