talk=t;run(1,1)
  PHOTON USE
     p



     up z
     con yore x 1 y 1 5 z 18 m fil;.001
     pause
     con ycok x 1 y 1 5 z 18 m fil;.001
     pause
     con temp x 1 y 1 5 z 1 16 fil;.001
     con temp x 1 y 1 5 z 18 m fil;.001
     *con yore x 1 y 1 5 z 18 m fil;.001
     *pause
     *con ycok x 1 y 1 5 z 18 m fil;.001
     *pause
     con temp x 1 y 7 m z 1 16 fil;.001
     con temp x 1 y 7 m z 18 m fil;.001
     gr ou x 1 y 1 5 z 18 m
     gr ou x 1 y 1 5 z 1 16 m
     gr ou x 1 y 7 m z 1 16
     gr ou x 1 y 7 m z 18 m
     msg
     msg
     msg   Liquid temperatures                         Solid tempera
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg   Fines temperatures                          Gas Temperatu
     *enduse
     pause;cl
     set con scale range on
     con flim x 1 y 7 m z 1 16 fil
     0.847e-04 0.958e-04
     .001
     set con scale range off
     con flim x 1 y 7 m z 18 m fil;.001
     con vpor  x 1 y 1 5 z 1 16 fil;.001
     con vpor x 1 y 1 5 z 18 m fil;.001
     gr ou x 1 y 7 m z 1 16
     gr ou x 1 y 7 m z 18 m
     gr ou x 1 y 1 5 z 1 16
     gr ou x 1 y 1 5 z 18 m
     msg
     msg
     msg   Solid fraction                            Solid + liquid
     msg   in solid/liquid                           volume
     msg   mixture                                   fraction
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg   Fines volume                                Gas volume
     msg   fraction                                    fraction
     pause;cl
     gr ou x 1 y 7 m z 1 16
     gr ou x 1 y 7 m z 18 m
     gr ou x 1 y 1 5 z 18 m
     gr ou x 1 y 1 5 z 1 16
     set vec ref 0.005
     vec x 1 y 1 5 z 18 m sh
     set vec ref 0.005
     vec x 1 y 7 m z 18 m sh
     set vec ref 800.
     vec x 1 y 1 5 z 1 15 sh
     set ve comp - v2 w2
     vec x 1 y 7 m z 1 15 sh
     set vec comp - v1 w1
     msg
     msg
     msg   Liquid velocities                           Solid velocit
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg
     msg   Fines velocities                            Gas velocitie
     pause;cl
     gr ou x 1 y 1 5 z 1 16
     con yo2 x 1 y 1 5 z 1 16 fil;.001
     msg   Oxygen
     pause;cl
     gr ou x 1 y 1 5 z 1 16
     con yn2 x 1 y 1 5 z 1 16 fil;.001
     msg   Nitrogen
     pause;cl
     gr ou x 1 y 1 5 z 1 16
     con yh2 x 1 y 1 5 z 1 16 fil;.001
     msg   Hydrogen
     pause;cl
     gr ou x 1 y 1 5 z 1 16
     con yco2 x 1 y 1 5 z 1 16 fil;.001
     msg Carbon dioxide
     pause;cl
     gr ou x 1 y 1 5 z 1 16
     con yh2o x 1 y 1 5 z 1 16 fil;.001
     msg Water vapour
     pause;cl
     gr ou x 1 y 1 5 z 1 16
     con yco x 1 y 1 5 z 1 16 fil;.001
     msg Carbon monoxide
     enduse
 ************************************************************
  Group 1. Run Title
 TEXT(Combustion-melting-driven ore/coke mixture flow)

  DISPLAY

  Combustion-melting-driven ore/coke  flow:

  One-phase, 2D, four spaces, scalar treatment of coal
  fines combustion, lump coke combustion, solid
  fusion(melting), provision for combustion-driven raceway,
  gas/solid and solid/liquid heat transfer)

  ENDDIS

 ************************************************************
  Group 2. Transience
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    F
 NX      =       1
 XULAST  = 1.000E-0
 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NY      =       11
 YVLAST  = 22.000E+00
 GRDPWR(Y,NY,YVLAST,1.)
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 PARAB   =    F
 NZ      =      33
 ZWLAST  = 5.700E+01
 GRDPWR(Z,NZ,ZWLAST,1.)
    * Cylindrical-polar grid
 CARTES=F
 ************************************************************
  Group 6. Body-Fitted coordinates
    * X-cyclic boundaries switched
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
    * Solved variables list
 SOLVE(P1  ,V1  ,W1 )
    * Additional solver options
 SOLUTN(P1  ,Y,Y,y,N,N,Y)
 SOLUTN(v1  ,Y,Y,p,N,N,Y)
 SOLUTN(w1  ,Y,Y,p,N,N,Y)

 ************************************************************
  Group 8. Terms & Devices
 DENPCO  =    T
 ************************************************************
  Group 9. Properties
       ** Introduce the densities
REAL(DEN1L,DEN1U)
    Gas
DEN1L=1.    ;  RG(1)=DEN1L
    Solid
DEN1U=1000. ;  RG(2)=DEN1U

STORE(RHO1)

RHO1=GRND
  PRPT01 DEN1=RG(1)
  REGION() 1 /ISWEEP.LE.50

   DEN1=1.e5/(RMIX*TEMP+tiny)
  REGION() 1 /ISWEEP.GE.1

     DEN1=RG(2)
  REGION() 2
     DEN1=RG(2)
  REGION() 3
       ** Introduce the viscosities
           VISL=1.e-05  for lower domain
           VISL=1.0     for upper domain
STORE(VISL)
ENUL=GRND
     VISL=1.e-05
  REGION() 1
     VISL=1.
  REGION() 2
     VISL=1.
  REGION() 3
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
INIADD=F
       ** The subdomains are isolated by one slab blockage
INTEGER(UPLOW,WESEA)
UPLOW=17;IG(1)=UPLOW
WESEA=6 ;IG(2)=WESEA
STORE(HPOR,NPOR,VPOR,SPOR,VVPO)
PATCH(ISOLAT1,INIVAL,1,1,1,NY,UPLOW,UPLOW,1,1)
INIT( ISOLAT1,VPOR, 0.000E+00, 0.0)
INIT( ISOLAT1,HPOR, 0.000E+00, 0.0)
INIT( ISOLAT1,NPOR, 0.000E+00, 0.0)
PATCH(ISOL1,INIVAL,1,1,1,NY,UPLOW-1,UPLOW-1,1,1)
INIT( ISOL1,HPOR, 0.000E+00, 0.0)

PATCH(ISOLAT2,INIVAL,1,1,WESEA,WESEA,1,NZ,1,1)
INIT( ISOLAT2,VPOR, 0.000E+00, 0.0)
INIT( ISOLAT2,HPOR, 0.000E+00, 0.0)
INIT( ISOLAT2,NPOR, 0.000E+00, 0.0)
PATCH(ISOL2,INIVAL,1,1,WESEA-1,WESEA-1,1,NZ,1,1)
INIT( ISOL2,NPOR, 0.000E+00, 0.0)

       ** Mark the upper and lower domains
           MARK=1 - gas space
           MARK=2 - solid space
           MARK=3 - liquid space
           MARK=4 - fines space
STORE(MARK)
PATCH(GASMARK,INIVAL,1,1,1,WESEA-1,1,UPLOW-1,1,1)
INIT( GASMARK,MARK, 0.000E+00, 1.)
PATCH(SOLMARK,INIVAL,1,1,1,WESEA-1,UPLOW+1,NZ,1,1)
INIT( SOLMARK,MARK, 0.000E+00, 2.)
PATCH(LIQMARK,INIVAL,1,1,WESEA+1,NY,UPLOW+1,NZ,1,1)
INIT( LIQMARK,MARK, 0.000E+00, 3.)
PATCH(FINMARK,INIVAL,1,1,WESEA+1,NY,1,UPLOW-1,1,1)
INIT( FINMARK,MARK, 0.000E+00, 4.)
       **  Initalisations
STORE(HPOR,VPOR,NPOR)
        * Gas space
PATCH(INIGAS,INIVAL,1,1,1,WESEA-1,1,UPLOW-1,1,1)
INIT( INIGAS,VPOR, 0.000E+00, 0.25)
INIT( INIGAS,HPOR, 0.000E+00, 0.25)
INIT( INIGAS,NPOR, 0.000E+00, 0.25)
INIT( INIGAS,SPOR, 0.000E+00, 0.25)
INIT( INIGAS,VVPO, 0.000E+00, 0.25)
INIT( INIGAS,RHO1, 0.000E+00, DEN1L)
        * Solid space
PATCH(INISOL,INIVAL,1,1,1,WESEA-1,UPLOW+1,NZ,1,1)
INIT( INISOL,VPOR, 0.000E+00, 0.75)
INIT( INISOL,HPOR, 0.000E+00, 0.75)
INIT( INISOL,NPOR, 0.000E+00, 0.75)
INIT( INISOL,SPOR, 0.000E+00, 0.75)
INIT( INISOL,VVPO, 0.000E+00, 0.75)
INIT( INISOL,RHO1, 0.000E+00, DEN1U)
        * Liquid space
PATCH(INILIQ,INIVAL,1,1,WESEA+1,NY,UPLOW+1,NZ,1,1)
INIT( INILIQ,VPOR, 0.000E+00, 0.75)
INIT( INILIQ,HPOR, 0.000E+00, 0.75)
INIT( INILIQ,NPOR, 0.000E+00, 0.75)
INIT( INILIQ,SPOR, 0.000E+00, 0.75)
INIT( INILIQ,VVPO, 0.000E+00, 0.75)
INIT( INILIQ,RHO1, 0.000E+00, DEN1U)
 ************************************************************
  Group 12. Convection and diffusion adjustments
   No PATCHes used for this Group
 ************************************************************
  Group 13. Boundary & Special Sources
REAL(VIN)
       ** Blast gas velocity
VIN=150.
       **  Mass inflow rates of gas
PATCH(GASINL   ,NORTH ,1,NX,WESEA-1,WESEA-1,2,2,1,1)
COVAL(GASINL   ,P1   , FIXFLU, DEN1L*VIN)
COVAL(GASINL   ,V1   , ONLYMS, -VIN)
       **  Fix pressure gas outlet
PATCH (TOPLOW     ,HIGH  ,1,NX,1,WESEA-1,UPLOW-1,UPLOW-1,1,1)
COVAL (TOPLOW     ,P1  , DEN1L*FIXP, 0.000E+00)
       **  Fix pressure coke inlet
PATCH(SOLINL   ,HIGH ,1,NX,1,WESEA-1,NZ,NZ,1,1)
COVAL(SOLINL   ,P1   , FIXP*5000.,.0)
       **  Frictional momentum loss for gas flow
PATCH (FRIC     ,VOLUME  ,1,NX,1,WESEA-1,1,UPLOW-1,1,1)
           CO=1.*3.e2*VPOR
COVAL (FRIC     ,V1  , GRND, 0.000E+00)
           CO=1.*3.e2*VPOR
COVAL (FRIC     ,W1  , GRND, 0.000E+00)
       **  Frictional momentum loss for liquid flow
PATCH (FRICLIQ     ,VOLUME  ,1,NX,WESEA+1,NY,UPLOW+1,NZ,1,1)
      CO=1.e5*FLIM
COVAL (FRICLIQ     ,V1  , GRND, 0.000E+00)
      CO=6.*FLIM
COVAL (FRICLIQ     ,W1  , GRND, 0.000E+00)
       ** 3D storage for Heat transfer
STORE(HTC)
  ***********************************************************
  Group 14. Downstream Pressure For PARAB
 ************************************************************
  Group 15. Terminate Sweeps
 LSWEEP  =     750
 ************************************************************
  Group 16. Terminate Iterations
 ************************************************************
  Group 17. Relaxation

 RELAX(P1  ,LINRLX, 0.75)
 RELAX(V1  ,FALSDT, 10.)
 RELAX(W1  ,FALSDT, 10.)
 RELAX(H1  ,FALSDT, 1000.)
 ************************************************************
  Group 22. Monitor Print-Out
 IXMON   =       1 ;IYMON  =       5 ;IZMON  =       6
 NPRMNT  =       1
 ************************************************************
  Group 23.Field Print-Out & Plot Control
 YZPR    =    T
   No PATCHes used for this Group
 ************************************************************
  Group 24. Dumps For Restarts

    ==========================================
NAMSAT=MOSG
nyprin=1;nzprin=1
iymon=ny/2;izmon=nz-2
lsweep=300
TSTSWP=-1
    ===========================================

REAL(HINCL,CINCL,NINCL,AINCL,GASCON)
REAL(AIRO2,AIRN2)
REAL(MN2,MC,MO2,MH2,MCO,MCO2,MH2O)

GASCON=8.3143e3

HINCL=0.05
CINCL=0.95
NINCL=0.0
AINCL=1.-CINCL-HINCL-NINCL

AIRO2=0.232
AIRN2=0.768

MN2=28.; MC=12.; MO2=32.; MH2=2.; MCO=28.; MCO2=44.; MH2O=18.

REAL(HCCO2,HCCO,HHH2O,HCHX,HCOCO2)
HCCO2 =3.279E7
HCCO  = 9.208E6
HHH2O = 1.209E6

REAL(FS,BURNRATE,LIQRATE,FINRATE,RFIN)
  ** FS is the mass of fuel per unit mass of air/fuel mixture
     to convert all carbon and oxygen to carbon monoxide.
FS=0.232/(0.232 + CINCL*16.0/12.0)
  ** The heat of combustion per unit mass of co is hcco2 minus hcco
     the mass of c per unit mass of co, ie 12/28
HCOCO2=(12.0/28.0)*(HCCO2-HCCO)
  ** The heat of coal combustion per unit mass of carbon
HCHX=(CINCL*HCCO2+HINCL*HHH2O)*(HINCL*MH2+CINCL*MC+NINCL*MN2)/MC
  ** The rate of coke burning
BURNRATE=0.4
  ** The rate of melting
LIQRATE=7.5e-04
  ** The rate of fines burning
FINRATE=1.e5

SOLVE(H1,FCL,POR,RF,YORE,YCOK)

REAL(CP,TFUEL,HGIN,TGIN,HSIN,TSIN)
REAL(TLIQ,TSOL,MIND,COKINSOL,OREINSOL)
TGIN = 600.0;TSIN=350.
CP=  1100.
HGIN  = CP*TGIN
HSIN  = CP*TSIN
RFIN=1.e-04
COKINSOL=0.4
OREINSOL=1.-COKINSOL

STORE(RMIX,HSUB,TEMP,YN2,YH2,YO2,YCO,YCO2,YH2O)
STORE(FLIM,FRAC,GO,GC,GH,GOFU,GOPA,RHO1)

SOLUTN(FCL ,Y,Y,y,P,P,P)
SOLUTN(POR ,Y,y,n,P,P,P)
SOLUTN(RF  ,Y,Y,y,P,P,P)
SOLUTN(H1  ,Y,Y,y,P,P,P)
SOLUTN(YCOK,Y,Y,y,P,P,P)
SOLUTN(YORE,Y,Y,y,P,P,P)

    GROUP 8. Terms (in differential equations) & devices
TERMS(FCL ,N,Y,N,P,P,P)
TERMS(POR ,N,Y,N,P,P,P)
TERMS(H1  ,N,Y,N,P,P,P)
TERMS(RF  ,N,Y,N,P,P,P)
TERMS(YCOK,N,Y,N,P,P,P)
TERMS(YORE,N,Y,N,P,P,P)

    GROUP 9. Properties of the medium (or media)
REAL(RHOIN1,WAIR)
PRESS0=1.e5
WAIR=32.
RHOIN1=PRESS0*WAIR/(8314.*TGIN)

    GROUP 11. Initialization of variable or porosity fields
FIINIT(FCL)=FS
FIINIT(POR)=1.0
FIINIT(POR)=RFIN
FIINIT(YORE)=OREINSOL
FIINIT(YCOK)=COKINSOL
INIT( INIGAS,H1, 0.000E+00,  HGIN)
INIT( INISOL,H1, 0.000E+00,  HSIN)
    GROUP 13. Boundary conditions and special sources
COVAL(GASINL,FCL , ONLYMS,0.0)
COVAL(GASINL,H1  , ONLYMS,HGIN)
COVAL(GASINL,POR , ONLYMS,0.0)
COVAL(GASINL,RF  , ONLYMS,RFIN)

COVAL(SOLINL,H1  , ONLYMS,HSIN)
COVAL(SOLINL,FCL , ONLYMS,SAME)
COVAL(SOLINL,POR , ONLYMS,SAME)
COVAL(SOLINL,RF  , ONLYMS,SAME)
COVAL(SOLINL,YORE, ONLYMS,OREINSOL)
COVAL(SOLINL,YCOK, ONLYMS,COKINSOL)

    * Liquid outlet from liquid space
PATCH(LIQOUT,LOW,1,NX,WESEA+1,NY,UPLOW+1,UPLOW+1,1,1)
COVAL(LIQOUT,P1 , FIXP*5000,0.0)


       **  Heat exchange between gas  and solid spaces
PATCH(SS001TEM,VOLUME,1,NX,1,NY,1,NZ,1,1)
   VAL=HTC*(TEMP[,,+IG(1)]-TEMP)*VPOR
COVAL(SS001TEM,H1,FIXFLU,GRND)
       **  Heat exchange between solid and gas spaces
PATCH(SS002TEM,PHASEM,1,NX,1,NY,1,NZ,1,1)
   VAL=HTC*(TEMP[,,-IG(1)]-TEMP)*VPOR
COVAL(SS002TEM,H1,FIXFLU,GRND)

REAL(HTLG)
HTLG=0.0
       **  Heat exchange between liquid  and solid spaces
PATCH(SS001TEM,VOLUME,1,NX,1,NY,1,NZ,1,1)
   VAL=:HTLG:*(TEMP[,+IG(2),]-TEMP)*VPOR
COVAL(SS001TEM,H1,FIXFLU,GRND)
       **  Heat exchange between solid and liquid spaces
PATCH(SS003TEM,VOLUME,1,NX,1,NY,1,NZ,1,1)
   VAL=:HTLG:*(TEMP[,-IG(2),]-TEMP)*VPOR
COVAL(SS003TEM,H1,FIXFLU,GRND)

  Carbon mass transfer related sources:
  ------------------------------------
      1. Gas from coke
PATCH(gasFcoke,VOLUME,1,NX,1,WESEA-1,1,UPLOW-1,1,1)
  (1) Transfer of mass leading to increase of gas flow rate:
       - VPOR is volume fraction of lump coal
     VAL=YCOK[,,+IG(1)]*BURN[,,+IG(1)]*VPOR[,,+IG(1)]/1.e-20
COVAL(gasFcoke,P1,1.e-20,GRND)
  (2) Transfer of carbon leading to increase of mixture
      fraction at the same rate:
       - CO=1. signifies that mass tarnsfer brings in
         material which is 100% carbon
COVAL(gasFcoke,FCL,ONLYMS,1.)
  (3) Transfer of enthalpy and heat leading to increase of
      gas enthalpy at the same rate:
       - Interphase gas temperature is assumed as TEMP.
       - HSUB = HCOCO2*YCO * HH2*YH2
     VAL=:CP:*TEMP+:HCHX:+HSUB
COVAL(gasFcoke,H1  ,ONLYMS,GRND)
COVAL(gasFcoke,POR ,ONLYMS,SAME)
COVAL(gasFcoke,RF  ,ONLYMS,SAME)
COVAL(gasFcoke,YORE,ONLYMS,SAME)
COVAL(gasFcoke,YCOK,ONLYMS,SAME)

      2. Coke to Gas
PATCH(coke2gas,VOLUME,1,NX,1,WESEA-1,UPLOW+1,NZ,1,1)
     VAL=-YCOK*BURN*VPOR
COVAL(coke2gas,P1  ,FIXFLU,GRND)
COVAL(coke2gas,RF  ,ONLYMS,SAME)
  CO = -YCOK*BURN*VPOR
  VAL=YORE*BURN[,+6,]*(1.-FLIM[,+6,])*YCOK/(BURN*YCOK+tiny)
COVAL(coke2gas,YORE,GRND,GRND)
  CO=-YORE*BURN[,+IG(2),]*VPOR*(1.-FLIM[,+IG(2),])
  VAL=YORE*YCOK*BURN/(BURN[,+6,]*(1.-FLIM[,+6,])*YORE+tiny)
COVAL(coke2gas,YCOK,GRND,GRND)

      3. Solid to Liquid
PATCH(sol2liq,VOLUME,1,NX,1,WESEA-1,UPLOW+1,NZ,1,1)
     VAL=-YORE*BURN[,+IG(2),]*VPOR*(1.-FLIM[,+IG(2),])
COVAL(sol2liq,P1,FIXFLU,GRND)

      4. Liquid From Coke
PATCH(liqFsol,VOLUME,1,NX,WESEA+1,NY,UPLOW+1,NZ,1,1)
  VAL=YORE[,-IG(2),]*BURN*VPOR[,-IG(2),]*(1.-FLIM)
COVAL(liqFsol,P1   ,FIXFLU,GRND)
COVAL(liqFsol,H1   ,ONLYMS,CP*TSOL)
COVAL(liqFsol,FCL  ,ONLYMS,SAME)
COVAL(liqFsol,POR  ,ONLYMS,SAME)
COVAL(liqFsol,YORE ,ONLYMS,SAME)
COVAL(liqFsol,YCOK,ONLYMS,SAME)

      5. Gas From Fines
PATCH(gasFfine,VOLUME,1,NX,1,WESEA-1,1,UPLOW-1,1,1)
     VAL=BURN*RF*VPOR
COVAL(gasFfine,P1 ,FIXFLU ,GRND)
COVAL(gasFfine,FCL,ONLYMS ,1.)
COVAL(gasFfine,RF ,ONLYMS ,0.)
     VAL=:CP:*TEMP+:HCHX:+HSUB
COVAL(gasFfine,H1  ,ONLYMS,GRND)
COVAL(gasFfine,POR ,ONLYMS,SAME)
COVAL(gasFfine,YORE,ONLYMS,SAME)
COVAL(gasFfine,YCOK,ONLYMS,SAME)

      6. Coke-combustion-driven raceway
PATCH(RACEWAY,VOLUME,1,NX,1,WESEA-1,1,UPLOW-1,1,1)
      CO=1.*BURN[,,+IG(1)]
      VAL=1./(1.*BURN[,,+IG(1)]+tiny)
COVAL(RACEWAY,POR,GRND,GRND)

    GROUP 16. Termination of iterations
LITHYD=10
VARMAX(FCL)=FS;VARMIN(FCL)=0.0
VARMIN(TEMP)=TGIN;VARMAX(TEMP)=3000.
VARMIN(RHO1)=0.001;VARMAX(RHO1)=DEN1U
VARMAX(POR)=1.0;VARMIN(POR)=0.0
VARMAX(RF) =RFIN;VARMIN(RF)=0.0
VARMAX(YORE)=1.;VARMIN(YORE)=0.0
VARMAX(YCOK)=1.;VARMIN(YCOK)=0.0
    GROUP 17. Under-relaxation devices
RELAX(P1  ,LINRLX,0.15)
RELAX(W1  ,FALSDT,0.01)
RELAX(V1  ,FALSDT,0.01)
RELAX(FCL ,FALSDT,0.01)
RELAX(POR ,FALSDT,0.01)
RELAX(H1  ,FALSDT,10.0)
RELAX(RHO1,LINRLX,0.15)
RELAX(RF  ,FALSDT,0.01)
RELAX(YORE,FALSDT,0.01)
RELAX(YCOK,FALSDT,0.01)

      Coal oxidation is presumed to proceed in two stages, viz:
        (1) to create CO2 and H2O, and then
        (2) to create CO and H2, as more fuel is added.

      The gas composition diagram, taking account of the elemental
      mass fractions of O, C and H, has got three regions:
      (1) Region 1  containing O2, CO2 & H2O
      (2) Region 2  containing CO2, H2O, H2 & CO
      (3) Region 3  containing H2 & CO.

      The values of oxygen fraction GO at which the formulae exbibit
      discontinuities of slope are called:-
        GOPA, where the oxygen has consumed part of the fuel, so as
                create CO and H2; and
        GOFU, where the products of combustion are CO2 and H2O.

           ** Cell-wise composition parameters
              --------------------------------
       FLIM=:AIRO2:/(:AIRO2:+:CINCL:*:MO2:/:MC:+$
                    :HINCL:*:MO2:/(2*:MH2:))
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       GO=:AIRO2:*(1-FCL)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       GC=:CINCL:*FCL
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       GH=:HINCL:*FCL
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       GOPA=GC*:MO2:/(2*:MC:)/(1-GO+GC*:MO2:/(2*:MC:)+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       GOFU=(GH*:MO2:/(2*:MH2:)+GC*:MO2:/:MC:)/$
               (1.-GO+GH*:MO2:/(2*:MH2:)+GC*:MO2:/:MC:+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       FRAC=(GO-GOPA)/(GOFU-GOPA+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)

           ** For all regions
              ---------------
       YN2=:NINCL:*FCL+:AIRN2:*(1.-FCL)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)

           ** Region 1
              --------
       YH2O=:HINCL:*FCL*:MH2O:/:MH2:
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YCO2=:CINCL:*FCL*:MCO2:/:MC:
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YO2 =:AIRO2:*(1-FCL)-:CINCL:*FCL*:MO2:/:MC:-$
                    :HINCL:*FCL*:MO2:/(2.*:MH2:)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YCO=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YH2=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)

           ** Region 2
              --------
      YH2O=:HINCL:*FCL*:MH2O:/2.*FRAC*(1-GOFU)/(1-GO+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
      YCO2=:CINCL:*FCL*:MCO2:/:MC:*FRAC*(1-GOFU)/(1-GO+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
       YO2=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
       YCO=:CINCL:*FCL*:MCO:/:MC:*(1-FRAC)*$
                    (1-GOPA)/(1-GO+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
       YH2=:HINCL:*FCL*(1-FRAC)*(1-GOPA)/(1-GO+TINY)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)

           ** Region 3
              --------
        YH2O=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YCO2=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YO2=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YCO=:AIRO2:*(1-FCL)*2*:MCO:/:MO2:
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YH2=:HINCL:*FCL
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)

     **********************************

    YH2 =AMAX1(0.,YH2 )
    YH2O=AMAX1(0.,YH2O)
    YCO =AMAX1(0.,YCO )
    YCO2=AMAX1(0.,YCO2)
    YN2 =AMAX1(0.,YN2 )
    YO2 =AMAX1(0.,YO2 )

       HSUB=0.0
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       RMIX=:GASCON:*(YO2/:MO2:+YH2O/:MH2O:+YCO2/:MCO2:+$
                    YN2/:MN2:)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.LE.FLIM)

       HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
       RMIX=:GASCON:*(YH2O/:MH2O:+YCO/:MCO:+YCO2/:MCO2:+$
                     YH2/:MH2:+YN2/:MN2:)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)

        HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        RMIX=:GASCON:*(YCO/:MCO:+YH2/:MH2:+YN2/:MN2:)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)

     ** Calculation of absolute gas and solid temperatures
        --------------------------------------------------
         * Solid
       TEMP=H1/:CP:.
   REGION(1,NX,1,IG(2)-1,IG(1)+1,NZ)
         * Gas
       TEMP=(H1-HSUB)/:CP:.
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       TEMP=AMIN1(4500.,AMAX1(350.,TEMP,350.))
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
         * Liquid
       TEMP=H1/:CP:.
   REGION(1,NX,IG(2)+1,NY,IG(1)+1,NZ)

     ** Auxiliary calculations
        ----------------------
STORE(YSUM,BURN)
         * Coke burning rate
       BURN=:BURNRATE:*(:FS:-FCL[,,-IG(1)])
   REGION(1,NX,1,IG(2)-1,IG(1)+1,NZ)
         * Fines burning rate
       BURN=:FINRATE:*(:FS:-FCL)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
         * Liquid rate
       BURN=:LIQRATE:*AMAX1(0.0,(TEMP[,-IG(2),]-:TSOL:))
   REGION(1,NX,IG(2)+1,NY,IG(1)+1,NZ)
         * Reference volume porosity in gas space
       VVPO=SPOR+(1.-POR)*(1.-SPOR)
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
   IF(LG(2))
         * Porosity limiters
       VVPO=1.
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
         IF(VVPO.GT.0.8.AND.LG(1))
       VVPO=SPOR
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
         IF(VVPO.LT.0.8.AND.LG(1))
         * Porosities in solid space
       VPOR=1.-VVPO[,,-IG(1)]
   REGION(1,NX,1,IG(2)-1,IG(1)+1,NZ)
       NPOR=1.-VVPO[,,-IG(1)]
   REGION(1,NX,1,IG(2)-1,IG(1)+1,NZ)
       HPOR=1.-VVPO[,,-IG(1)]
   REGION(1,NX,1,IG(2)-1,IG(1)+1,NZ)
         * Porosities in gas space
       VPOR=VVPO
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       NPOR=VVPO
   REGION(1,NX,1,IG(2)-1,1,IG(1)-1)
       HPOR=VVPO
   REGION(1,NX,1,IG(2)-1,1,IG(1)-2)
         * Volume fraction of melted solid
   Parameter           T liquidus        T solidus
+    MIND=1.       ;    TLIQ=1450.    ;   TSOL=1050.
     FLIM=AMIN1(1.,AMAX1(((:TLIQ:-TEMP)$
             /(:TLIQ:-:TSOL:))**:MIND:,0.0))
    REGION(1,NX,IG(2)+1,NY,IG(1)+1,NZ) /ISWEEP.GT.10
         ====================================
              Post-processing
STORE(V2,W2)
      * Fines velocities
     V2=V1[,-IG(2),]
    REGION(1,NX,IG(2)+1,NY,1,IG(1)-1) /ISWEEP.GT.200
     W2=W1[,-IG(2),]
    REGION(1,NX,IG(2)+1,NY,1,IG(1)-1) /ISWEEP.GT.200
      * Fines volume fraction
     FLIM=RF[,-IG(2),]
    REGION(1,NX,IG(2)+1,NY,1,IG(1)-1) /ISWEEP.GT.200
       * Fines temperature
     TEMP=TEMP[,-IG(2),]
    REGION(1,NX,IG(2)+1,NY,1,IG(1)-1) /isweep.gt.200
        * Sum of gas components mass fractions
     YSUM=YN2+YO2+YCO+YCO2+YH2O+YH2
    REGION(1,NX,1,IG(2)-1,1,IG(1)-1) /isweep.gt.200
         * Total solid rate
     BURN=BURN[,-IG(2),+IG(2)]+BURN[,,+IG(1)]
    REGION(1,NX,IG(2)+1,NY,1,IG(1)-1) /isweep.gt.200

fiinit(temp)=tgin
  restrt(all)
lsweep=1700

   ** Introduce the sharpness in raceway treatment
lg(1)=f
   ** Introduce the solid-coke-combustion-driven raceway
lg(2)=t

cartes=t

FIINIT(HTC)=25.e-1
HTLG=5.e2

RELAX(P1  ,LINRLX,0.50)
RELAX(RHO1,LINRLX,0.15)
RELAX(W1  ,FALSDT,0.01)
RELAX(V1  ,FALSDT,0.01)
RELAX(FCL ,FALSDT,0.01)
RELAX(POR ,FALSDT,0.01)
RELAX(H1  ,FALSDT,40.5)
RELAX(YORE,FALSDT,1000.01)
RELAX(YCOK,FALSDT,1000.01)

  ** The rate of coke burning
BURNRATE=4.0
  ** The rate of melting
LIQRATE=0.1*2.0*7.5e-04
  ** The rate of fines burning
FINRATE=1.e5

inifld=t
STOP