anchor.diff:      p(:)=p0(:)
       u(:,:)=u0(:,:)
       temp(:)=temp0(:)
  -    Y(:,:)=Y0(:,:)
  +    massFrac(:,:)=massFrac0(:,:)
     end subroutine
     
     !> derive primitive state {p,u,T} from conserved state {rho,rhou,rhoE}
  -  subroutine deriveState(rhoi,rhoui,rhoEi,Yi,po,uo,tempo)
  +  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)::Yi(:,:) !< mass fraction of species
  +    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
  
positive.code:     use modFileIO
      integer,parameter::FID=10
      double precision::pInit,uInit(DIMS),TInit,pE(DIMS),tmpD
      double precision,allocatable::YInit(:)
      integer::udfIc,udfBc,iUdf(5)
      
      ! read inputs
  
negative.code:     allocate(u0(DIMS,grid%nE))
      allocate(temp(grid%nE))
      allocate(temp0(grid%nE))
      allocate(massFrac(1,grid%nE))! FIXME fix nSp
      allocate(massFrac0(1,grid%nE))! FIXME fix nSp
      allocate(massFracInit(1))! FIXME fix nSp
      allocate(visc(grid%nE))
      allocate(cond(grid%nE))
      allocate(flowRho(grid%nP))
  