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:     allocate(u0(DIMS,grid%nE))
      allocate(temp(grid%nE))
      allocate(temp0(grid%nE))
      allocate(Y(1,grid%nE))! FIXME fix nSp
      allocate(Y0(1,grid%nE))! FIXME fix nSp
      allocate(YInit(1))! FIXME fix nSp
      allocate(visc(grid%nE))
      allocate(cond(grid%nE))
      allocate(flowRho(grid%nP))
  
negative.code:     call grid%clear()
    end subroutine
    
    !> record {rho,rhoU,rhoE,p,u,temp,massFrac} in {rho0,rhoU0,rhoE0,p0,u0,temp0,massFrac0}
    subroutine recordState0()
      rho0(:)=rho(:)
      rhou0(:,:)=rhou(:,:)
  