C Adjust units of input parametersC      IF (DELTAH.GT.0.) DELTAH_FUEL = DELTAH*1000.      EPUMO2  = EPUMO2*1000.      TMP_LOWER = TMP_LOWER + TMPM      TMP_CRIT = CRITICAL_FLAME_TEMPERATURE + TMPM      DTHRR   = DTSAM      DTMINT  = DTSAM      Y_O2_LL = X_O2_LL*32./(X_O2_LL*32. + (1.-X_O2_LL)*28.)      EC_LL   = VISIBILITY_FACTOR/MAXIMUM_VISIBILITY      E_GAS   = E*1000.CC Set up Mixture Fraction variablesC      IF (MIXTURE_FRACTION) CALL STATE_RELATIONSHIPSCC Set chi_rC      CHI_R = RADIATIVE_FRACTIONCC Get molecular weightsC      DO N=0,NSPEC      CALL GAS_PROPS(SPECIES_ID(N),SIG(N),EPSK(N),MWN(N))      ENDDOCC Compute average molecular weightC      YY0(0)  = 1.      DO N=0,NSPEC      RCON(N) = R0/MWN(N)      ENDDOC      IF (MIXTURE_FRACTION) THENC      RCON_STATE(1) = R0/MW_FUEL      RCON_STATE(2) = R0/MW_O2      RCON_STATE(3) = R0/MW_N2      RCON_STATE(4) = R0/MW_H2O      RCON_STATE(5) = R0/MW_CO2      RCON_STATE(6) = R0/MW_CO      RCON_STATE(7) = R0/MW_H2      RCON_STATE(8) = R0/MW_C      RSUM0         = R0/MW_AVG(0)C      MW_MIN = MW_AVG(0)      MW_MAX = MW_AVG(0)      DO I=1,10000      MW_MIN = MIN(MW_MIN,MW_AVG(I))      MW_MAX = MAX(MW_MAX,MW_AVG(I))      ENDDOC      IF (NSPEC.GT.1) THEN         RCON_MF =  RSUM0         EXTRA_SPECIES_LOOP: DO N=1,NSPEC         IF (N.EQ.IFUEL) CYCLE EXTRA_SPECIES_LOOP         RSUM0 = RSUM0 + YY0(N)*(RCON(N)-RCON_MF)         MW_MIN = MIN(MW_MIN,MWN(N))         MW_MAX = MAX(MW_MAX,MWN(N))         ENDDO EXTRA_SPECIES_LOOP      ENDIFC      RHOA = PINF/(TMPA*RSUM0)C      ENDIFC      IF (.NOT.MIXTURE_FRACTION) THENC      DO N=1,NSPEC      YY0(0) = YY0(0) - YY0(N)      ENDDOC      MW_MIN = MWN(0)      MW_MAX = MWN(0)      RSUM0  = YY0(0)*R0/MWN(0)      DO N=1,NSPEC      MW_MIN = MIN(MW_MIN,MWN(N))      MW_MAX = MAX(MW_MAX,MWN(N))      RSUM0 = RSUM0 + YY0(N)*R0/MWN(N)      ENDDOC      RHOA = PINF/(TMPA*RSUM0)C      ENDIFCC Compute constant-temperature specific heatsC      CP_GAMMA = RCON(0)*GAMMA/(GAMMA-1.)      IF (DT0DZ.EQ.-1000000.) DT0DZ= -GRAV/CP_GAMMA      CPOPR = CP_GAMMA/PRCC Compute variable-temperature specific heats for specific speciesC      IF (.NOT.MIXTURE_FRACTION) THENC      SMOKE3D = .FALSE.C      SPECIES_LOOP: DO N=0,NSPECC      HH(N,0) = 0.C      TEMPERATURE_LOOP: DO J=1,500      TE = MAX(J*0.01_EB,0.301_EB)C      SELECT CASE(SPECIES_ID(N))      CASE DEFAULT   !  Use Molecular Weight      CP(N,J) = 0.001*RCON(N)*GAMMA/(GAMMA-1.)c     CASE('AIR')  ! Nitrogenc     CP(N,J) = 0.931857+0.293529*TE-0.0705765*TE**2+0.00568836*TE**3+c    .          0.001589693/TE**2      CASE('NITROGEN')      CP(N,J) = 0.931857+0.293529*TE-0.0705765*TE**2+0.00568836*TE**3+     .          0.001589693/TE**2      CASE('OXYGEN')      CP(N,J) = 0.926844+0.191789*TE-0.0370788*TE**2+0.00299313*TE**3-     .          0.00686447/TE**2      CASE('METHANE')      IF (J.LE.130) THEN      CP(N,J) = -0.0439393+6.77983*TE-2.6576*TE**2+0.366424*TE**3+     .          0.0424103/TE**2      ELSE      CP(N,J) = 5.36326+0.704042*TE-0.132134*TE**2+0.00863688*TE**3-     .          1.65139/TE**2      ENDIF      CASE('CARBON DIOXIDE')      IF (J.LE.120) THEN      CP(N,J) = 0.568122+1.25425*TE-0.765713*TE**2+0.180645*TE**3-     .          0.00310541/TE**2      ELSE      CP(N,J) = 1.32196+0.0618199*TE-0.0111884*TE**2+0.00082818*TE**3-     .          0.146529/TE**2      ENDIF      CASE('WATER VAPOR')      IF (J.LE.170) THEN      CP(N,J) = 1.67178+0.379583*TE+0.377413*TE**2-0.140804*TE**3+     .          0.00456328/TE**2      ELSE      CP(N,J) = 2.33135+0.479003*TE-0.00833212*TE**2+0.00545106*TE**3-     .          0.619869/TE**2      ENDIF      CASE('HELIUM')      CP(N,J) = 0.001*RCON(N)*GAMMA/(GAMMA-1.)      END SELECT      CP(N,J) = CP(N,J)*1000.            ! J/kg/K      HH(N,J) = HH(N,J-1) + CP(N,J)*10.  ! J/kg      RCP(N,J) = 1./CP(N,J)      ENDDO TEMPERATURE_LOOPC      ENDDO SPECIES_LOOPC      ENDIFCC Compute temperature-dependent specific heat based on mixture fractionC      VARIABLE_CP: IF (MIXTURE_FRACTION) THENC      HH(:,0) = 0.      T_LOOP: DO J=1,500C      TE  = MAX(0.01_EB*J,0.301_EB)C      CP_O2 = 0.926844+0.191789*TE-0.0370788*TE**2+0.00299313*TE**3-     .        0.00686447/TE**2C      CP_N2 = 0.931857+0.293529*TE-0.0705765*TE**2+0.00568836*TE**3+     .        0.001589693/TE**2C      IF (J.LE.130) THEN   ! Ethylene      CP_F  = -0.0439393 + 6.77983*TE - 2.6576*TE**2 + 0.366424*TE**3 +     .         0.0424103/TE**2      ELSE      CP_F  = 5.36326+0.704042*TE-0.132134*TE**2+0.00863688*TE**3-     .        1.65139/TE**2      ENDIFC      IF (J.LE.120) THEN      CP_CO2 = 0.568122+1.25425*TE-0.765713*TE**2+0.180645*TE**3-     .         0.00310541/TE**2      ELSE      CP_CO2 = 1.32196+0.0618199*TE-0.0111884*TE**2+0.00082818*TE**3-     .         0.146529/TE**2      ENDIF      IF (J.LE.170) THEN      CP_H2O = 1.67178+0.379583*TE+0.377413*TE**2-0.140804*TE**3+     .         0.00456328/TE**2      ELSE      CP_H2O = 2.33135+0.479003*TE-0.00833212*TE**2+0.00545106*TE**3-     .         0.619869/TE**2      ENDIFC      Z_LOOP: DO I=0,100      IYY = 100*I      CP(I,J) = ( Y_STATE(IYY,1)*CP_F   +      .            Y_STATE(IYY,2)*CP_O2  +     .            Y_STATE(IYY,3)*CP_N2  +      .            Y_STATE(IYY,4)*CP_H2O +     .            Y_STATE(IYY,5)*CP_CO2 )*1000.      RCP(I,J) = 1./CP(I,J)C      HH(I,J) = HH(I,J-1) + CP(I,J)*10.C      ENDDO Z_LOOPC      ENDDO T_LOOPC      ENDIF VARIABLE_CPCC For mixture fraction model, get the fuel heat of combustionC      IF (MIXTURE_FRACTION) THEN         COMBUSTION_MODEL = 2         DELTAH_FUEL = EPUMO2*MW_O2*NU_O2/(MW_FUEL*NU_FUEL)         DO IPC=1,NPC         LP=>LAGRANGIAN(IPC)         IF (LP%DELTAH.GT.0.)      .   LP%ADJUST_EVAPORATION = LP%DELTAH/DELTAH_FUEL         ENDDO         ENDIFCC For finite rate reaction, set parametersC      DO N=1,NSPEC      IF (FUEL.EQ.SPECIES_ID(N)) IFUEL = N      ENDDOC      IF (IFUEL.GT.0 .AND. IOXYGEN.GT.0) THEN         COMBUSTION_MODEL = 3         MW_FUEL          = MWN(IFUEL)         IF (NUN(IFUEL)  .EQ.0.) NUN(IFUEL)   = -1.         IF (NUN(IFUEL)  .GT.0.) NUN(IFUEL)   = -NUN(IFUEL)         IF (NUN(IOXYGEN).GT.0.) NUN(IOXYGEN) = -NUN(IOXYGEN)         IF (NUN(IOXYGEN).EQ.0.) THEN            WRITE(MESSAGE,'(A)')      .     'ERROR: Specify a stoichiometric coefficient for oxygen'            CALL SHUTDOWN(MESSAGE)            ENDIF         EPUMO2 = DELTAH_FUEL*MWN(IFUEL)*NUN(IFUEL)/     .                       (MWN(IOXYGEN)*NUN(IOXYGEN))      ENDIFCC Compute viscosities and thermal conductivities for species 0-NC      IF (.NOT.MIXTURE_FRACTION) THENC      ALLOCATE(MU_SPEC(0:NSPEC,1:500),STAT=IZERO)      CALL ChkMemErr('READ','MU_SPEC',IZERO)      ALLOCATE(K_SPEC(0:NSPEC,1:500),STAT=IZERO)      CALL ChkMemErr('READ','K_SPEC',IZERO)C      DO N=0,NSPEC      SIGMA2 = SIG(N)**2      DO ITMP=1,500      TSTAR = ITMP*10/EPSK(N)      OMEGA = 1.16145*TSTAR**(-0.14874) +      .        0.52487*EXP(-0.77320*TSTAR) +     .        2.16178*EXP(-2.43787*TSTAR)      MU_N = 26.69E-7*(MWN(N)*ITMP*10)**0.5/(SIGMA2*OMEGA)      IF (MU_USER(N).GE.0.) MU_N = MU_USER(N)*(ITMP*10/TMPA)**0.75      K_N = MU_N*CP(N,ITMP)/PR      IF (K_USER(N).GE.0.)  K_N = K_USER(N)*(ITMP*10/TMPA)**0.75      MU_SPEC(N,ITMP) = MU_N      K_SPEC(N,ITMP)  = K_N      ENDDO      ENDDOC      ALLOCATE(D_SPEC(0:NSPEC,1:500),STAT=IZERO)      CALL ChkMemErr('READ','D_SPEC',IZERO)C      DO N=1,NSPEC      SIGMA2 = (0.5*(SIG(N)+SIG(0)))**2      EPSIJ  = SQRT(EPSK(N)*EPSK(0))      AMW    = SQRT( (MWN(N)+MWN(0))/(2.*MWN(N)*MWN(0)) )      DO ITMP=1,500      TSTAR = ITMP*10/EPSIJ      OMEGA = 1.06036/TSTAR**(0.15610) +     .        0.19300/EXP(0.47635*TSTAR) +     .        1.03587/EXP(1.52996*TSTAR) +     .        1.76474/EXP(3.89411*TSTAR)      D_N0 = 0.00266E-4*AMW*(10*ITMP)**1.5/(SIGMA2*OMEGA)      IF (D_USER(N).GE.0.) D_N0 = D_USER(N)*(ITMP*10/TMPA)**1.75      D_SPEC(N,ITMP) = D_N0      ENDDO      ENDDOC      ENDIFC      IF (MIXTURE_FRACTION) THENC      ALLOCATE(MU_SPEC(0:100,1:500),STAT=IZERO)      CALL ChkMemErr('READ','MU_SPEC',IZERO) ; MU_SPEC = 0.      ALLOCATE(K_SPEC(0:100,1:500),STAT=IZERO)      CALL ChkMemErr('READ','K_SPEC',IZERO)  ; K_SPEC = 0.      ALLOCATE(D_SPEC(0:100,1:500),STAT=IZERO)      CALL ChkMemErr('READ','D_SPEC',IZERO)    ; D_SPEC = 0.C      STATE_SPECIES(1) = 'ETHYLENE'      STATE_SPECIES(2) = 'OXYGEN'      STATE_SPECIES(3) = 'NITROGEN'      STATE_SPECIES(4) = 'WATER VAPOR'      STATE_SPECIES(5) = 'CARBON DIOXIDE'C      DO N=1,5      SIG_N = -1. ; EPSK_N = -1. ; MW_N = -.1      CALL GAS_PROPS(STATE_SPECIES(N),SIG_N,EPSK_N,MW_N)      SIGMA2 = SIG_N**2      DO ITMP=1,500      TSTAR = ITMP*10/EPSK_N      OMEGA = 1.16145*TSTAR**(-0.14874) +     .        0.52487*EXP(-0.77320*TSTAR) +     .        2.16178*EXP(-2.43787*TSTAR)      MU_N = 26.69E-7*(MW_N*ITMP*10)**0.5/(SIGMA2*OMEGA)      DO IYY=0,100      K_N  = MU_N*CP(IYY,ITMP)/PR      WGT  = Y_STATE(100*IYY,N)*MW_AVG(100*IYY)/MW_N      MU_SPEC(IYY,ITMP) = MU_SPEC(IYY,ITMP) + WGT*MU_N      K_SPEC(IYY,ITMP)  = K_SPEC(IYY,ITMP)  + WGT*K_N      ENDDO      ENDDOC      SIGMA2 = (0.5*(SIG_N+SIG(0)))**2      EPSIJ  = SQRT(EPSK_N*EPSK(0))      AMW    = SQRT( (MW_N+MWN(0))/(2.*MW_N*MWN(0)) )      DO ITMP=1,500      TSTAR = ITMP*10/EPSIJ      OMEGA = 1.06036/TSTAR**(0.15610) +     .        0.19300/EXP(0.47635*TSTAR) +     .        1.03587/EXP(1.52996*TSTAR) +     .        1.76474/EXP(3.89411*TSTAR)      D_N0 = 0.00266E-4*AMW*(10*ITMP)**1.5/(SIGMA2*OMEGA)      DO IYY=0,100      D_SPEC(IYY,ITMP) = D_SPEC(IYY,ITMP) + Y_STATE(100*IYY,N)*D_N0      ENDDO      ENDDO      ENDDOC      ENDIFCC Assign Mass Per Unit EnergyC      IF (IOXYGEN.GT.0) THEN      ALLOCATE(MPUE(NSPEC))      ALLOCATE(EPUM(NSPEC))      CALL ChkMemErr('READ','MPUE',IZERO)      MPUE(1:NSPEC) = NUN(1:NSPEC)*MWN(1:NSPEC)/     .          ABS(EPUMO2*MWN(IOXYGEN)*NUN(IOXYGEN))      EPUM = 1./MPUE      ENDIFCC Define all possible output quantitiesC      CALL DEFINE_OUTPUT_QUANTITIESC      CONTAINSC C      SUBROUTINE GAS_PROPS(SPECIES,SIGMA,EPSOK,MW)CC Molecular weight (g/mol) and Lennard-Jones propertiesC      REAL(EB) SIGMA,EPSOK,MW,SIGMAIN,EPSOKIN,MWIN      CHARACTER(26) SPECIESC      SIGMAIN = SIGMA      EPSOKIN = EPSOK      MWIN    = MWC      SELECT CASE(SPECIES)      CASE('AIR')             ; SIGMA=3.711 ; EPSOK= 78.6  ; MW=29.      CASE('CARBON MONOXIDE') ; SIGMA=3.690 ; EPSOK= 91.7  ; MW=28.      CASE('CARBON DIOXIDE')  ; SIGMA=3.941 ; EPSOK=195.2  ; MW=44.      CASE('ETHYLENE')        ; SIGMA=4.163 ; EPSOK=224.7  ; MW=28.      CASE('HELIUM')          ; SIGMA=2.551 ; EPSOK= 10.22 ; MW= 4.      CASE('HYDROGEN')        ; SIGMA=2.827 ; EPSOK= 59.7  ; MW= 2.      CASE('METHANE')         ; SIGMA=3.758 ; EPSOK=148.6  ; MW=16.      CASE('NITROGEN')        ; SIGMA=3.798 ; EPSOK= 71.4  ; MW=28.      CASE('OXYGEN')          ; SIGMA=3.467 ; EPSOK=106.7  ; MW=32.      CASE('PROPANE')         ; SIGMA=5.118 ; EPSOK=237.1  ; MW=44.      CASE('WATER VAPOR')     ; SIGMA=2.641 ; EPSOK=809.1  ; MW=18.      CASE DEFAULT            ; SIGMA=3.711 ; EPSOK= 78.6  ; MW=29.      END SELECT  C      IF (SIGMAIN.GT.0.) SIGMA = SIGMAIN      IF (EPSOKIN.GT.0.) EPSOK = EPSOKIN      IF (MWIN   .GT.0.) MW    = MWIN C      IF (SPECIES.EQ.'MIXTURE_FRACTION') MW = MW_FUELC      END SUBROUTINE GAS_PROPSCC      SUBROUTINE STATE_RELATIONSHIPSCC Calculate state relations, Y_STATE(Z), and average molecular weight.C Y_STATE(0:10000,I): I=1-Fuel,2-O2,3-N2,4-H2O,5-CO2,6-CO,7-H2,8-SOOTC      REAL(EB) Z_TERM,S_TERM,TOTAL_MASS,NU_N2_AIR,NU_N2_DILUENT,ETA,CHI,     .         Y_N2_INLETCC Normalize ideal stoichiometric coefficientsC      IF (NU_FUEL.NE.1.) THEN         NU_O2   = NU_O2/NU_FUEL         NU_H2O  = NU_H2O/NU_FUEL         NU_CO2  = NU_CO2/NU_FUEL         NU_FUEL = 1.         ENDIFC      Y_N2_INFTY = 1.-Y_O2_INFTYCC If nitrogen is in the fuel, correct MW_FUEL and Y_F_INLETC      Y_F_INLET  = Y_F_INLET*(1.-MW_N2*FUEL_N2/MW_FUEL)      MW_FUEL    = MW_FUEL - MW_N2*FUEL_N2      Y_N2_INLET = 1.-Y_F_INLET CC Apply Faeth's correlation to get CO yield from soot yieldC      IF (CO_YIELD.LT.0.) THEN      GEN_TO_YLD = NU_CO2*MW_C/MW_FUEL      CO_YIELD   = GEN_TO_YLD*(0.0014 + 0.37*SOOT_YIELD/GEN_TO_YLD)      ENDIFCC Adjust ideal stoichiometric coefficients to account for soot, CO, etc.C      NU_SOOT = SOOT_YIELD*MW_FUEL/MW_C      NU_H2   = H2_YIELD*MW_FUEL/MW_H2