TEXT(CH4 Statoil Furnace
TITLE
  DISPLAY
  The case considered is the turbulent combustion of CH4 with
  O2 in an industrial furnace. A swirling central jet of
  'oxidiser'enters the furnace at 552K with the composition of
  80% O2 and 20% H2O, and mixes and burns with an annular 'fuel'
  jet of 25% CH4, 12%CO, 40% H2O & 23% CO2 at a temperature of
  1055 K.
 
  The fast-chemistry assumption is invoked for main combustion
  calculation in terms of two competing single-step infinite-rate
  reactions: 2H2 + O2 > 2H2O and CO + O2 > CO2. The kinetic
  calculation is peformed as a post-processing run with a frozen
  hydrodynamic and temperature fields.
  ENDDIS
 IRUNN   =       1 ;LIBREF =       0
 CHAR(ANS2,CMOD);BOOLEAN(THRAD,KINET,TWOPNT)
 REAL(PI,RV,RQ,RF,Z1,Z2,Z3,NCZ1,NCZ2,NCZ3,NCY1,NCY2,NCZALL)
 REAL(ALPHA,SCALE,KEIN,KESWIRL,EPIN,EPSWIRL)
 REAL(MAXV,MINCELL,RELX); RELX=1.0E-3
 PI     = 3.1415
 * Radius vortex generator
 RV     = 0.095
 * Radius quarl
 RQ     = 0.190
 * Radius furnace
 RF     = 0.220
 * Axial distance from vortex generator to quarl
 Z1     = 0.142
 * Axial distance over quarl
 Z2     = 0.261
 * Axial distance over furnace
 Z3     = 2.097
 * Number of cells in radial and axial dir.
 NCY1=9;NCY2=2
 NCZ1=4;NCZ2=8;NCZ3=15
 NCZALL=NCZ1+NCZ2+NCZ3
 * U0 er aksiell hastighet ved innloepet
 * Re = U0*RV/ENUL = 4.7*0.095/1.78E-06 = 250 000
  U0     = 4.7
  * WMAX er tagentiell hastighet ved innloept, r=Rmax
  * WMAX = S0*2*U0 hvor S0=0.7 gir WMAX = 6.6
  WMAX   = 6.6
 ALPHA = 55.0*3.1416/180.0
 REAL(WIN,WSWIRL,USWIRL)
 WIN = 7.0;WSWIRL = 3.0;USWIRL = WSWIRL*TAN(ALPHA)
 
 KEIN = 5*(WIN**2)/100
 EPIN = 0.09**0.75*KEIN**1.5/0.007*RV/2
 KESWIRL = 5*(WSWIRL**2)/100
 EPSWIRL = 0.09**0.75*KESWIRL**1.5/0.007*RV/2
 
  *MIB
 REAL(GASCON,TFUEL,TOX,RHOOX,RHOFU,TWAL)
 REAL(YCH4F,YH2F,YO2F,YH2OF,YCOF,YCO2F,YN2F)
 REAL(YCH4O,YH2O,YO2O,YH2OO,YCOO,YCO2O,YN2O)
 GASCON=8314.0;PRESS0=38.E5;TWAL=800.
   ** Fuel and oxidiser stream inlet values
 TOX=552;TFUEL=1055
 YCH4F=0.25;YH2F=0.0;YO2F=0.0;YH2OF=0.4;YCOF=0.12
 YCO2F=0.23;YN2F=0.
 YCH4O=0.0;YH2O=0.0;YO2O=0.8;YH2OO=0.2;YCOO=0.0
 YCO2O=0.0;YN2O=0.
MESG( Do you want a KINETics run? (default=N)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ KINET=T
+ MESG( Kinetics run selected.
ELSE
+ KINET=F
+ MESG( Main combustion run selected.
ENDIF
  Group 1. Run Title
  Group 2. Transience
 STEADY  =    T
  Groups 3, 4, 5  Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
 RSET(M,1,NCY1+NCY2,NCZ1+NCZ2+NCZ3)
    * Set overall domain extent:
    *        xulast  yvlast  zwlast    name
 XSI= 1.000E+00;YSI= RF;ZSI= Z1+Z2+Z3;
 RSET(D,CHAM,XSI,YSI,ZSI)
  Group 6. Body-Fitted coordinates
 BFC=T
    * Set points
 XPO= 0.0000E+00;YPO= 0.0000E+00;ZPO= 0.0000E+00;
 GSET(P,P1,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RV        ;ZPO= 0.0000E+00;
 GSET(P,P2,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RF        ;ZPO= 0.0000E+00;
 GSET(P,P3,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= 0.0000E+00;ZPO= Z1        ;
 GSET(P,P4,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RV        ;ZPO= Z1        ;
 GSET(P,P5,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RF        ;ZPO= Z1        ;
 GSET(P,P6,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= 0.0000E+00;ZPO= Z1+Z2     ;
 GSET(P,P7,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RQ        ;ZPO= Z1+Z2     ;
 GSET(P,P8,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RF        ;ZPO= Z1+Z2     ;
 GSET(P,P9,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= 0.0000E+00;ZPO= Z1+Z2+Z3  ;
 GSET(P,P10,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RQ        ;ZPO= Z1+Z2+Z3  ;
 GSET(P,P11,XPO,YPO,ZPO)
 XPO= 0.0000E+00;YPO= RF        ;ZPO= Z1+Z2+Z3  ;
 GSET(P,P12,XPO,YPO,ZPO)
    * Set lines/arcs
 GSET(L,L1,P1,P2,NCY1,1.0)
 GSET(L,L2,P2,P3,NCY2,1.0)
 GSET(L,L3,P4,P5,NCY1,1.0)
 GSET(L,L4,P5,P6,NCY2,1.0)
 GSET(L,L5,P7,P8,NCY1,1.0)
 GSET(L,L6,P8,P9,NCY2,1.0)
 GSET(L,L7,P10,P11,NCY1,1.0)
 GSET(L,L8,P11,P12,NCY2,1.0)
 
 GSET(L,L9,P1,P4,NCZ1,1.0)
 GSET(L,L10,P4,P7,NCZ2,1.0)
 GSET(L,L11,P7,P10,NCZ3,1.2)
 GSET(L,L12,P2,P5,NCZ1,1.0)
 GSET(L,L13,P5,P8,NCZ2,1.0)
 GSET(L,L14,P8,P11,NCZ3,1.2)
 GSET(L,L15,P3,P6,NCZ1,1.0)
 GSET(L,L16,P6,P9,NCZ2,1.0)
 GSET(L,L17,P9,P12,NCZ3,1.2)
    * Set frames
 GSET(F,F1,P1,-,P2,-,P5,-,P4,-)
 GSET(F,F2,P2,-,P3,-,P6,-,P5,-)
 GSET(F,F3,P4,-,P5,-,P8,-,P7,-)
 GSET(F,F4,P5,-,P6,-,P9,-,P8,-)
 GSET(F,F5,P7,-,P8,-,P11,-,P10,-)
 GSET(F,F6,P8,-,P9,-,P12,-,P11,-)
    * Match a grid mesh
 GSET(M,F1,+J+K,1,1,1,TRANS)
 GSET(M,F2,+J+K,1,NCY1+1,1,TRANS)
 GSET(M,F3,+J+K,1,1,NCZ1+1,TRANS)
 GSET(M,F4,+J+K,1,NCY1+1,NCZ1+1,TRANS)
 GSET(M,F5,+J+K,1,1,NCZ1+NCZ2+1,TRANS)
 GSET(M,F6,+J+K,1,NCY1+1,NCZ1+NCZ2+1,TRANS)
    * Copy/Transfer/Block grid planes
    * rotasjon
 GSET(C,I2,F,I1,1,NCY1+NCY2,1,NCZALL,RZ,-1.0000E-01,0,0,INC,1)
    * translasjon
    GSET(C,I4,F,I1,1,NCY1+NCY2,1,NCZALL,+,1.0000E-01,0,0,INC,1)
NONORT  =    T
  Group 7. Variables: STOREd,SOLVEd,NAMEd
IF(KINET) THEN
+ STORE(P1,U1,V1,W1,KE,EP)
ELSE
+ SOLVE(P1,U1,V1,W1);SOLUTN(P1,P,P,Y,P,P,P)
+ TURMOD(KECHEN)
ENDIF
STORE(ENUT,LEN1);STORE(H1,LEN1,RHO1,TMP1)
  Group 9. Properties
ENUL    = 1.780E-06
IF(KINET) THEN
+ CHEMKIN(SPECIES,CH4,O2,H2,CO,H2O,CO2,CH3,CH2,CH,CH2O,HCO,H,O,OH,H$
O2,H2O2,N2)
+ STORE(F,SPH1,MMWT,SRAD)
ELSE
  *** START OF EXTENDED SCRS MODEL SETTINGS
   CFUEL=CH4 COXID=O2 CFP1=CO2 CFP2=H2O NRSYS=-2
INTEGER(NSPEC,NELEM);NSPEC=7;NELEM=4
INTEGER(ICOMB,NCSTEP,NCREAC)
MESG( Enter required combustion model
MESG( Default: 1 step finite rate EBU
MESG( The options are:
MESG(  EBU1  - 1 step finite-rate EBU
MESG(  EBU2  - 2 step 2 reaction finite-rate EBU
MESG(  EBU3  - 2 step 3 reaction finite-rate EBU
MESG(  BURN  - Infinite-rate model
MESG(  DDEL  - Infinite-rate model with Double-Delta PDF
READVDU(CMOD,CHAR,EBU1)
CASE :CMOD: OF
WHEN EBU1,4
+ MESG(1 step finite-rate EBU model
+ MESG(2CH4 +  4O2 > 2CO2 + 4H2O
+ NCSTEP=1;NCREAC=1;ICOMB=1
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FRATE*)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCH4)
+ STORE(P1RS)
WHEN EBU2,4
+ MESG(2 step 2 reactions finite-rate EBU model
+ MESG( 2CH4 + 3O2 > 2CO  + 4H2O
+ MESG( 2CO  + O2  > 2CO2
+ NCSTEP=2;NCREAC=2;ICOMB=2
  ** remove * from FRATE for weak linearisation
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FRATE)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,STWO)
+ SCRS(EBU,P,2.0);SCRS(EBU,S1,0.00)
   + STORE(P1RS,S1RS)
WHEN EBU3,4
+ MESG(2 step 3 reactions finite-rate EBU model
+ MESG( 2CH4 + O2  > 2CO  + 4H2
+ MESG( 2CO  + O2  > 2CO2
+ MESG( 2H2  + O2  > 2H2O
+ NCSTEP=2;NCREAC=3;ICOMB=3
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FRATE*)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCRS)
+ SCRS(EBU,P,2.0);SCRS(EBU,S1,0.0);SCRS(EBU,S2,1.0)
   + STORE(P1RS,S1RS)
WHEN BURN,4
+ MESG(Infinite rate model
+ MESG(2CH4 +  4O2 > 2CO2 + 4H2O
+ NCSTEP=1;NCREAC=1;ICOMB=4
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FASTC)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCH4)
WHEN DDEL,4
+ MESG(Infinite-rate model with Double-Delta PDF
+ MESG(2CH4 +  4O2 > 2CO2 + 4H2O
+ NCSTEP=1;NCREAC=1;ICOMB=5
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FASTC)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCH4)
+ SCRS(PDF,DDELTA)
ENDCASE
STORE(MMWT)
IF(ICOMB.EQ.3) THEN
+ YCH4F=0.2;YH2F=0.05
ENDIF
   ** Define fuel & oxidiser composition & temperature
SCRS(FUIN,YCH4F,YO2F,YH2F,YCOF,YH2OF,YCO2F,YN2F,TFUEL)
SCRS(OXIN,YCH4O,YO2O,YH2O,YCOO,YH2OO,YCO2O,YN2O,TOX)
  *** END OF EXTENDED SCRS MODEL SETTINGS
ENDIF
IF(KINET) THEN
+ THRAD=F
ELSE
+ MESG( Thermal radiation required ? (default=N)
+ READVDU(ANS2,CHAR,N)
IF(:ANS2:.EQ.Y) THEN
+ THRAD=T
ELSE
+ THRAD=F
ENDIF
IF(THRAD) THEN
+ REAL(ABSORB,SCAT,SIGMA,EMPW,EMISW,EMISG,EMPG)
+ ABSORB=0.5;SCAT=0.; EMISG=0.07
+ SIGMA=5.6697E-8; EMISW=0.4
+ EMPW=SIGMA*TWAL**4; EMPG=SIGMA*EMISG
+ RADIAT(RADI,ABSORB,SCAT,H1)
+ SOLUTN(SRAD,P,P,Y,P,P,P);SOLUTN(H1,P,P,Y,P,P,P)
ENDIF
ENDIF
  Group 11.Initialise Var/Porosity Fields
CONPOR(QUARL,0.00,VOLUME,#1,#1,#2,#2,#1,#2)
CONPOR(LIP,0.00,CELL,1,1,5,5,1,4)
IF(THRAD) THEN
+ REAL(TGUESS);TGUESS=1000.; FIINIT(SRAD)=0.07*SIGMA*TGUESS**4
+ FIINIT(H1)=-2.392E+06
ENDIF
 
IF(KINET) THEN
+ RESTRT(ALL)
+ DO KK=CHSOB+6,CHSOB+16
+   FIINIT(:KK:)=1.E-10
+ ENDDO
MESG( Do you want the CHEMKIN solver? (default=N)
READVDU(ANS,CHAR,N)
  
IF(:ANS:.EQ.Y) THEN
+ TWOPNT=T;CHEMKIN(REACT,TWOPNT)
  ** the following statement prevent termination of EARTH run
     by the TWOPNT solver
+ SPEDAT(SET,CHEM,BYPASS,L,T)
+ MESG( CHEMKIN  solver selected.
ELSE
+ TWOPNT=F;CHEMKIN(REACT,PHOENICS)
+ MESG( PHOENICS solver selected.
ENDIF
+ CHSOC=37.5
+ CHEMKIN(PROP,STAT,CHSOC)
+ RHO1=GRND9;TMP1=0.0
  TMP1B=13
ELSE
+ FIINIT(U1)=1.E-5; FIINIT(V1)=1.E-5; FIINIT(W1)=WIN
+ FIINIT(KE)=KEIN; FIINIT(EP)=EPIN
ENDIF
  Group 13. Boundary & Special Sources
 
   * deaktivert fast legeme rotasjon
   COVAL(SWIRL:JJ:, UC1 , ONLYMS, USWIRL*JJ/10)
IF(KINET) THEN
DO JJ =  1,4
+ INLET(NOCPCKO:JJ:,LOW,1,1,JJ,JJ,1,1,1,LSTEP)
+ VALUE(NOCPCKO:JJ:,P1,GRND9); VALUE(NOCPCKO:JJ:,W1,WSWIRL)
+ VALUE(NOCPCKO:JJ:,UCRT,0.0); VALUE(NOCPCKO:JJ:,VCRT,0.0)
+ VALUE(NOCPCKO:JJ:,TMP1,TOX)
+ VALUE(NOCPCKO:JJ:,O2,YO2O); VALUE(NOCPCKO:JJ:,H2O,YH2OO)
ENDDO
+ INLET(NOCPCKF,LOW,1,1,6,9,1,1,1,LSTEP)
+ VALUE(NOCPCKF,P1,GRND9); VALUE(NOCPCKF,W1,WIN)
+ VALUE(NOCPCKF,TMP1,TFUEL)
+ VALUE(NOCPCKF,H2O,YH2OF); VALUE(NOCPCKF,CO,YCOF)
+ VALUE(NOCPCKF,CH4,YCH4F); VALUE(NOCPCKF,CO2,YCO2F)
ELSE
DO JJ =  1,4
+ INLET(SCRSO:JJ:,LOW,1,1,JJ,JJ,1,1,1,LSTEP)
+ VALUE(SCRSO:JJ:,P1,GRND3); VALUE(SCRSO:JJ:,U1,USWIRL)
+ VALUE(SCRSO:JJ:,V1,GRND3); VALUE(SCRSO:JJ:,W1,GRND3)
+ VALUE(SCRSO:JJ:,UCRT,0.0); VALUE(SCRSO:JJ:,VCRT,0.0)
+ VALUE(SCRSO:JJ:,WCRT,WSWIRL)
+ VALUE(SCRSO:JJ:,KE,KESWIRL); VALUE(SCRSO:JJ:,EP,EPSWIRL)
+ VALUE(SCRSO:JJ:,H1,GRND3)
+ VALUE(SCRSO:JJ:,F,0.0); VALUE(SCRSO:JJ:,CH4,YCH4O)
+ VALUE(SCRSO:JJ:,CO,YCOO); VALUE(SCRSO:JJ:,H2,YH2O)
ENDDO
+ INLET(SCRSF,LOW,1,1,6,9,1,1,1,LSTEP)
+ VALUE(SCRSF,P1,GRND3); VALUE(SCRSF,KE,KEIN)
+ VALUE(SCRSF,EP,EPIN); VALUE(SCRSF,H1,GRND3)
+ VALUE(SCRSF,F,1.0); VALUE(SCRSF,CH4,YCH4F); VALUE(SCRSF,CO,YCOF)
+ VALUE(SCRSF,U1,0.0); VALUE(SCRSF,V1,GRND3); VALUE(SCRSF,W1,GRND3)
+ VALUE(SCRSF,UCRT,0.0); VALUE(SCRSF,VCRT,0.0);VALUE(SCRSF,WCRT,WIN)
+ VALUE(SCRSF,H2,YH2F)
ENDIF
  ## Coefficient at outlet controls pressure distribution
  ** Exit Boundary Conditions
PATCH(OUT,HIGH,#1,#1,#1,#2,#3,#3,#1,#1)
REAL(CMOUT);CMOUT=1.0
IF(KINET) THEN
  ** convert CMOUT from s/m**3 to g*s/(kg*cm**3)
+ COVAL(OUT,P1,CMOUT*1.E-3,0.0)
+ DO KK=CHSOB,CHSOB+16
+ VALUE(OUT,:KK:,0.0)
+ ENDDO
IF(TWOPNT) THEN
+ DO KK=CHSOB,CHSOB+16
+ RESREF(:KK:)=-1.E-12
+ ENDIT(:KK:)=1.E-6
+ ENDDO
ENDIF
ELSE
+ VALUE(OUT,F,0.0);COVAL(OUT,P1,CMOUT,0.0)
+ DO KK=CHSOB,CHSOB+12
+ VALUE(OUT,:KK:,0.0)
+ ENDDO
ENDIF
 
PATCH(OUTER,NWALL,1,NX,NY,NY,#3,#NREGZ,1,1)
COVAL(OUTER,KE,GRND2,GRND2);COVAL(OUTER,EP,GRND2,GRND2)
COVAL(OUTER,U1,GRND2,0);COVAL(OUTER,W1,GRND2,0)
 
PATCH(QUARLN,NWALL,#1,#1,#1,#1,#1,#2,1,1)
COVAL(QUARLN,KE,GRND2,GRND2);COVAL(QUARLN,EP,GRND2,GRND2)
COVAL(QUARLN,U1,GRND2,0);COVAL(QUARLN,W1,GRND2,0)
 
PATCH(QUARLL,LWALL,#1,#1,#2,#2,#3,#3,1,1)
COVAL(QUARLL,KE,GRND2,GRND2);COVAL(QUARLL,EP,GRND2,GRND2)
COVAL(QUARLL,U1,GRND2,0);COVAL(QUARLL,V1,GRND2,0)
 
PATCH(I_W_S,SWALL,1,1,5,5,1,4,1,LSTEP)
COVAL(I_W_S,KE,GRND2,GRND2);COVAL(I_W_S,EP,GRND2,GRND2)
COVAL(I_W_S,U1,GRND2,0);COVAL(I_W_S,V1,GRND2,0)
COVAL(I_W_S,W1,GRND2,0)
 
PATCH(S_W_N,NWALL,1,1,5,5,1,4,1,LSTEP)
COVAL(S_W_N,KE,GRND2,GRND2);COVAL(S_W_N,EP,GRND2,GRND2)
COVAL(S_W_N,U1,GRND2,0);COVAL(S_W_N,V1,GRND2,0)
COVAL(S_W_N,W1,GRND2,0)
 
PATCH(L_W_L,LWALL,1,1,5,5,4,4,1,LSTEP)
COVAL(L_W_L,KE,GRND2,GRND2);COVAL(L_W_L,EP,GRND2,GRND2)
COVAL(L_W_L,U1,GRND2,0);COVAL(L_W_L,V1,GRND2,0)
COVAL(L_W_L,W1,GRND2,0)
 
IF(THRAD) THEN
+ PATCH(NWALL3R,NORTH,1,NX,#NREGY,#NREGY,#3,#NREGZ,#1,#NREGT)
+ COVAL(NWALL3R,SRAD,EMISW/(2.0-EMISW),EMPW)
ENDIF
 
UCONV = T
PATCH(SWRC,CELL,1,NX,1,NY,1,NZ,1,1);COVAL(SWRC,U1,FIXFLU,GRND)
PATCH(SWRV,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(SWRV,V1,FIXFLU,GRND);COVAL(SWRV,W1,FIXFLU,GRND)
  Group 15. Terminate Sweeps
LSWEEP  =  2000
  Group 17. Relaxation
IF(KINET) THEN
IF(TWOPNT) THEN
+ CHEMKIN(RELAX,1.E-4)
ELSE
+ CHEMKIN(RELAX,1.E-2)
+ SELREF=T
ENDIF
  ** LG(94)=T so as to convert density from kg/m**3 to g/cm**3 (by
     multiplying by 1.E-3) for use in CHEMKIN on 1st sweep of
     kinetics run. However, on kinetics restart runs you must set
     LG(94)=F
+ RELAX(DEN1,LINRLX,1.E-10)
MESG( RE-STARTING from SCRS combustion solution? (default=Y)
READVDU(ANS,CHAR,Y)
IF(:ANS:.EQ.Y) THEN
+ LG(94)=T
+ MESG( RESTRT run from SCRS solution
ELSE
+ LG(94)=F; RESTRT(CH3,CH2,CH,CH2O,HCO,H,O,OH,HO2,H2O2,N2)
+ MESG( RESTRT run from KINETICS solution
ENDIF
ELSE
+ RELAX(P1,LINRLX,0.3); RELAX(U1,FALSDT,RELX); RELAX(V1,FALSDT,RELX)
+ RELAX(W1,FALSDT,RELX); RELAX(KE,FALSDT,RELX);RELAX(EP,FALSDT,RELX)
+ RELAX(F,LINRLX,0.5); RELAX(RHO1,LINRLX,0.5); RELAX(H1,LINRLX,0.5)
IF(THRAD) THEN
+ RELAX(SRAD,FALSDT,1.0)
ENDIF
+ KELIN=1
IF(ICOMB.LE.3) THEN
+ RELAX(CH4,FALSDT,RELX);RELAX(CO,FALSDT,RELX);RELAX(H2,FALSDT,RELX)
ENDIF
IF(ICOMB.EQ.3) THEN
+ RELAX(RHO1,LINRLX,0.05)
ENDIF
ENDIF
  Group 18. Limits
  Group 19. EARTH Calls To GROUND Station
  Group 20. Preliminary Printout
  OUTPUT(SPH1,P,P,Y,P,P,P)
  Group 22. Monitor Print-Out
IXMON=1;IYMON=4;IZMON=24;TSTSWP=-1;ITABL=3
  Group 23.Field Print-Out & Plot Control
  Group 24. Dumps For Restarts
  ** Activate thermal radiation after solution has
     developed; following commands required
IF(THRAD) THEN
IF(KINET) THEN
ELSE
+ MESG( Enter mode of radiation calculation
+ MESG( Default: Restart from non-radiative SCRS run
+ MESG( The options are:
+ MESG(  RAD1  - Re-starting from non-radiative SCRS run
+ MESG(  RAD2  - Re-starting from radiative SCRS run
+ MESG(  RAD3  - Starting from scratch
+ READVDU(CMOD,CHAR,RAD1)
CASE :CMOD: OF
WHEN RAD1,4
+ MESG( RE-STARTING from non-radiative SCRS combustion run
+ RESTRT(P1,U1,V1,W1,KE,EP,RHO1,H1,CH4,CO2,H2,N2)
+ RESTRT(H2O,CO2,CO,O2,F,ENUT,UCRT,VCRT,WCRT,TMP1)
+ FIINIT(SRAD)=0.07*SIGMA*TGUESS**4
WHEN RAD2,4
+ MESG( RE-STARTING from radiative SCRS combustion run
+ RESTRT(P1,U1,V1,W1,KE,EP,RHO1,H1,CH4,CO2,H2,N2)
+ RESTRT(H2O,CO2,CO,O2,F,ENUT,UCRT,VCRT,WCRT,TMP1,SRAD)
WHEN RAD3,4
+ MESG( Starting from scratch
+ FIINIT(SRAD)=0.07*SIGMA*TGUESS**4
ENDCASE
ENDIF
ENDIF
NZPRIN=1;NYPRIN=1;IZPRF=12;IZPRL=16;NPLT=100
   IDISPA=10000;CSG1=SW