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



  GR OU X    21 Y     1     3 Z     1     1
  GR OU X     1 Y     1     3 Z     1     1
  GR OU Y     4 X     1    20 Z     1     1
  GR OU Y     1 X     1    20 Z     1     1
  GR OU X     6 Y     4    10 Z     1     1
  GR OU X     1 Y     4    10 Z     1     1
  GR OU Y    11 X     1     5 Z     1     1
  GR OU Y     4 X     1     5 Z     1     1
  GR OU X     4 Y    11    14 Z     1     1
  GR OU X     1 Y    11    14 Z     1     1
  GR OU Y    15 X     1     3 Z     1     1
  GR OU Y    11 X     1     3 Z     1     1
  GR OU X     6 Y    15    22 Z     1     1
  GR OU X     1 Y    15    22 Z     1     1
  GR OU Y    23 X     1     5 Z     1     1
  GR OU Y    15 X     1     5 Z     1     1
  GR OU X    21 Y    15    22 Z     1     1
  GR OU X    14 Y    15    22 Z     1     1
  GR OU Y    23 X    14    20 Z     1     1
  GR OU Y    15 X    14    20 Z     1     1
  GR OU X    21 Y    11    14 Z     1     1
  GR OU X    16 Y    11    14 Z     1     1
  GR OU Y    15 X    16    20 Z     1     1
  GR OU Y    11 X    16    20 Z     1     1
  GR OU X    21 Y     9    10 Z     1     1
  GR OU X    14 Y     9    10 Z     1     1
  GR OU Y    11 X    14    20 Z     1     1
  GR OU Y     9 X    14    20 Z     1     1
  GR OU X    16 Y    11    14 Z     1     1 COL    14
  GR OU X     4 Y    11    14 Z     1     1 COL    14
  GR OU Y    15 X     4    15 Z     1     1 COL    14
  GR OU Y    11 X     4    15 Z     1     1 COL    14
  MSG                     COOLER GEOMETRY
  msg Press Enter to continue
  PAUSE;red
  VEC Z 1 Y 1 28 SH
  MSG                     Velocity vectors
  msg Press Enter to continue
  pause;vec off;red
  con tgas z 1 Y 1 28 fil;0.01
  msg                    Temperature contours
  msg Press Enter to continue
  pause;con off;red
  con yn2 z 1 Y 1 28 fil;0.01
  msg                  Nitrogen
  msg Press Enter to continue
  pause;con off;red
  con yo2 z 1 Y 1 28 fil;0.01
  msg                  Oxygen
  msg Press Enter to continue
  pause;con off;red
  con yh2 z 1 Y 1 28 fil;0.01
  msg                  Hydrogen
  msg Press Enter to continue
  pause;con off;red
  con yco z 1 Y 1 28 fil;0.01
  msg                   CO
  msg Press Enter to continue
  pause;con off;red
  con yco2 z 1 Y 1 28 fil;0.01
  msg                   CO2
  msg Press Enter to continue
  pause;con off;red
  con yh2o z 1 Y 1 28 fil;0.01
  msg                   H2O
  msg Press Enter to end
  ENDUSE

   ************************************************************
   *  GROUP 1.  Run identifiers and other preliminaries.
TEXT(Environmental level coke combustion simulation

  DISPLAY

  Environmental level coke combustion simulation

  The packed bed of heat generating materials is cooled by air
  flow generated by buyoancy.  The pressure drop across a packed
  bed is obtained from Ergun equation. The material tempeartures
  are also calculated.

  ENDDIS
              ==================
REAL(TENV,PORI,QDVOL,VVOL); TENV=-40.; PORI=0.8
VVOL=5.2*4.62*2.345
QDVOL=20.E6/VVOL*(1-PORI)
REAL(DIAM,LENGTH);DIAM=0.02
LENGTH=DIAM*PORI/(1.-PORI)
   ************************************************************
   *  GROUP 2.  Time-dependence and related parameters.
   ************************************************************
   *  GROUP 3.  x-direction grid specification.
NX      =      25
XULAST  = 14.4
XFRAC (  1) =  4.630E-02 ;XFRAC (  2) =  9.259E-02
XFRAC (  3) =  1.389E-01 ;XFRAC (  4) =  1.649E-01
XFRAC (  5) =  1.910E-01 ;XFRAC (  6) =  2.179E-01
XFRAC (  7) =  2.448E-01 ;XFRAC (  8) =  2.717E-01
XFRAC (  9) =  2.986E-01 ;XFRAC ( 10) =  3.255E-01
XFRAC ( 11) =  3.524E-01 ;XFRAC ( 12) =  3.793E-01
XFRAC ( 13) =  4.063E-01 ;XFRAC ( 14) =  4.323E-01
XFRAC ( 15) =  4.583E-01 ;XFRAC ( 16) =  4.972E-01
XFRAC ( 17) =  5.361E-01 ;XFRAC ( 18) =  5.750E-01
XFRAC ( 19) =  6.139E-01 ;XFRAC ( 20) =  6.528E-01
XFRAC ( 21) =  6.838E-01 ;XFRAC ( 22) =  7.406E-01
XFRAC ( 23) =  8.142E-01 ;XFRAC ( 24) =  9.012E-01
XFRAC ( 25) =  1.000E+00
   ************************************************************
   *  GROUP 4.  y-direction grid specification.
NY      =      30
YVLAST  = 27.0
YFRAC (  1) =  6.173E-02 ;YFRAC (  2) =  1.235E-01
YFRAC (  3) =  1.852E-01 ;YFRAC (  4) =  2.022E-01
YFRAC (  5) =  2.193E-01 ;YFRAC (  6) =  2.363E-01
YFRAC (  7) =  2.533E-01 ;YFRAC (  8) =  2.704E-01
YFRAC (  9) =  2.778E-01 ;YFRAC ( 10) =  2.852E-01
YFRAC ( 11) =  3.065E-01 ;YFRAC ( 12) =  3.278E-01
YFRAC ( 13) =  3.491E-01 ;YFRAC ( 14) =  3.704E-01
YFRAC ( 15) =  4.028E-01 ;YFRAC ( 16) =  4.352E-01
YFRAC ( 17) =  4.676E-01 ;YFRAC ( 18) =  5.000E-01
YFRAC ( 19) =  5.324E-01 ;YFRAC ( 20) =  5.648E-01
YFRAC ( 21) =  5.972E-01 ;YFRAC ( 22) =  6.296E-01
YFRAC ( 23) =  6.759E-01 ;YFRAC ( 24) =  7.222E-01
YFRAC ( 25) =  7.685E-01 ;YFRAC ( 26) =  8.148E-01
YFRAC ( 27) =  8.611E-01 ;YFRAC ( 28) =  9.074E-01
YFRAC ( 29) =  9.537E-01 ;YFRAC ( 30) =  1.000E+00
   *  GROUP 7.  Variables (including porosities) named,
   *            stored & solved.
SOLVE(P1,U1,V1,H1,METL); SOLUTN(P1,Y,Y,Y,N,N,N)
STORE(TMP1)
STORE(ENUT,EPOR,NPOR,VPOR)
   ************************************************************
   *  GROUP 8.  Terms (in differential equations) and devices.
TERMS(METL,N,N,N,N,Y,N)
   ************************************************************
   *  GROUP 9.  Properties of the medium (or media).
PRESS0=1.0000E+05
   *  CP1   = Specific heat (J/kg K)
REAL(CAP);CAP=996.
   *  TREF1 = Reference temperature (K)
REAL(TREF1);TREF1=273.
   *  HREF1 = Reference enthalpy (J/kg)
REAL(HREF1);HREF1=CAP*(TREF1+TENV)
   *  VEXP1 = Volume expansivity coefficient (1/K)
REAL(VEXP1);VEXP1=3.6720E-03
   *  Density (kg/m^3)
RHO1=1.2200E+00
   *  Set TMP1 = TREF + H1/CP1
TMP1=GRND2;TMP1A=0.0;TMP1B=1/CAP; CP1=CAP
   *  Prandtl number of the fluid/media
PRNDTL(H1)=7.3600E-01
   *  Laminar kinematic viscosity (m^2/s)
ENUL=1.4650E-05
TURMOD(KEMODL)
   *  Turbulent Prandtl number of the fluid/media
PRT(H1)=1
   ************************************************************
   *  GROUP 10. Interphase-transfer processes and properties.
   ************************************************************
   *  GROUP 11. Initialization of fields of variables,
   *            porosities, etc.
CONPOR(0,CELL,1,-20,1,-3,1,1)
CONPOR(0,CELL,1,-5,4,10,1,1)
CONPOR(0,CELL,1,-3,11,14,1,1)
CONPOR(0,CELL,1,-5,15,-22,1,1)
CONPOR(0,CELL,-14,-20,15,-22,1,1)
CONPOR(0,CELL,-16,-20,11,14,1,1)
CONPOR(0,CELL,-14,-20,-9,10,1,1)
CONPOR(PORI,VOLUME,4,15,11,14,1,1)
CONPOR(PORI,EAST,4,15,11,14,1,1)
CONPOR(PORI,NORTH,4,15,11,14,1,1)

FIINIT(U1)=0.0
FIINIT(V1)=0.1
FIINIT(H1)=HREF1
FIINIT(KE)=5.0E-02
FIINIT(EP)=2.4573E-01
   *
   ************************************************************
   *  GROUP 12. Convection and diffusion adjustments
   ************************************************************
   *  GROUP 13. Boundary conditions and special sources
OUTLET(OUTR,EAST,NX,NX,1,NY,1,NZ,1,1)
VALUE(OUTR,U1,SAME)
VALUE(OUTR,V1,SAME)
VALUE(OUTR,P1,0.)
VALUE(OUTR,H1,CAP*(TREF1+TENV))
VALUE(OUTR,KE,SAME)
VALUE(OUTR,EP,SAME)

OUTLET(OUTL,WEST,1,1,23,30,1,NZ,1,1)
VALUE(OUTL,U1,SAME)
VALUE(OUTL,V1,SAME)
VALUE(OUTL,P1,0.)
VALUE(OUTL,H1,CAP*(TREF1+TENV))
VALUE(OUTL,KE,SAME)
VALUE(OUTL,EP,SAME)
   * OUTLET boundary condition, name OUTUP
OUTLET(OUTUP,NORTH,1,NX,NY,NY,1,NZ,1,1)
VALUE(OUTUP,U1,SAME)
VALUE(OUTUP,V1,SAME)
VALUE(OUTUP,P1,0.)
VALUE(OUTUP,H1,CAP*(TREF1+TENV))
VALUE(OUTUP,KE,SAME)
VALUE(OUTUP,EP,SAME)
   * HEAT-SOURCE boundary condition, name HOTHEAT
PATCH(HEAT,VOLUME,4,15,11,14,1,1,1,1)
COVAL(HEAT,H1,FIXFLU,QDVOL)
COVAL(HEAT,METL,FIXFLU,QDVOL)
   ** Ergun's resitance law
      * Linear resistance term
PATCH(ERLIN,PHASEM,4,15,11,14,1,1,1,1)
COVAL(ERLIN,U1,0.1*ENUL/LENGTH**2,0.0)
COVAL(ERLIN,V1,0.1*ENUL/LENGTH**2,0.0)
      * Quadratic resistance term
PATCH(ERQUA,PHASEM,4,15,11,14,1,1,1,1)
    CO=RG(2)*ABS(U1)
COVAL(ERQUA,U1,GRND,0.0)
    CO=RG(2)*ABS(V1)
COVAL(ERQUA,V1,GRND,0.0)
   * Metal temperature
PATCH(METEMP,VOLUME,4,15,11,14,1,1,1,1)
    VAL=TMP1
COVAL(METEMP,METL,100.,GRND)
   * GRAVITY boundary condition, name BUOYANCY
PATCH(BUOYANCY,PHASEM,1,NX,1,NY,1,1,1,1)
COVAL(BUOYANCY,U1,FIXFLU,GRND3)
COVAL(BUOYANCY,V1,FIXFLU,GRND3)
   *  Set gravity resolutes
BUOYA=0.;BUOYB=-9.8100E+00;BUOYC=0.
   *  Set contants for Boussinesq approximation
BUOYD=VEXP1/CAP;BUOYE=HREF1
   *  GROUP 15. Termination criteria for sweeps and
LSWEEP=300
   *  GROUP 17. Under-relaxation and related devices.
REAL(MAXV,MINL,RELX)
MAXV=4.
MINL=5.7600E-01
RELX=1
RELAX(P1,LINRLX,0.8)
RELAX(U1,FALSDT,MINL/MAXV*RELX)
RELAX(V1,FALSDT,MINL/MAXV*RELX)
RELAX(H1,FALSDT,MINL/MAXV*1.E3)
RELAX(KE,FALSDT,MINL/MAXV*RELX/10)
RELAX(EP,FALSDT,MINL/MAXV*RELX/10)
   *  GROUP 19. Data communicated by SATELLITE to GROUND
NAMSAT=MOSG
RG(1)=150.*ENUL/LENGTH**2
RG(2)=0.0175/LENGTH
   *  GROUP 22. Location of spot-value & frequency of
   *            residual printout.
IXMON=14;  IYMON=13

tstswp=-1

   Coal-combustion model.

   Reactions: C (s) + 0.5 O2  > CO
              CO    + 0.5 O2  > CO2
              C(s)  + CO2     > 2CO
              C(s)  + H2O     > CO  + H2
              H2    + 0.5 O2  > H2O

      The numerical values which appear below are:-
      * the mass fractions of O2 and N2 in air, namely: 0.232 and 0.
      * the molecular weights of O2, N2, CO2, CO, H2O, H2, -
                  32.0, 28.0, 44.0, 28.0, 18.0 and 2.0
      * the Universal Gas Constant,                    GASCON: 8.314
      * the heat of reaction for    C  + O2     -> CO2, HCCO2: 3.279
      * the heat of reaction for    C  + 0.5*O2 -> CO , HCCO : 9.208
      * the heat of reaction for    H2 + 0.5*O2 -> H2O, HHH2O: 1.209
      * the specific heat at constant pressure,         CP   : 1.100
                                                        all in SI un

      CINCL, the mass fraction of C in the coal;
      HINCL, the mass fraction of H in the fuel;
      NINCL, the mass fraction of N in the fuel;
      AINCL, the mass fraction of ash in the fuel.


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

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, WIN)
  ** 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=3.

  ** FCL - carbon element mass fraction;
  ** The specific heat at constant pressure, CP =1.100E3

    H = CP*T + HCHX*YCHX + HCOCO2*YCO * HHH2*YH2

SOLVE(H11,FCL)

REAL(CP,TFUEL,HGIN,TGIN)
TGIN = TENV+273.
HGIN  = CAP*TGIN

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

SOLUTN(FCL,Y,Y,Y,P,P,P)
SOLUTN(H11  ,Y,Y,Y,P,P,P)

    GROUP 8. Terms (in differential equations) & devices
TERMS(FCL,N,Y,N,P,P,P)
TERMS(H11,N,Y,N,P,P,P)

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

    GROUP 11. Initialization of variable or porosity fields
FIINIT(FCL)=1.
    GROUP 13. Boundary conditions and special sources
  ** Inlet Boundaries
VALUE(OUTL,FCL, 0.0)
VALUE(OUTL,H11, HGIN)

VALUE(OUTR,FCL , 0.0)
VALUE(OUTR,H11, HGIN)

VALUE(OUTUP,FCL, 0.0)
VALUE(OUTUP,H11, HGIN)

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

PATCH(CARGAS,VOLUME,1,NX,1,NY,1,NZ,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(CARGAS,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(CARGAS,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 Tgas.
       - HSUB = HCOCO2*YCO * HH2*YH2
     VAL=:CP:*TGAS+:HCHX:+HSUB
COVAL(CARGAS,H11,ONLYMS,GRND)


    GROUP 15. Termination of sweeps
LSWEEP=350

    GROUP 16. Termination of iterations
LITHYD=10
VARMAX(FCL)=FS;VARMIN(FCL)=0.0
VARMIN(TGAS)=TGIN
    GROUP 17. Under-relaxation devices
RELAX(FCL,FALSDT,0.25)
RELAX(H11,FALSDT,0.25)


    GROUP 23. Field print-out and plot control
NPLT=1;NYPRIN=1;NZPRIN=1
NYPRIN=1;IYPRF=1;IYPRL=30
TSTSWP=-1

      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:))
    GO=:AIRO2:*(1-FCL)
    GC=:CINCL:*FCL
    GH=:HINCL:*FCL
    GOPA=GC*:MO2:/(2*:MC:)/(1-GO+GC*:MO2:/(2*:MC:)+TINY)
    GOFU=(GH*:MO2:/(2*:MH2:)+GC*:MO2:/:MC:)/$
            (1.-GO+GH*:MO2:/(2*:MH2:)+GC*:MO2:/:MC:+TINY)
    FRAC=(GO-GOPA)/(GOFU-GOPA+TINY)

           ** For all regions
              ---------------
    YN2=:NINCL:*FCL+:AIRN2:*(1.-FCL)

           ** Region 1
              --------
    YH2O=:HINCL:*FCL*:MH2O:/:MH2:
   IF(FCL.LE.FLIM)
    YCO2=:CINCL:*FCL*:MCO2:/:MC:
   IF(FCL.LE.FLIM)
    YO2 =:AIRO2:*(1-FCL)-:CINCL:*FCL*:MO2:/:MC:-$
                 :HINCL:*FCL*:MO2:/(2.*:MH2:)
   IF(FCL.LE.FLIM)
    YCO=0.0
   IF(FCL.LE.FLIM)
    YH2=0.0
   IF(FCL.LE.FLIM)
    HSUB=0.0
   IF(FCL.LE.FLIM)
  RMIX=8314.3*(YO2/32.+YH2O/18.+YCO2/44.+YN2/28.)
   IF(FCL.LE.FLIM)

           ** Region 2
              --------
    YH2O=:HINCL:*FCL*:MH2O:/:MH2:*FRAC*(1-GOFU)/(1-GO+TINY)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
    YCO2=:CINCL:*FCL*:MCO2:/:MC:*FRAC*(1-GOFU)/(1-GO+TINY)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
    YO2=0.0
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
    YCO=:CINCL:*FCL*:MCO:/:MC:*(1-FRAC)*$
                 (1-GOPA)/(1-GO+TINY)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
    YH2=:HINCL:*FCL*(1-FRAC)*(1-GOPA)/(1-GO+TINY)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
    HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)
  RMIX=8314.3*(YH2O/18.+YCO/28.+YCO2/44.+YH2/2.+YN2/28.)
   IF(FCL.GT.FLIM.AND.FRAC.GE.0.)

           ** Region 3
              --------
    YH2O=0.0
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
    YCO2=0.0
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
    YO2=0.0
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
    YCO=:AIRO2:*(1-FCL)*2*:MCO:/:MO2:
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
    YH2=:HINCL:*FCL
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
    HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)
   RMIX=8314.3*(YCO/28.+YH2/2.+YN2/28.)
   IF(FCL.GT.FLIM.AND.FRAC.LT.0.)

     ** Calculation of absolute gas temperature
        --------------------------------------
    TGAS=AMAX1(:TGIN:.,(H11-HSUB)/:CAP:.)
store(ysum)
    YSUM=YN2+YO2+YCO+YCO2+YH2O+YH2
store(DEN1)
    DEN1=1.e5/(RMIX*TGAS+tiny)


FIINIT(YO2)=0.232;FIINIT(YN2)=0.768

MAXV=4.
MINL=5.7600E-01
RELX=1
RELAX(P1,LINRLX,0.8)
RELAX(U1,FALSDT,MINL/MAXV*RELX)
RELAX(V1,FALSDT,MINL/MAXV*RELX)
RELAX(H1,FALSDT,MINL/MAXV*1.E3)
RELAX(KE,FALSDT,MINL/MAXV*RELX/10)
RELAX(EP,FALSDT,MINL/MAXV*RELX/10)
burnrate=1.
dmpstk=t
 LIBREF=262
STOP