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



     up z
     con temp x 1 z 1 16 fil;.001
     con burn x 1 z 18 m fil;.001
     gr ou x 1 z 1 16
     gr ou x 1 z 18 m
     msg SAFIR sub-model: Combustion of descending coke
     msg ---------------------------------------------
     msg Coke burning rate in upper domain:
     msg
     msg Temperature of gas in lower domain
     pause;cl
     con vpor x 1 z 1 16 fil;.001
     con vpor x 1 z 18 m fil;.001
     gr ou x 1 z 1 16
     gr ou x 1 z 18 m
     msg SAFIR sub-model: Combustion of descending coke
     msg ---------------------------------------------
     msg Volume fraction of coke in upper domain:
     msg
     msg Volume fraction  of gas in lower domain
     pause;cl
     set vec ref 0.01
     vec x 1 z 18 m sh
     set vec ref 700.
     vec x 1 z 1 16 sh
     gr ou x 1 z 1 16
     gr ou x 1 z 18 m
     msg   Solid velocities in upper domain
     msg   and
     msg   Gas velocities in lower domain
     enduse
 ************************************************************
  Group 1. Run Title
 TEXT(Combustion-driven coke flow)

  DISPLAY

    Combustion-driven coke flow:

    One-phase, 2D, two spaces, provision for raceway.

  ENDDIS

 ************************************************************
  Group 2. Transience
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    F
 NX      =       1
 XULAST  = 1.000E-01
 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NY      =       15
 YVLAST  = 8.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)

 ************************************************************
  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
     DEN1=RG(1)
  REGION() 1
  IF(ISWEEP.LE.50)
     DEN1=1.e5/(RMIX*TEMP+tiny)
  REGION() 1
  IF(ISWEEP.GT.50)
     DEN1=RG(2)
  REGION() 2
       ** 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
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
 INIADD  =    F
       ** The subdomains are isolated by one slab blockage
INTEGER(UPLOW)
UPLOW=17;IG(1)=UPLOW
STORE(HPOR,NPOR,VPOR,SPOR)
PATCH(ISOLAT,INIVAL,1,1,1,NY,UPLOW,UPLOW,1,1)
INIT( ISOLAT,VPOR, 0.000E+00, 0.0)
INIT( ISOLAT,HPOR, 0.000E+00, 0.0)
INIT( ISOLAT,NPOR, 0.000E+00, 0.0)
PATCH(ISOL,INIVAL,1,1,1,NY,UPLOW-1,UPLOW-1,1,1)
INIT( ISOL,HPOR, 0.000E+00, 0.0)

       ** Mark the upper and lower domains
           MARK=1 - lower domain
           MARK=2 - upper domain
STORE(MARK)
PATCH(LOWMARK,INIVAL,1,1,1,NY,1,UPLOW-1,1,1)
INIT( LOWMARK,MARK, 0.000E+00, 1.)
PATCH(UPPMARK,INIVAL,1,1,1,NY,UPLOW+1,NZ,1,1)
INIT( UPPMARK,MARK, 0.000E+00, 2.)
       **  Initalisations
STORE(HPOR,VPOR,NPOR)
    Lower domain
PATCH(INILOW,INIVAL,1,1,1,NY,1,UPLOW-1,1,1)
         *   Porosities
INIT( INILOW,VPOR, 0.000E+00, 0.25)
INIT( INILOW,HPOR, 0.000E+00, 0.25)
INIT( INILOW,NPOR, 0.000E+00, 0.25)
INIT( INILOW,SPOR, 0.000E+00, 0.25)
INIT( INILOW,RHO1, 0.000E+00, DEN1L)
    Upper domain
PATCH(INIUP,INIVAL,1,1,1,NY,UPLOW+1,NZ,1,1)
         *   Porosities
INIT( INIUP,VPOR, 0.000E+00, 0.75)
INIT( INIUP,HPOR, 0.000E+00, 0.75)
INIT( INIUP,NPOR, 0.000E+00, 0.75)
INIT( INIUP,SPOR, 0.000E+00, 0.75)
INIT( INIUP,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=20.
       **  Mass inflow rates of gas
PATCH(GASINL   ,NORTH ,1,NX,NY,NY,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,NY,UPLOW-1,UPLOW-1,1,1)
COVAL (TOPLOW     ,P1  , DEN1L*FIXP, 0.000E+00)
       **  Fix pressure coke inlet
PATCH(SOLINL   ,HIGH ,1,NX,1,NY,NZ,NZ,1,1)
COVAL(SOLINL   ,P1   , FIXP*5000.,.0)
       **  Frictional momentum loss for gas flow
PATCH (FRIC     ,VOLUME  ,1,NX,1,NY,1,UPLOW-1,1,1)
      CO=100.*(1.-VPOR)
COVAL (FRIC     ,V1  , GRND, 0.000E+00)
      CO=6.*(1.-VPOR)
COVAL (FRIC     ,W1  , GRND, 0.000E+00)
       **  gravity
  PATCH (GRAV     ,PHASEM  ,1,1,1,NY,UPLOW+1,NZ,1,1)
  COVAL (GRAV     ,W1  , FIXFLU, -9.8)
       ** 3D storage for Heat transfer
STORE(HTC)
FIINIT(HTC)=10.
       **  Heat exchange between lower  and upper sub-domains
  PATCH(SS001TEM,VOLUME,1,NX,1,NY,1,NZ,1,1)
  SORC10 CO =HTC
  SORC10 VAL=H1[,,+IG(1)]
  COVAL(SS001TEM,H1,GRND,GRND)
       **  Heat exchange between upper and lower sub-domains
  PATCH(SS002TEM,VOLUME,1,NX,1,NY,1,NZ,1,1)
  SORC11 CO =HTC
  SORC11 VAL=H1[,,-IG(1)]
  COVAL(SS002TEM,H1,GRND,GRND)
  ***********************************************************
  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=uplow-1
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)
  ** 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 burning
burnrate=5.0

SOLVE(H1,FCL,POR)

REAL(CP,TFUEL,HGIN,TGIN)
TGIN = 350.0
CP=  1100.
HGIN  = CP*TGIN

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,Y,P,P,P)
SOLUTN(H1  ,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)

    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
    GROUP 13. Boundary conditions and special sources
COVAL(GASINL,FCL , ONLYMS,0.0)
COVAL(GASINL,H1  , ONLYMS,HGIN)
COVAL(GASINL,POR , ONLYMS,0.0)

  Carbon mass transfer related sources:
  ------------------------------------

PATCH(gasFcoke,VOLUME,1,NX,1,NY,1,UPLOW-1,1,1)

  (1) Transfer of mass leading to increase of gas flow rate:
       - VPOR is volume fraction of lump coal
     VAL=:BURNRATE:*(1.-VPOR)*(:FS:-FCL)
COVAL(gasFcoke,P1,FIXFLU,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)

PATCH(coke2gas,VOLUME,1,NX,1,NY,UPLOW+1,NZ,1,1)
     VAL=(-:BURNRATE:*VPOR*(:FS:-FCL[,,-IG(1)]))/1.e-20
COVAL(coke2gas,P1,1.e-20,GRND)


PATCH(RACEWAY,VOLUME,1,NX,1,NY,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
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.15)
RELAX(W1,FALSDT,.01)
RELAX(V1,FALSDT,.01)
RELAX(FCL,FALSDT,.01)
RELAX(POR,FALSDT,.01)
RELAX(H1,FALSDT,.01)
RELAX(RHO1,LINRLX,0.15)

      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,NY,1,IG(1)-1)
       GO=:AIRO2:*(1-FCL)
   REGION(1,NX,1,NY,1,IG(1)-1)
       GC=:CINCL:*FCL
   REGION(1,NX,1,NY,1,IG(1)-1)
       GH=:HINCL:*FCL
   REGION(1,NX,1,NY,1,IG(1)-1)
       GOPA=GC*:MO2:/(2*:MC:)/(1-GO+GC*:MO2:/(2*:MC:)+TINY)
   REGION(1,NX,1,NY,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,NY,1,IG(1)-1)
       FRAC=(GO-GOPA)/(GOFU-GOPA+TINY)
   REGION(1,NX,1,NY,1,IG(1)-1)

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

           ** Region 1
              --------
       YH2O=:HINCL:*FCL*:MH2O:/:MH2:
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YCO2=:CINCL:*FCL*:MCO2:/:MC:
   REGION(1,NX,1,NY,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,NY,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YCO=0.0
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       YH2=0.0
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       HSUB=0.0
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.LE.FLIM)
       RMIX=:GASCON:*(YO2/:MO2:+YH2O/:MH2O:+YCO2/:MCO2:+$
                    YN2/:MN2:)
   REGION(1,NX,1,NY,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,NY,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,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
       YO2=0.0
   REGION(1,NX,1,NY,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,NY,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,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
       HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   REGION(1,NX,1,NY,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,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)

           ** Region 3
              --------
        YH2O=0.0
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YCO2=0.0
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YO2=0.0
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YCO=:AIRO2:*(1-FCL)*2*:MCO:/:MO2:
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        YH2=:HINCL:*FCL
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
        RMIX=:GASCON:*(YCO/:MCO:+YH2/:MH2:+YN2/:MN2:)
   REGION(1,NX,1,NY,1,IG(1)-1)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)

     ** Calculation of absolute gas temperature
        --------------------------------------
       TEMP=(H1-HSUB)/:CP:.
   REGION(1,NX,1,NY,1,IG(1)-1)
       TEMP=AMIN1(4500.,AMAX1(350.,TEMP,350.))
   REGION(1,NX,1,NY,1,IG(1)-1)
store(ysum,burn,VVPO)
       YSUM=YN2+YO2+YCO+YCO2+YH2O+YH2
   REGION(1,NX,1,NY,1,IG(1)-1)
       BURN=:BURNRATE:*(1.-VPOR)*(:FS:-FCL[,,-IG(1)])
   REGION(1,NX,1,NY,IG(1)+1,NZ)
       VVPO=SPOR+(1.-POR)*(1.-SPOR)
   REGION(1,NX,1,NY,1,IG(1)-1)

   VVPO=1.
       REGION(1,NX,1,NY,1,IG(1)-1)
       IF(VVPO.GT.0.8.AND.LG(1))

   VVPO=SPOR
       REGION(1,NX,1,NY,1,IG(1)-1)
       IF(VVPO.LT.0.8.AND.LG(1))

   VPOR=VVPO
       REGION(1,NX,1,NY,1,IG(1)-1)
   VPOR=1.-VVPO[,,-IG(1)]
       REGION(1,NX,1,NY,IG(1)+1,NZ)

    NPOR=VVPO
       REGION(1,NX,1,NY,1,IG(1)-1)
    NPOR=1.-VVPO[,,-IG(1)]
       REGION(1,NX,1,NY,IG(1)+1,NZ)

    HPOR=VVPO
       REGION(1,NX,1,NY,1,IG(1)-2)
    HPOR=1.-VVPO[,,-IG(1)]
       REGION(1,NX,1,NY,IG(1)+1,NZ)

inifld=t
fiinit(temp)=tgin
lsweep=700
lg(1)=f
STOP