************************************************************
  Group 1. Run Title
TEXT(Low Load, Coal and Wood)
BOOLEAN(SLOPES,WWOOD)
     **** A special subroutine is needed to be linked into
          EARTH executable if SLOPES=T
SLOPES=F
WWOOD=T
    *** Burning rates
REAL(WOODBURN,CHARBURN,COALBURN)
WOODBURN=1.0E1;CHARBURN=1.0E1;COALBURN=1.0E4

REAL(DG2RD)
DG2RD=3.14159/180.0
    
    *** Total air flow
    *** Biasing of Total Primary Air Flow, kg/sec
    ---  Input Mill to Mill PA bias
REAL(abpa,bbpa,cbpa,dbpa)
REAL(apa,bpa,cpa,dpa)
    --- Input Corner to Corner Primary Air Biases
REAL(a1bpa,a2bpa,a3bpa,a4bpa)
REAL(b1bpa,b2bpa,b3bpa,b4bpa)
REAL(c1bpa,c2bpa,c3bpa,c4bpa)
REAL(d1bpa,d2bpa,d3bpa,d4bpa)
REAL(B1p1,B2p1,B3p1,B4p1)
REAL(B1p2,B2p2,B3p2,B4p2)
REAL(B1p3,B2p3,B3p3,B4p3)
REAL(B1p4,B2p4,B3p4,B4p4)
REAL(B1X1,B2X1,B3X1,B4X1)
REAL(B1X2,B2X2,B3X2,B4X2)
REAL(B1X3,B2X3,B3X3,B4X3)
REAL(B1X4,B2X4,B3X4,B4X4)
REAL(B1Z1,B2Z1,B3Z1,B4Z1)
REAL(B1Z2,B2Z2,B3Z2,B4Z2)
REAL(B1Z3,B2Z3,B3Z3,B4Z3)
REAL(B1Z4,B2Z4,B3Z4,B4Z4)
    ***  Input windbox flow biases, mass rate basis
REAL(wbx1bs,wbx2bs,wbx3bs,wbx4bs)
    ***  Calculate Individual Windbox Flows
REAL(wbx1sa,wbx2sa,wbx3sa,wbx4sa)
    ***  Fractional openings of fuel air 
REAL(f1abs,f2abs,f3abs,f4abs) 
REAL(f1bbs,f2bbs,f3bbs,f4bbs) 
REAL(f1cbs,f2cbs,f3cbs,f4cbs) 
REAL(f1dbs,f2dbs,f3dbs,f4dbs) 
    ***  Fractional openings of auxiliary air
REAL(a1bs1, a2bs1, a3bs1, a4bs1) 
REAL(a1bs2, a2bs2, a3bs2, a4bs2) 
REAL(a1bs3, a2bs3, a3bs3, a4bs3) 
REAL(a1bs4, a2bs4, a3bs4, a4bs4)
REAL(A1BS5, A2BS5, A3BS5, A4BS5)
REAL(a1bs6, a2bs6, a3bs6, a4bs6) 
REAL(a1bs7, a2bs7, a3bs7, a4bs7) 
REAL(a1bs8, a2bs8, a3bs8, a4bs8) 

REAL(a1Obs1, a2Obs1, a3Obs1, a4Obs1)
REAL(A1OBS2, A2OBS2, A3OBS2, A4OBS2)
REAL(a1Obs3, a2Obs3, a3Obs3, a4Obs3) 
    ***  Fractional openings of OFA 
REAL(A1OFB1,A2OFB1,A3OFB1,A4OFB1)
REAL(A1OFB2,A2OFB2,A3OFB2,A4OFB2)
    ***   Total air calculations
REAL(t1bs,t2bs,t3bs,t4bs)
    ***  Calculate Individual FA compartment flows, kg/sec
REAL(f1a1,f1a2,f1a3,f1a4,f1a5,f1a6,f1a7,f1a8)
REAL(f2a1,f2a2,f2a3,f2a4,f2a5,f2a6,f2a7,f2a8)
REAL(f3a1,f3a2,f3a3,f3a4,f3a5,f3a6,f3a7,f3a8)
REAL(f4a1,f4a2,f4a3,f4a4,f4a5,f4a6,f4a7,f4a8)
    ***  Calculate Auxiliary Air Flows
REAL(a1a1,A1A2,A1A3,A1A4,A1A5,A1A6,A1A7,A1A8)
REAL(a2a1,A2A2,A2A3,A2A4,A2A5,A2A6,A2A7,A2A8)
REAL(a3a1,A3A2,A3A3,A3A4,A3A5,A3A6,A3A7,A3A8)
REAL(a4a1,A4A2,A4A3,A4A4,A4A5,A4A6,A4A7,A4A8)

REAL(a1O1,A1O2,A1O3)
REAL(a2O1,A2O2,A2O3)
REAL(a3O1,A3O2,A3O3)
REAL(a4O1,A4O2,A4O3)

REAL(A1OFA1,A2OFA1,A3OFA1,A4OFA1)
REAL(A1OFA2,A2OFA2,A3OFA2,A4OFA2)
    *** Slope Variables to GROUND
REAL(SLTYPF,SLTYPB,SLTYPN)
REAL(PNT1X,PNT1Y,PNT2X,PNT2Y,PNT3X,PNT3Y)
    *** Temporary variable
REAL(FRAC)

    *** Overall Dimensions, X: Back to Front; Z: Right to Left
REAL(XWIDTH,YHIGHT,ZDEPTH)
XWIDTH=7.494
YHIGHT=20.18
ZDEPTH=9.875
    *** Bottom slopes
REAL(SLPANG,XSLPLN,YSLPLN)
XSLPLN=3.334
YSLPLN=4.375
    SLPANG=52.68 deg
SLPANG=0.9195
    *** Wind box dimentions
REAL(INLETX,INLETZ)
REAL(YINLET)
REAL(YBAXLN,YOSALN,YOILLN,YCOLLN,YFLALN)
REAL(YH,X0TEMP,Z0TEMP)
INLETX=0.826
INLETZ=0.856
YINLET=7.497
YBAXLN=0.3797
YOSALN=0.2427
YOILLN=0.2614
YCOLLN=0.2081
YFLALN=0.1
    *** OFA dimensions
REAL(YOFAL1,YOFAL2)
YOFAL1=0.5094
YOFAL2=0.5047
    *** Nose
REAL(XNOSLN,YNOSLN,NOSANG)
NOSANG=20.0*DG2RD
XNOSLN=3.4
YNOSLN=XNOSLN/COS(NOSANG)*SIN(NOSANG)

   *** Coal Properties **************************************

REAL(TWALL,NPAT,KEAX,KEFA,KEPA,KEWD,TWBX,TPIPE,RHOWBX)
NPAT=1
KEAX=0.12;KEFA=0.12;KEPA=0.12;KEWD=0.12
PRESS0=1.0000E+05
TWBX=544.
TPIPE=350.
RHOWBX=PRESS0/(287.41*TWBX)
TWALL=615.
   
REAL(FC,FH2,FN2,TNA,ALFC,ALFH2,ALFN2)
fh2 = 4.49;fc  =71.74;fn2 = 1.18
fc=fc/100.;fh2=fh2/100.;fn2=fn2/100.
  -- fraction of essential elements
tna=fc+fh2+fn2
  -- normalized mass fraction of carbon
alfc=fc/tna
  -- normalized mass fraction of hydrogen
alfh2=fh2/tna
  -- normalized mass fraction of nitrogen
alfn2=fn2/tna
REAL(HHVRT,HR1,HR2,HR3)
hhvrt=0.573
hr1=32.792E6*hhvrt;hr2=9.208E6;hr3=120.9E6*hhvrt

REAL(GFS,GFF,GHA,GHB,GHF,GHS,GHO,GHA1,GHA2,GALPHA,HINWBX)
gfs=0.232/(0.232+alfc*16.0/12.0)
gff=2.*0.232/18./(alfh2+2.*(alfn2-0.768+44.*alfc/12.)/18.)
gha =gff*(alfc*hr1+alfh2*hr3)
ghb =alfc*gfs*hr2
ghf =0.0
gho =ghb/(1.-gfs)
gha1=ghb*(1.-gff)/(1.-gfs)
gha2=gha-gha1
ghs =ghf+gha2

REAL(CP,HINCOAL,rhopipe,hinpipe)
cp=1.8e3
hinwbx = cp*(twbx-temp0)
hincoal= ghs+cp*(tpipe-temp0)
rhopipe=press0/(287.41*tpipe)
hinpipe= cp*(tpipe-temp0)

REAL(TC,TW,TA,TPA,TSA,TWPA)
   
   Total coal when full load with wood
TC=6.829
TC
   Total wood
TW=2.333
TW
   Total wood primary air
TWPA=8.182

REAL(TWOOD,MWOOD,TWA,RHOWA)
TWOOD=273+25;TWA=TW+TWPA
MWOOD=TW/TWA
RHOWA=PRESS0/(287.41*TWOOD)

TW=TW/4
TWPA=TWPA/4

REAL(ACF,BCF,CCF,DCF)
ACF=TC/4.0;BCF=TC/4.0;CCF=TC/4.0;DCF=TC/4.0

    ***  Total Air for full load with wood 
TA=93.
    
    *** Total Primary Air for Coal with full load wood
TPA=18.

    *** Wood properties
real(hhh2o,hcco2)
hhh2o=120.9e6;hcco2=32.792e6

real(hwood,hchar,hvol,hinwd,vlinwd,hinvl,cinvl,mlwtvl,moinwd)
vlinwd=0.5;hinvl=0.02;cinvl=1.0-hinvl;mlwtvl=100;moinwd=0.20                        
  
hvol =hhh2o*hinvl+hcco2*cinvl; hchar=hcco2
hinwd=(hhh2o*vlinwd*hinvl+hcco2*(1-vlinwd+vlinwd*cinvl))
hinwd=cp*twood+hinwd
   *  Input Corner to Corner Pipe Split Biases
    *** Coal Input Corner to Corner Pipe Split Biases
REAL(C1B1,C1B2,C1B3,C1B4)
REAL(C2B1,C2B2,C2B3,C2B4)
REAL(C3B1,C3B2,C3B3,C3B4)
REAL(C4B1,C4B2,C4B3,C4B4)
REAL(C1F1,C1F2,C1F3,C1F4)
REAL(C2F1,C2F2,C2F3,C2F4)
REAL(C3F1,C3F2,C3F3,C3F4)
REAL(C4F1,C4F2,C4F3,C4F4)

C1B1=1.0;C1B2=1.0;C1B3=1.0;C1B4=1.0
C2B1=1.0;C2B2=1.0;C2B3=1.0;C2B4=1.0
C3B1=1.0;C3B2=1.0;C3B3=1.0;C3B4=1.0
C4B1=1.0;C4B2=1.0;C4B3=1.0;C4B4=1.0
FRAC=acf/(C1B1+C1B2+C1B3+C1B4)
C1F1=C1B1*FRAC;C1F2=C1B2*FRAC;C1F3=C1B3*FRAC;C1F4=C1B4*FRAC
FRAC=Bcf/(C2B1+C2B2+C2B3+C2B4)
C2F1=C2B1*FRAC;C2F2=C2B2*FRAC;C2F3=C2B3*FRAC;C2F4=C2B4*FRAC
FRAC=Ccf/(C3B1+C3B2+C3B3+C3B4)
C3F1=C3B1*FRAC;C3F2=C3B2*FRAC;C3F3=C3B3*FRAC;C3F4=C3B4*FRAC
FRAC=Dcf/(C4B1+C4B2+C4B3+C4B4)
C4F1=C4B1*FRAC;C4F2=C4B2*FRAC;C4F3=C4B3*FRAC;C4F4=C4B4*FRAC

   -- Biasing of Total Primary Air Flow, kg/sec

   *  Input Mill to Mill PA bias

abpa=1.0;bbpa=1.0;cbpa=1.0;dbpa=1.0
frac=tpa/(abpa+bbpa+cbpa+dbpa)
apa=abpa*frac;bpa=bbpa*frac
cpa=cbpa*frac;dpa=dbpa*frac

   -- Input Corner to Corner Primary Air Biases

a1bpa=1. ;a2bpa=1. ;a3bpa=1. ;a4bpa=1.
b1bpa=1. ;b2bpa=1. ;b3bpa=1. ;b4bpa=1.
c1bpa=1. ;c2bpa=1. ;c3bpa=1. ;c4bpa=1.
d1bpa=1. ;d2bpa=1. ;d3bpa=1. ;d4bpa=1.
frac=apa/(a1bpa+a2bpa+a3bpa+a4bpa)
B1p1=a1bpa*frac;B2p1=a2bpa*frac
B3p1=a3bpa*frac;B4p1=a4bpa*frac
frac=bpa/(b1bpa+B2bpa+b3bpa+b4bpa)
B1p2=b1bpa*frac;B2p2=b2bpa*frac
B3p2=b3bpa*frac;B4p2=b4bpa*frac
frac=cpa/(c1bpa+B2bpa+c3bpa+c4bpa)
B1p3=c1bpa*frac;B2p3=c2bpa*frac
B3p3=c3bpa*frac;B4p3=c4bpa*frac
frac=dpa/(d1bpa+B2bpa+d3bpa+d4bpa)
B1p4=d1bpa*frac;B2p4=d2bpa*frac
B3p4=d3bpa*frac;B4p4=d4bpa*frac

   *  Input windbox flow biases, mass rate basis
wbx1bs=1.0;wbx2bs=1.0;wbx3bs=1.0;wbx4bs=1.0

   *  Calculate Individual Windbox Flows

TSA=TA-TPA-TWPA
FRAC=tsa/(wbx1bs+wbx2bs+wbx3bs+wbx4bs)
wbx1sa=wbx1bs*FRAC;wbx2sa=wbx2bs*FRAC
wbx3sa=wbx3bs*FRAC;wbx4sa=wbx4bs*FRAC

   *  Fractional openings of fuel and auxiliary air 
    -- Fuel Air when full load
f1abs=0.08; f2abs=0.08; f3abs=0.08; f4abs=0.08
f1bbs=0.08; f2bbs=0.08; f3bbs=0.08; f4bbs=0.08
f1cbs=0.08; f2cbs=0.08; f3cbs=0.08; f4cbs=0.08
f1dbs=0.08; f2dbs=0.08; f3dbs=0.08; f4dbs=0.08
    -- Auxiliary Air when full load
a1bs1=0.06; a2bs1=0.06; a3bs1=0.06; a4bs1=0.06
a1bs2=0.03; a2bs2=0.03; a3bs2=0.03; a4bs2=0.03
a1bs3=0.03; a2bs3=0.03; a3bs3=0.03; a4bs3=0.03
a1bs4=0.03; a2bs4=0.03; a3bs4=0.03; a4bs4=0.03
A1BS5=0.03; A2BS5=0.03; A3BS5=0.03; A4BS5=0.03
a1bs6=0.03; a2bs6=0.03; a3bs6=0.03; a4bs6=0.03
a1bs7=0.03; a2bs7=0.03; a3bs7=0.03; a4bs7=0.03
a1bs8=0.05; a2bs8=0.05; a3bs8=0.05; a4bs8=0.05
    -- Oil Gun
a1Obs1=0.02; a2Obs1=0.02; a3Obs1=0.02; a4Obs1=0.02
A1OBS2=0.0 ; A2OBS2=0.0 ; A3OBS2=0.0 ; A4OBS2=0.0
a1Obs3=0.02; a2Obs3=0.02; a3Obs3=0.02; a4Obs3=0.02 
    -- OFA
      0.11 when full load
a1OFB1=0.16; a2OFB1=0.16; a3OFB1=0.16; a4OFB1=0.16
      0.77 when full load
A1OFB2=0.16; A2OFB2=0.16; A3OFB2=0.16; A4OFB2=0.16

REAL(TOTAUX,TOTFA,TOTOFA)
TOTAUX=0.0;TOTFA=0.0;TOTOFA=0.0

    *** Rear-right wind box: #1
    *** Rear-left wind box: #2
    *** Front-left wind box: #3
    *** Front-right wind box: #4

DO ii=1, 4

   *   Total air calculations
t:ii:bs=YBAXLN*a:ii:bs1+2.0*YFLALN*F:ii:abs+YOSALN*a:ii:bs2
t:ii:bs=t:ii:bs+YOILLN*a:ii:Obs1+YOSALN*a:ii:bs3+2.0*YFLALN*f:ii:bbs
t:ii:bs=t:ii:bs+YOSALN*a:ii:bs4+YOILLN*a:ii:Obs2+YOSALN*a:ii:bs5
t:ii:bs=t:ii:bs+2.0*YFLALN*f:ii:cbs+YOSALN*a:ii:bs6+YOILLN*a:ii:Obs3
t:ii:bs=t:ii:bs+YOSALN*a:ii:bs7+2.0*YFLALN*f:ii:dbs
t:ii:bs=t:ii:bs+YOSALN*a:ii:bs8+YOFAL1*A:ii:OFB1+YOFAL2*A:ii:OFB2

   *  Calculate Individual FA compartment flows, kg/sec
FRAC=YFLALN*f:ii:abs/t:ii:bs
f:ii:a1=wbx:ii:sa*FRAC;f:ii:a2=wbx:ii:sa*FRAC
FRAC=YFLALN*f:ii:bbs/t:ii:bs
f:ii:a3=wbx:ii:sa*FRAC;f:ii:a4=wbx:ii:sa*FRAC
FRAC=YFLALN*f:ii:cbs/t:ii:bs
f:ii:a5=wbx:ii:sa*FRAC;f:ii:a6=wbx:ii:sa*FRAC
FRAC=YFLALN*f:ii:dbs/t:ii:bs
f:ii:a7=wbx:ii:sa*FRAC;f:ii:a8=wbx:ii:sa*FRAC

TOTFA=TOTFA+f:ii:a1+f:ii:a2+f:ii:a3+f:ii:a4
TOTFA=TOTFA+f:ii:a5+f:ii:a6+f:ii:a7+f:ii:a8

   *  Calculate Auxiliary Air Flows
a:ii:a1=wbx:ii:sa*YBAXLN*a:ii:bs1/t:ii:bs
a:ii:a2=wbx:ii:sa*YOSALN*a:ii:bs2/t:ii:bs
a:ii:a3=wbx:ii:sa*YOSALN*a:ii:bs3/t:ii:bs
a:ii:a4=wbx:ii:sa*YOSALN*a:ii:bs4/t:ii:bs
a:ii:a5=wbx:ii:sa*YOSALN*a:ii:bs5/t:ii:bs
a:ii:a6=wbx:ii:sa*YOSALN*a:ii:bs6/t:ii:bs
a:ii:a7=wbx:ii:sa*YOSALN*a:ii:bs7/t:ii:bs
a:ii:A8=wbx:ii:sa*YOSALN*a:ii:BS8/t:ii:bs

a:ii:O1=wbx:ii:sa*YOILLN*a:ii:Obs1/t:ii:bs
a:ii:O2=wbx:ii:sa*YOILLN*a:ii:Obs2/t:ii:bs
a:ii:O3=wbx:ii:sa*YOILLN*a:ii:OBS3/t:ii:bs

a:ii:OFA1=wbx:ii:sa*YOFAL1*a:ii:OFB1/t:ii:bs
a:ii:OFA2=wbx:ii:sa*YOFAL2*a:ii:OFB2/t:ii:bs

TOTAUX=TOTAUX+a:ii:a1+a:ii:a2+a:ii:a3+a:ii:a4+a:ii:a5+a:ii:a6
TOTAUX=TOTAUX+a:ii:a7+a:ii:a8+a:ii:O1+a:ii:O2+a:ii:O3

TOTOFA=TOTOFA+a:ii:OFA1+a:ii:OFA2

ENDDO
  ************************************************************
  Print out Air Flows
  ************************************************************
TA
TPA
TOTFA
TOTAUX
TOTOFA
 ************************************************************
  Groups 3, 4, 5  Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
    * Overall domain extent, RSET(D,name,XULAST,YVLAST,ZWLAST)
RSET(D,CHAM,XWIDTH,YHIGHT,ZDEPTH)
    * Set objects: name  x0       y0        z0
    *                    dx       dy        dz
    *** Bottom slopes
RSET(B,BCKSLP , 0.000E+00, 0.000E+00, 0.000E+00                 , $
XSLPLN, YSLPLN, ZWLAST)                                    
RSET(B,FNTSLP , XULAST-XSLPLN, 0.000E+00, 0.000E+00             , $
XSLPLN, YSLPLN, ZWLAST)
                          
    *** Wind boxes
    --- Rear-right wind box: #1
    --- Rear-left wind box: #2
    --- Front-left wind box: #3
    --- Front-right wind box: #4
DO ii=1, 4
CASE ii OF
WHEN 1
+   X0TEMP=0.0
+   Z0TEMP=0.0
WHEN 2
+   X0TEMP=0.0
+   Z0TEMP=ZWLAST-INLETZ
WHEN 3
+   X0TEMPTEMP=XULAST-INLETX
+   Z0TEMP=ZWLAST-INLETZ
WHEN 4
+   X0TEMP=XULAST-INLETX
+   Z0TEMP=0.0
ENDCASE
YH=YINLET
RSET(B,I:ii:LAX1,X0TEMP,YH,Z0TEMP,INLETX,YBAXLN,INLETZ)
YH=YH+YBAXLN
RSET(B,I:ii:LFA1,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:COL1,X0TEMP,YH,Z0TEMP,INLETX,YCOLLN,INLETZ)
YH=YH+YCOLLN
RSET(B,I:ii:UFA1,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:OSA1,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
RSET(B,I:ii:OIL1,X0TEMP,YH,Z0TEMP,INLETX,YOILLN,INLETZ)
YH=YH+YOILLN
RSET(B,I:ii:OSA2,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
RSET(B,I:ii:LFA2,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:COL2,X0TEMP,YH,Z0TEMP,INLETX,YCOLLN,INLETZ)
YH=YH+YCOLLN
RSET(B,I:ii:UFA2,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:OSA3,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
RSET(B,I:ii:OIL1,X0TEMP,YH,Z0TEMP,INLETX,YOILLN,INLETZ)
YH=YH+YOILLN
RSET(B,I:ii:OSA4,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
RSET(B,I:ii:LFA3,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:COL3,X0TEMP,YH,Z0TEMP,INLETX,YCOLLN,INLETZ)
YH=YH+YCOLLN
RSET(B,I:ii:UFA3,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:OSA5,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
RSET(B,I:ii:OIL1,X0TEMP,YH,Z0TEMP,INLETX,YOILLN,INLETZ)
YH=YH+YOILLN
RSET(B,I:ii:OSA6,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
RSET(B,I:ii:LFA4,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:COL4,X0TEMP,YH,Z0TEMP,INLETX,YCOLLN,INLETZ)
YH=YH+YCOLLN
RSET(B,I:ii:UFA4,X0TEMP,YH,Z0TEMP,INLETX,YFLALN,INLETZ)
YH=YH+YFLALN
RSET(B,I:ii:OSA7,X0TEMP,YH,Z0TEMP,INLETX,YOSALN,INLETZ)
YH=YH+YOSALN
    *** OFA
RSET(B,I:ii:OFA1,X0TEMP,YH,Z0TEMP,INLETX,YOFAL1,INLETZ)
YH=YH+YOFAL1
RSET(B,I:ii:OFA2,X0TEMP,YH,Z0TEMP,INLETX,YOFAL2,INLETZ)
ENDDO

    *** Nose
RSET(B,NOSESLP,0.0,YVLAST-YNOSLN,0.0,XNOSLN,YNOSLN,ZWLAST)

    *** Set grid spacing in each region
RSET(X,1, 7,1.0);RSET(X,2, 5,1.0)
RSET(X,3, 1,1.0);RSET(X,4, 2,1.0)
RSET(X,5, 6,1.0);RSET(X,6, 7,1.0)
RSET(Y,1, 7,1.0);RSET(Y,2, 5,1.0)
RSET(Y,28,5,1.0);RSET(Y,29,7,1.0)
RSET(Z,1, 6,1.0);RSET(Z,2,10,1.0)
RSET(Z,3, 6,1.0)
  
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd

onephs=f
solve(p1,u1,u2,v1,v2,w1,w2,r1,r2,rs)
name(r1)=gas;name(r2)=fue;name(rs)=shad
solutn(p1,y,y,y,p,p,p)

    ** provide storage for inter-phase mass transfer etc.
store(mdot,cfip,rho1,prps)

    ** Solve additionally for the mixture fraction, i.e. the
       quantity of phase-2 material which has entered phase 1.

store(yco,yo2,yco2,yn2,yh2,yh2o)

IF(WWOOD) THEN
solve(c1);name(c1)=wood
solve(c3);name(c3)=fwd
solve(char)
ENDIF

solve(c5);name(c5)=mixf
store(yvol)
  
    ** solve for enthalpy & store temperature
    
solve(h1,h2);store(tmp1,tmp2)
 
    ** k-e turbulence model
    
turmod(kemodl);store(enut);kelin=1

    GROUP 8. Terms (in differential equations) & devices
    
eqdvdp=f;denpco=t;isolx=0;isoly=0;isolz=0
terms(h1,n,p,y,p,p,p);terms(h2,n,p,y,p,p,p)
if(wwood) then
terms(fwd,n,p,y,p,p,y)
terms(wood,p,p,n,p,y,n);terms(char,p,p,p,p,y,n)
endif

STORE(VPOR,EPOR,NPOR,HPOR)
STORE(PRPS)
 ************************************************************
  Group 9. Properties
REAL (CINCL,HINCL)
CP1     = CP;     CP2     = CP1
RHO1    = 7GASES; RHO2    = 1.E03
PRESS0  = 1.E05
TEMP0   = 0.0
CINCL   = 0.7075; HINCL   = 0.0552
RHO1A   = CINCL;  RHO1B   = HINCL

HINCOAL=CP2*TPIPE+CINCL*32.792E6+HINCL*120.9E6
HINCOAL
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
  
CFIPS=GRND1 ; cfipc=1.e5
CMDOT=grnd3; CMDTA=COALBURN; CMDTC=gfs
CINT(MIXF)=0.0; CINT(FWD)=0.0
PHINT(MIXF)=1.0;PHINT(FWD)=0.0
PHINT(WOOD)=0.0;PHINT(CHAR)=0.0
PHINT(H1) = 7GASES; PHINT(H2) = 7GASES
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
FIINIT(VPOR) =  1.000E+00 ;FIINIT(EPOR) =  1.000E+00
FIINIT(NPOR) =  1.000E+00 ;FIINIT(HPOR) =  1.000E+00
fiinit(u1)=0.0;fiinit(v1)=3.0;fiinit(w1)=0.0
fiinit(u2)=0.0;fiinit(v2)=3.0;fiinit(w2)=0.0

FIINIT(MIXF) = 0.1
FIINIT(GAS)=.9999;FIINIT(FUE)=1.E-04;FIINIT(SHAD)=1.E-04
FIINIT(MDOT)=0.01

fiinit(h1)=hinwbx;fiinit(h2)=hinCOAL
fiinit(rho1)=rhowbx

   *  Assumes 5% turbulence intensity
REAL(KEINIT,EPINIT)
KEINIT=0.0025*10.*10.
fiinit(ke)=KEINIT
epinit=0.1643*keinit**1.5/YOSALN/10.
fiinit(ep)= epinit

INIADD  =    F

IF (SLOPES) THEN
     ** Settings for sloped hopper and furnace nose
     *  Front Half of Bottom Slope
sltypF=4
pnt1x=0.0; pnt1y=yslpLN
PATCH(slp1,inival,#1,#2,#1,#1,#1,nz,1,1)
coval(slp1,vpor,0.0,grnd)
coval(slp1,epor,0.0,grnd)
coval(slp1,npor,0.0,grnd)
coval(slp1,hpor,0.0,grnd)
SPEDAT(SET,FRNTSLP,TYPE,I,4)
SPEDAT(SET,FRNTSLP,ANGLE,R,SLPANG/DG2RD)
SPEDAT(SET,FRNTSLP,PNTX,R,0.0)
SPEDAT(SET,FRNTSLP,PNTY,R,YSLPLN)
   *  Back Half of Bottom Slope
sltypB=2
pnt2x=xULAST-xslpLN; pnt2y=0.0
PATCH(slp2,inival,#5,#6,#1,#1,#1,nz,1,1)
coval(slp2,vpor,0.0,grnd)
coval(slp2,epor,0.0,grnd)
coval(slp2,npor,0.0,grnd)
coval(slp2,hpor,0.0,grnd)
SPEDAT(SET,BACKSLP,TYPE,I,2)
SPEDAT(SET,BACKSLP,ANGLE,R,SLPANG/DG2RD)
SPEDAT(SET,BACKSLP,PNTX,R,xULAST-xslpLN)
SPEDAT(SET,BACKSLP,PNTY,R,0.0)
   *  Nose
sltypN=2
pnt3x=0.0; pnt3y=YVLAST-ynosLN
PATCH(slp3,inival,#1,#4,#29,#29,#1,nz,1,1)
coval(slp3,vpor,0.0,grnd)
coval(slp3,epor,0.0,grnd)
coval(slp3,npor,0.0,grnd)
coval(slp3,hpor,0.0,grnd)
SPEDAT(SET,NOSESLP,TYPE,I,2)
SPEDAT(SET,NOSESLP,ANGLE,R,NOSANG/DG2RD)
SPEDAT(SET,NOSESLP,PNTX,R,0.0)
SPEDAT(SET,NOSESLP,PNTY,R,YVLAST-ynosLN)
ENDIF

     *** Patches to block behind windboxes ***
        
      ** south-west corner
      
PATCH(WBX1B1,INIVAL,1,1,#3,#27,6,6,1,1)
INIT(WBX1B1,PRPS,0.0,198.0)
coval(WBX1B1,npor,0.0,0.0)
coval(WBX1B1,hpor,0.0,0.0)
coval(WBX1B1,epor,0.0,0.0)
PATCH(WBX1B2,INIVAL,1,7,#3,#27,1,1,1,1)
INIT(WBX1B2,PRPS,0.0,198.0)
coval(WBX1B2,npor,0.0,0.0)
coval(WBX1B2,hpor,0.0,0.0)
coval(WBX1B2,epor,0.0,0.0)
PATCH(WBX1B3,INIVAL,1,6,#3,#27,2,2,1,1)
INIT(WBX1B3,PRPS,0.0,198.0)
coval(WBX1B3,npor,0.0,0.0)
coval(WBX1B3,hpor,0.0,0.0)
coval(WBX1B3,epor,0.0,0.0)
PATCH(WBX1B4,INIVAL,5,5,#3,#27,3,3,1,1)
INIT(WBX1B4,PRPS,0.0,198.0)
coval(WBX1B4,npor,0.0,0.0)
coval(WBX1B4,hpor,0.0,0.0)
coval(WBX1B4,epor,0.0,0.0)
PATCH(WBX1B9,INIVAL,1,1,#3,#27,3,3,1,1)
INIT(WBX1B9,PRPS,0.0,198.0)

PATCH(WBX1B5,INIVAL,1,1,%2,%2,1,5,1,1)
coval(WBX1B5,npor,0.0,0.0)
PATCH(WBX1B6,INIVAL,1,1,%27,%27,1,5,1,1)
coval(WBX1B6,npor,0.0,0.0)
PATCH(WBX1B7,INIVAL,2,4,%2,%2,1,3,1,1)
coval(WBX1B7,npor,0.0,0.0)
PATCH(WBX1B8,INIVAL,2,4,%27,%27,1,3,1,1)
coval(WBX1B8,npor,0.0,0.0)

      ** south-east corner
      
PATCH(WBX2B1,INIVAL,1,1,#3,#27,NZ-5,(NZ-5),1,1)
INIT(WBX2B1,PRPS,0.0,198.0)
coval(WBX2B1,npor,0.0,0.0)
coval(WBX2B1,hpor,0.0,0.0)
coval(WBX2B1,epor,0.0,0.0)
PATCH(WBX2B2,INIVAL,1,7,#3,#27,NZ,NZ,1,1)
INIT(WBX2B2,PRPS,0.0,198.0)
coval(WBX2B2,npor,0.0,0.0)
coval(WBX2B2,hpor,0.0,0.0)
coval(WBX2B2,epor,0.0,0.0)
PATCH(WBX2B3,INIVAL,1,6,#3,#27,NZ-1,NZ-1,1,1)
INIT(WBX2B3,PRPS,0.0,198.0)
coval(WBX2B3,npor,0.0,0.0)
coval(WBX2B3,hpor,0.0,0.0)
coval(WBX2B3,epor,0.0,0.0)
PATCH(WBX2B4,INIVAL,5,5,#3,#27,NZ-2,NZ-2,1,1)
INIT(WBX2B4,PRPS,0.0,198.0)
coval(WBX2B4,npor,0.0,0.0)
coval(WBX2B4,hpor,0.0,0.0)
coval(WBX2B4,epor,0.0,0.0)
PATCH(WBX2B9,INIVAL,1,1,#3,#27,NZ-2,NZ-2,1,1)
INIT(WBX2B9,PRPS,0.0,198.0)

PATCH(WBX2B5,INIVAL,1,1,%2,%2,NZ-4,NZ,1,1)
coval(WBX2B5,npor,0.0,0.0)
PATCH(WBX2B6,INIVAL,1,1,%27,%27,NZ-4,NZ,1,1)
coval(WBX2B6,npor,0.0,0.0)
PATCH(WBX2B7,INIVAL,2,4,%2,%2,NZ-2,NZ,1,1)
coval(WBX2B7,npor,0.0,0.0)
PATCH(WBX2B8,INIVAL,2,4,%27,%27,NZ-2,NZ,1,1)
coval(WBX2B8,npor,0.0,0.0)

      ** north-east corner
      
PATCH(WBX3B1,INIVAL,NX,NX,#3,#27,NZ-5,NZ-5,1,1)
INIT(WBX3B1,PRPS,0.0,198.0)
coval(WBX3B1,npor,0.0,0.0)
coval(WBX3B1,hpor,0.0,0.0)
coval(WBX3B1,epor,0.0,0.0)
PATCH(WBX3B2,INIVAL,NX-6,NX  ,#3,#27,NZ,NZ,1,1)
INIT(WBX3B2,PRPS,0.0,198.0)
coval(WBX3B2,npor,0.0,0.0)
coval(WBX3B2,hpor,0.0,0.0)
coval(WBX3B2,epor,0.0,0.0)
PATCH(WBX3B3,INIVAL,NX-5,NX  ,#3,#27,NZ-1,NZ-1,1,1)
INIT(WBX3B3,PRPS,0.0,198.0)
coval(WBX3B3,npor,0.0,0.0)
coval(WBX3B3,hpor,0.0,0.0)
coval(WBX3B3,epor,0.0,0.0)
PATCH(WBX3B4,INIVAL,NX-4,NX-4,#3,#27,NZ-2,NZ-2,1,1)
INIT(WBX3B4,PRPS,0.0,198.0)
coval(WBX3B4,npor,0.0,0.0)
coval(WBX3B4,hpor,0.0,0.0)
coval(WBX3B4,epor,0.0,0.0)
PATCH(WBX3B9,INIVAL,NX,NX,#3,#27,NZ-2,NZ-2,1,1)
INIT(WBX3B9,PRPS,0.0,198.0)

PATCH(WBX3B5,INIVAL,NX,NX,%2,%2,NZ-4,NZ,1,1)
coval(WBX3B5,npor,0.0,0.0)
PATCH(WBX3B6,INIVAL,NX,NX,%27,%27,NZ-4,NZ,1,1)
coval(WBX3B6,npor,0.0,0.0)
PATCH(WBX3B7,INIVAL,NX-3,NX-1,%2,%2,NZ-2,NZ,1,1)
coval(WBX3B7,npor,0.0,0.0)
PATCH(WBX3B8,INIVAL,NX-3,NX-1,%27,%27,NZ-2,NZ,1,1)
coval(WBX3B8,npor,0.0,0.0)

      ** north-west corner
      
PATCH(WBX4B1,INIVAL,NX,NX,#3,#27,6,6,1,1)
INIT(WBX4B1,PRPS,0.0,198.0)
coval(WBX4B1,npor,0.0,0.0)
coval(WBX4B1,hpor,0.0,0.0)
coval(WBX4B1,epor,0.0,0.0)
PATCH(WBX4B2,INIVAL,NX-6,NX  ,#3,#27,1,1,1,1)
INIT(WBX4B2,PRPS,0.0,198.0)
coval(WBX4B2,npor,0.0,0.0)
coval(WBX4B2,hpor,0.0,0.0)
coval(WBX4B2,epor,0.0,0.0)
PATCH(WBX4B3,INIVAL,NX-5,NX  ,#3,#27,2,2,1,1)
INIT(WBX4B3,PRPS,0.0,198.0)
coval(WBX4B3,npor,0.0,0.0)
coval(WBX4B3,hpor,0.0,0.0)
coval(WBX4B3,epor,0.0,0.0)
PATCH(WBX4B4,INIVAL,NX-4,NX-4,#3,#27,3,3,1,1)
INIT(WBX4B4,PRPS,0.0,198.0)
coval(WBX4B4,npor,0.0,0.0)
coval(WBX4B4,hpor,0.0,0.0)
coval(WBX4B4,epor,0.0,0.0)
PATCH(WBX4B9,INIVAL,NX,NX,#3,#27,3,3,1,1)
INIT(WBX4B9,PRPS,0.0,198.0)
PATCH(WBX4B5,INIVAL,NX,NX,%2,%2,1,5,1,1)
coval(WBX4B5,npor,0.0,0.0)
PATCH(WBX4B6,INIVAL,NX,NX,%27,%27,1,5,1,1)
coval(WBX4B6,npor,0.0,0.0)
PATCH(WBX4B7,INIVAL,NX-3,NX-1,%2,%2,1,3,1,1)
coval(WBX4B7,npor,0.0,0.0)
PATCH(WBX4B8,INIVAL,NX-3,NX-1,%27,%27,1,3,1,1)
coval(WBX4B8,npor,0.0,0.0)

 ************************************************************
  Group 12. Convection and diffusion adjustments
 ************************************************************
  Group 13. Boundary & Special Sources

    ************** Wall patches ***********
    --- Bottom water
WALL(W01BOT,SOUTH,#3,#4,#1, #1, #1,#3,#1,#NREGT)
    --- Back slope
WALL(W02WLO,WEST, #1,#1,#1, #1, #1,#3,#1,#NREGT)
    --- Between back slope and windbox
WALL(W03WLO,WEST, #1,#1,#2, #2, #1,#3,#1,#NREGT)
    --- Front slope
WALL(W04ELO,EAST, #6,#6,#1, #1, #1,#3,#1,#NREGT)
    --- Between front slope and windbox
WALL(W05ELO,EAST, #6,#6,#2, #2, #1,#3,#1,#NREGT)
WALL(W06LLO,LOW,  #1,#6,#1, #2, #1,#1,#1,#NREGT)
WALL(W07HLO,HIGH, #1,#6,#1, #2, #3,#3,#1,#NREGT)
    --- Windbox region
WALL(W12LMD,LOW,  #2,#5,#3, #27,#1,#1,#1,#NREGT)
WALL(W13HMD,HIGH, #2,#5,#3, #27,#3,#3,#1,#NREGT)
WALL(W14WMD,WEST, #1,#1,#3, #27,#2,#2,#1,#NREGT)
WALL(W15EMD,EAST, #6,#6,#3, #27,#2,#2,#1,#NREGT)
    --- Back between windbox and nose
WALL(W08WHI,WEST, #1,#1,#28,#28,#1,#3,#1,#NREGT)
WALL(W09EHI,EAST, #6,#6,#28,#29,#1,#3,#1,#NREGT)
WALL(W10LHI,LOW,  #1,#6,#28,#29,#1,#1,#1,#NREGT)
WALL(W11HHI,HIGH, #1,#6,#28,#29,#3,#3,#1,#NREGT)
    
   * WALL friction/heat transfer boundary condition

COVAL(W01BOT,H1,grnd2,cp*300)
COVAL(W02WLO,H1,grnd2,cp*twall)
COVAL(W03WLO,H1,grnd2,cp*twall)
COVAL(W04ELO,H1,grnd2,cp*twall)
COVAL(W05ELO,H1,grnd2,cp*twall)
COVAL(W06LLO,H1,grnd2,cp*twall)
COVAL(W07HLO,H1,grnd2,cp*twall)
COVAL(W08WHI,H1,grnd2,cp*twall)
COVAL(W09EHI,H1,grnd2,cp*twall)
COVAL(W10LHI,H1,grnd2,cp*twall)
COVAL(W11HHI,H1,grnd2,cp*twall)
COVAL(W12LMD,H1,grnd2,cp*twall)
COVAL(W13HMD,H1,grnd2,cp*twall)
COVAL(W14WMD,H1,grnd2,cp*twall)
COVAL(W15EMD,H1,grnd2,cp*twall)

    ******************** COMBUSTION patches *************
REAL(XCOMP,ZCOMP)
XCOMP=0.5086
ZCOMP=0.4914

REAL(XBURN,ZBURN)
XBURN=.311
ZBURN=.301

REAL(AANG1,AANG2,AANG3,AANG4)
REAL(BANG1,BANG2,BANG3,BANG4)
REAL(a1u,a1w,a2u,a2w,a3u,a3w,a4u,a4w)
REAL(B1u,B1w,B2u,B2w,B3u,B3w,B4u,B4w)
REAL(VWBX1,VWBX2,VWBX3,VWBX4)

REAL(EPFAC1,EPFAC2,INLA1,INLA2,YAHT,DUMVEL,DUMKE,DUMEP,DUMFLX)
REAL(XRGTMP,ZRGTMP,USIGN,WSIGN)
REAL(A1X1,A1X2,A1X3,A1X4,A1X5,A1X6,A1X7,A1X8)
REAL(A2X1,A2X2,A2X3,A2X4,A2X5,A2X6,A2X7,A2X8)
REAL(A3X1,A3X2,A3X3,A3X4,A3X5,A3X6,A3X7,A3X8)
REAL(A4X1,A4X2,A4X3,A4X4,A4X5,A4X6,A4X7,A4X8)
REAL(A1Z1,A1Z2,A1Z3,A1Z4,A1Z5,A1Z6,A1Z7,A1Z8)
REAL(A2Z1,A2Z2,A2Z3,A2Z4,A2Z5,A2Z6,A2Z7,A2Z8)
REAL(A3Z1,A3Z2,A3Z3,A3Z4,A3Z5,A3Z6,A3Z7,A3Z8)
REAL(A4Z1,A4Z2,A4Z3,A4Z4,A4Z5,A4Z6,A4Z7,A4Z8)

REAL(F1X1,F1X2,F1X3,F1X4,F1X5,F1X6,F1X7,F1X8)
REAL(F2X1,F2X2,F2X3,F2X4,F2X5,F2X6,F2X7,F2X8)
REAL(F3X1,F3X2,F3X3,F3X4,F3X5,F3X6,F3X7,F3X8)
REAL(F4X1,F4X2,F4X3,F4X4,F4X5,F4X6,F4X7,F4X8)
REAL(F1Z1,F1Z2,F1Z3,F1Z4,F1Z5,F1Z6,F1Z7,F1Z8)
REAL(F2Z1,F2Z2,F2Z3,F2Z4,F2Z5,F2Z6,F2Z7,F2Z8)
REAL(F3Z1,F3Z2,F3Z3,F3Z4,F3Z5,F3Z6,F3Z7,F3Z8)
REAL(F4Z1,F4Z2,F4Z3,F4Z4,F4Z5,F4Z6,F4Z7,F4Z8)

BANG1=33.5*DG2RD
BANG2=41.5*DG2RD
BANG3=BANG1
BANG4=BANG2

B1U=SIN(BANG1);B1W=COS(BANG1)
B2U=SIN(BANG2);B2W=COS(BANG2)
B3U=SIN(BANG3);B3W=COS(BANG3)
B4U=SIN(BANG4);B4W=COS(BANG4)

AANG1=BANG1-23.0*DG2RD
AANG2=BANG2+23.0*DG2RD
AANG3=AANG1
AANG4=AANG2

A1U=SIN(AANG1);A1W=COS(AANG1)
A2U=SIN(AANG2);A2W=COS(AANG2)
A3U=SIN(AANG3);A3W=COS(AANG3)
A4U=SIN(AANG4);A4W=COS(AANG4)

real(tltbx1,tltbx2,tltbx3,tltbx4)
tltbx1=14.0;tltbx2=14.0;tltbx3=14.0;tltbx4=14.0
tltbx1=tltbx1*dg2rd;tltbx2=tltbx2*dg2rd
tltbx3=tltbx3*dg2rd;tltbx4=tltbx4*dg2rd
VWBX1=sin(tltbx1);VWBX2=sin(tltbx2)
VWBX3=sin(tltbx3);VWBX4=sin(tltbx4)

EPFAC1 = 0.1643/ZBURN
EPFAC2 = 0.1643/XBURN

    **** FUEL INLET PATCHES ****

INLA1 = YCOLLN * ZBURN
INLA2 = YCOLLN * XBURN

DO jj=1,4

CASE jj OF
WHEN 1
+ YAHT = 5
WHEN 2
+ YAHT = 11
WHEN 3
+ YAHT = 17
WHEN 4
+ YAHT = 23
ENDCASE

  * Coal Flows

FRAC=XCOMP/INLA1

  * Primary Air Flows

b1x:jj:=    b1p:jj:*XCOMP
b2x:jj:=    b2p:jj:*XCOMP
b3x:jj:=    b3p:jj:*XCOMP
b4x:jj:=    b4p:jj:*XCOMP

b1z:jj:=    b1p:jj:*ZCOMP
b2z:jj:=    b2p:jj:*ZCOMP
b3z:jj:=    b3p:jj:*ZCOMP
b4z:jj:=    b4p:jj:*ZCOMP

patch(B1X:jj:,WEST,2,2,#YAHT,#YAHT,4,5,1,1)
DUMFLX=B1X:jj:/INLA1
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(B1X:jj:,P1,FIXFLU, DUMFLX)
COVAL(B1X:jj:,P2,FIXFLU, C1F:jj:*FRAC)
COVAL(B1X:jj:,U1,onlyms, DUMVEL*b1u)
COVAL(B1X:jj:,V1,onlyms, DUMVEL*VWBX1)
COVAL(B1X:jj:,W1,ONLYMS, DUMVEL*b1w)
COVAL(B1X:jj:,u2,onlyms, DUMVEL*b1u)
COVAL(B1X:jj:,V2,onlyms, DUMVEL*VWBX1)
COVAL(B1X:jj:,W2,ONLYMS, DUMVEL*b1w)
COVAL(B1X:jj:,H1,ONLYMS, HINPIPE)
COVAL(B1X:jj:,H2,ONLYMS, HINCOAL)
COVAL(B1X:jj:,KE,ONLYMS, DUMKE)
COVAL(B1X:jj:,EP,ONLYMS, DUMEP)

patch(B2X:jj:,WEST,2,2,#YAHT,#YAHT,NZ-4,NZ-3,1,1)
DUMFLX=B2X:jj:/INLA1
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(B2X:jj:,P1,FIXFLU, DUMFLX)
COVAL(B2X:jj:,P2,FIXFLU, C2F:jj:*FRAC)
COVAL(B2X:jj:,U1,onlyms, DUMVEL*b2u)
COVAL(B2X:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(B2X:jj:,W1,ONLYMS,-DUMVEL*b2w)
COVAL(B2X:jj:,u2,onlyms, DUMVEL*b2u)
COVAL(B2X:jj:,V2,ONLYMS, DUMVEL*VWBX2)
COVAL(B2X:jj:,W2,ONLYMS,-DUMVEL*b2w)
COVAL(B2X:jj:,H1,ONLYMS, HINPIPE)
COVAL(B2X:jj:,H2,ONLYMS, HINCOAL)
COVAL(B2X:jj:,KE,ONLYMS, DUMKE)
COVAL(B2X:jj:,EP,ONLYMS, DUMEP)

patch(B3X:jj:,EAST,NX-1,NX-1,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=B3X:jj:/INLA1
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(B3X:jj:,P1,FIXFLU, DUMFLX)
COVAL(B3X:jj:,P2,FIXFLU, C3F:jj:*FRAC)
COVAL(B3X:jj:,U1,onlyms,-DUMVEL*b3u)
COVAL(B3X:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(B3X:jj:,W1,ONLYMS,-DUMVEL*b3w)
COVAL(B3X:jj:,u2,onlyms,-DUMVEL*b3u)
COVAL(B3X:jj:,V2,ONLYMS, DUMVEL*VWBX3)
COVAL(B3X:jj:,W2,ONLYMS,-DUMVEL*b3w)
COVAL(B3X:jj:,H1,ONLYMS, HINPIPE)
COVAL(B3X:jj:,H2,ONLYMS, HINCOAL)
COVAL(B3X:jj:,KE,ONLYMS, DUMKE)
COVAL(B3X:jj:,EP,ONLYMS, DUMEP)

patch(B4X:jj:,EAST,NX-1,NX-1,#YAHT,#YAHT,4,5,1,1)
DUMFLX=B4X:jj:/INLA1
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(B4X:jj:,P1,FIXFLU, DUMFLX)
COVAL(B4X:jj:,P2,FIXFLU, C4F:jj:*FRAC)
COVAL(B4X:jj:,U1,onlyms,-DUMVEL*b4u)
COVAL(B4X:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(B4X:jj:,W1,ONLYMS, DUMVEL*b4w)
COVAL(B4X:jj:,u2,onlyms,-DUMVEL*B4u)
COVAL(B4X:jj:,V2,ONLYMS, DUMVEL*VWBX4)
COVAL(B4X:jj:,W2,ONLYMS, DUMVEL*b4w)
COVAL(B4X:jj:,H1,ONLYMS, HINPIPE)
COVAL(B4X:jj:,H2,ONLYMS, HINCOAL)
COVAL(B4X:jj:,KE,ONLYMS, DUMKE)
COVAL(B4X:jj:,EP,ONLYMS, DUMEP)

FRAC=ZCOMP/INLA1

patch(B1Z:jj:,LOW,2,4,#YAHT,#YAHT,4,4,1,1)
DUMFLX=B1Z:jj:/INLA2
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(B1Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(B1Z:jj:,P2,FIXFLU, C1F:jj:*FRAC)
COVAL(B1Z:jj:,U1,ONLYMS, DUMVEL*b1u)
COVAL(B1Z:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(B1Z:jj:,W1,onlyms, DUMVEL*b1w)
COVAL(B1Z:jj:,u2,ONLYMS, DUMVEL*B1u)
COVAL(B1Z:jj:,V2,ONLYMS, DUMVEL*VWBX1)
COVAL(B1Z:jj:,W2,onlyms, DUMVEL*b1w)
COVAL(B1Z:jj:,H1,ONLYMS, HINPIPE)
COVAL(B1Z:jj:,H2,ONLYMS, HINCOAL)
COVAL(B1Z:jj:,KE,ONLYMS, DUMKE)
COVAL(B1Z:jj:,EP,ONLYMS, DUMEP)

patch(B2Z:jj:,HIGH,2,4,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=B2Z:jj:/INLA2
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(B2Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(B2Z:jj:,P2,FIXFLU, C2F:jj:*FRAC)
COVAL(B2Z:jj:,U1,ONLYMS, DUMVEL*b2u)
COVAL(B2Z:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(B2Z:jj:,W1,onlyms,-DUMVEL*b2w)
COVAL(B2Z:jj:,u2,ONLYMS, DUMVEL*B2u)
COVAL(B2Z:jj:,V2,ONLYMS, DUMVEL*VWBX2)
COVAL(B2Z:jj:,W2,onlyms,-DUMVEL*b2w)
COVAL(B2Z:jj:,H1,ONLYMS, HINPIPE)
COVAL(B2Z:jj:,H2,ONLYMS, HINCOAL)
COVAL(B2Z:jj:,KE,ONLYMS, DUMKE)
COVAL(B2Z:jj:,EP,ONLYMS, DUMEP)
        
patch(B3Z:jj:,HIGH,nx-3,nx-1,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=B3Z:jj:/INLA2
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(B3Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(B3Z:jj:,P2,FIXFLU, C3F:jj:*FRAC)
COVAL(B3Z:jj:,U1,ONLYMS,-DUMVEL*b3u)
COVAL(B3Z:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(B3Z:jj:,W1,onlyms,-DUMVEL*b3w)
COVAL(B3Z:jj:,u2,ONLYMS,-DUMVEL*B3u)
COVAL(B3Z:jj:,V2,ONLYMS, DUMVEL*VWBX3)
COVAL(B3Z:jj:,W2,onlyms,-DUMVEL*b3w)
COVAL(B3Z:jj:,H1,ONLYMS, HINPIPE)
COVAL(B3Z:jj:,H2,ONLYMS, HINCOAL)
COVAL(B3Z:jj:,KE,ONLYMS, DUMKE)
COVAL(B3Z:jj:,EP,ONLYMS, DUMEP)

patch(B4Z:jj:,LOW,nx-3,nx-1,#YAHT,#YAHT,4,4,1,1)
DUMFLX=B4Z:jj:/INLA2
DUMVEL=DUMFLX/RHOPIPE
DUMKE=KEPA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(B4Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(B4Z:jj:,P2,FIXFLU, C4F:jj:*FRAC)
COVAL(B4Z:jj:,U1,ONLYMS,-DUMVEL*b4u)
COVAL(B4Z:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(B4Z:jj:,W1,onlyms, DUMVEL*b4w)
COVAL(B4Z:jj:,u2,ONLYMS,-DUMVEL*B4u)
COVAL(B4Z:jj:,V2,ONLYMS, DUMVEL*VWBX4)
COVAL(B4Z:jj:,W2,onlyms, DUMVEL*b4w)
COVAL(B4Z:jj:,H1,ONLYMS, HINPIPE)
COVAL(B4Z:jj:,H2,ONLYMS, HINCOAL)
COVAL(B4Z:jj:,KE,ONLYMS, DUMKE)
COVAL(B4Z:jj:,EP,ONLYMS, DUMEP)

ENDDO

   **** FUEL AIR INLET PATCHES ****

INLA1 = YFLALN * ZBURN
INLA2 = YFLALN * XBURN

DO jj=1,8

CASE jj OF
WHEN 1
+ YAHT = 4
WHEN 2
+ YAHT = 6
WHEN 3
+ YAHT = 10
WHEN 4
+ YAHT = 12
WHEN 5
+ YAHT = 16
WHEN 6
+ YAHT = 18
WHEN 7
+ YAHT = 22
WHEN 8
+ YAHT = 24
ENDCASE

f1x:jj:=f1a:jj:*XCOMP
f2x:jj:=f2a:jj:*XCOMP
f3x:jj:=f3a:jj:*XCOMP
f4x:jj:=f4a:jj:*XCOMP

f1z:jj:=f1a:jj:*ZCOMP
f2z:jj:=f2a:jj:*ZCOMP
f3z:jj:=f3a:jj:*ZCOMP
f4z:jj:=f4a:jj:*ZCOMP

patch(F1X:jj:,WEST,2,2,#YAHT,#YAHT,4,5,1,1)
DUMFLX=F1X:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(F1X:jj:,P1,FIXFLU, DUMFLX)
COVAL(F1X:jj:,U1,onlyms, DUMVEL*b1u)
COVAL(F1X:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(F1X:jj:,W1,ONLYMS, DUMVEL*b1w)
COVAL(F1X:jj:,H1,ONLYMS, HINWBX)
COVAL(F1X:jj:,KE,ONLYMS, DUMKE)
COVAL(F1X:jj:,EP,ONLYMS, DUMEP)

patch(F2X:jj:,WEST,2,2,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=F2X:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(F2X:jj:,P1,FIXFLU, DUMFLX)
COVAL(F2X:jj:,U1,onlyms, DUMVEL*b2u)
COVAL(F2X:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(F2X:jj:,W1,ONLYMS,-DUMVEL*b2w)
COVAL(F2X:jj:,H1,ONLYMS, HINWBX)
COVAL(F2X:jj:,KE,ONLYMS, DUMKE)
COVAL(F2X:jj:,EP,ONLYMS, DUMEP)

patch(F3X:jj:,east,NX-1,NX-1,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=F3X:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(F3X:jj:,P1,FIXFLU, DUMFLX)
COVAL(F3X:jj:,U1,onlyms,-DUMVEL*b3u)
COVAL(F3X:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(F3X:jj:,W1,ONLYMS,-DUMVEL*b3w)
COVAL(F3X:jj:,H1,ONLYMS, HINWBX)
COVAL(F3X:jj:,KE,ONLYMS, DUMKE)
COVAL(F3X:jj:,EP,ONLYMS, DUMEP)

patch(F4X:jj:,EAST,NX-1,NX-1,#YAHT,#YAHT,4,5,1,1)
DUMFLX=F4X:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(F4X:jj:,P1,FIXFLU, DUMFLX)
COVAL(F4X:jj:,U1,onlyms,-DUMVEL*b4u)
COVAL(F4X:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(F4X:jj:,W1,ONLYMS, DUMVEL*b4w)
COVAL(F4X:jj:,H1,ONLYMS, HINWBX)
COVAL(F4X:jj:,KE,ONLYMS, DUMKE)
COVAL(F4X:jj:,EP,ONLYMS, DUMEP)

patch(F1Z:jj:,LOW,2,4,#YAHT,#YAHT,4,4,1,1)
DUMFLX=F1Z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(F1Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(F1Z:jj:,U1,ONLYMS, DUMVEL*b1u)
COVAL(F1Z:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(F1Z:jj:,W1,onlyms, DUMVEL*b1w)
COVAL(F1Z:jj:,H1,ONLYMS, HINWBX)
COVAL(F1Z:jj:,KE,ONLYMS, DUMKE)
COVAL(F1Z:jj:,EP,ONLYMS, DUMEP)

patch(F2Z:jj:,HIGH,2,4,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=F2Z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(F2Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(F2Z:jj:,U1,ONLYMS, DUMVEL*b2u)
COVAL(F2Z:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(F2Z:jj:,W1,onlyms,-DUMVEL*b2w)
COVAL(F2Z:jj:,H1,ONLYMS, HINWBX)
COVAL(F2Z:jj:,KE,ONLYMS, DUMKE)
COVAL(F2Z:jj:,EP,ONLYMS, DUMEP)

patch(F3Z:jj:,HIGH,nx-3,nx-1,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=F3Z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(F3Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(F3Z:jj:,U1,ONLYMS,-DUMVEL*b3u)
COVAL(F3Z:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(F3Z:jj:,W1,onlyms,-DUMVEL*b3w)
COVAL(F3Z:jj:,H1,ONLYMS, HINWBX)
COVAL(F3Z:jj:,KE,ONLYMS, DUMKE)
COVAL(F3Z:jj:,EP,ONLYMS, DUMEP)

patch(F4Z:jj:,LOW,nx-3,nx-1,#YAHT,#YAHT,4,4,1,1)
DUMFLX=F4Z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEFA**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(F4Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(F4Z:jj:,U1,ONLYMS,-DUMVEL*b4u)
COVAL(F4Z:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(F4Z:jj:,W1,onlyms, DUMVEL*b4w)
COVAL(F4Z:jj:,H1,ONLYMS, HINWBX)
COVAL(F4Z:jj:,KE,ONLYMS, DUMKE)
COVAL(F4Z:jj:,EP,ONLYMS, DUMEP)

ENDDO

   **** Auxiliary AIR INLET PATCHES ****

DO jj=1,8

CASE jj OF
WHEN 1
+ YAHT = 3
+ INLA1 = XBURN * yBaxLN
+ INLA2 = ZBURN * yBaxLN
WHEN 2
+ YAHT = 7
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
WHEN 3
+ YAHT = 9
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
WHEN 4
+ YAHT = 13
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
WHEN 5
+ YAHT = 15
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
WHEN 6
+ YAHT = 19
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
WHEN 7
+ YAHT = 21
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
WHEN 8
+ YAHT = 25
+ INLA1 = XBURN * YOSALN
+ INLA2 = ZBURN * YOSALN
ENDCASE

a1x:jj:=a1a:jj:*XCOMP
a2x:jj:=a2a:jj:*XCOMP
a3x:jj:=a3a:jj:*XCOMP
a4x:jj:=a4a:jj:*XCOMP

a1z:jj:=a1a:jj:*ZCOMP
a2z:jj:=a2a:jj:*ZCOMP
a3z:jj:=a3a:jj:*ZCOMP
a4z:jj:=a4a:jj:*ZCOMP

patch(A1X:jj:,WEST,2,2,#YAHT,#YAHT,4,5,1,1)
DUMFLX=A1x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(A1X:jj:,P1,FIXFLU, DUMFLX)
COVAL(A1X:jj:,U1,onlyms, DUMVEL*a1u)
COVAL(A1X:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(A1X:jj:,W1,ONLYMS, DUMVEL*a1w)
COVAL(A1X:jj:,H1,ONLYMS, HINWBX)
COVAL(A1X:jj:,KE,ONLYMS, DUMKE)
COVAL(A1X:jj:,EP,ONLYMS, DUMEP)

patch(A2X:jj:,WEST,2,2,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=A2x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(A2X:jj:,P1,FIXFLU, DUMFLX)
COVAL(A2X:jj:,U1,onlyms, DUMVEL*a2u)
COVAL(A2X:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(A2X:jj:,W1,ONLYMS,-DUMVEL*a2w)
COVAL(A2X:jj:,H1,ONLYMS, HINWBX)
COVAL(A2X:jj:,KE,ONLYMS, DUMKE)
COVAL(A2X:jj:,EP,ONLYMS, DUMEP)

patch(A3X:jj:,east,NX-1,NX-1,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=A3x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(A3X:jj:,P1,FIXFLU, DUMFLX)
COVAL(A3X:jj:,U1,onlyms,-DUMVEL*a3u)
COVAL(A3X:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(A3X:jj:,W1,ONLYMS,-DUMVEL*a3w)
COVAL(A3X:jj:,H1,ONLYMS, HINWBX)
COVAL(A3X:jj:,KE,ONLYMS, DUMKE)
COVAL(A3X:jj:,EP,ONLYMS, DUMEP)

patch(A4X:jj:,EAST,NX-1,NX-1,#YAHT,#YAHT,4,5,1,1)
DUMFLX=A4x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(A4X:jj:,P1,FIXFLU, DUMFLX)
COVAL(A4X:jj:,U1,onlyms,-DUMVEL*a4u)
COVAL(A4X:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(A4X:jj:,W1,ONLYMS, DUMVEL*a4w)
COVAL(A4X:jj:,H1,ONLYMS, HINWBX)
COVAL(A4X:jj:,KE,ONLYMS, DUMKE)
COVAL(A4X:jj:,EP,ONLYMS, DUMEP)

patch(A1Z:jj:,LOW,2,4,#YAHT,#YAHT,4,4,1,1)
DUMFLX=A1z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(A1Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(A1Z:jj:,U1,ONLYMS, DUMVEL*a1u)
COVAL(A1Z:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(A1Z:jj:,W1,onlyms, DUMVEL*a1w)
COVAL(A1Z:jj:,H1,ONLYMS, HINWBX)
COVAL(A1Z:jj:,KE,ONLYMS, DUMKE)
COVAL(A1Z:jj:,EP,ONLYMS, DUMEP)

patch(A2Z:jj:,HIGH,2,4,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=A2z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(A2Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(A2Z:jj:,U1,ONLYMS, DUMVEL*a2u)
COVAL(A2Z:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(A2Z:jj:,W1,onlyms,-DUMVEL*a2w)
COVAL(A2Z:jj:,H1,ONLYMS, HINWBX)
COVAL(A2Z:jj:,KE,ONLYMS, DUMKE)
COVAL(A2Z:jj:,EP,ONLYMS, DUMEP)

patch(A3Z:jj:,HIGH,nx-3,nx-1,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=A3z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(A3Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(A3Z:jj:,U1,ONLYMS,-DUMVEL*a3u)
COVAL(A3Z:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(A3Z:jj:,W1,onlyms,-DUMVEL*a3w)
COVAL(A3Z:jj:,H1,ONLYMS, HINWBX)
COVAL(A3Z:jj:,KE,ONLYMS, DUMKE)
COVAL(A3Z:jj:,EP,ONLYMS, DUMEP)

patch(A4Z:jj:,LOW,nx-3,nx-1,#YAHT,#YAHT,4,4,1,1)
DUMFLX=A4z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(A4Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(A4Z:jj:,U1,ONLYMS,-DUMVEL*a4u)
COVAL(A4Z:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(A4Z:jj:,W1,onlyms, DUMVEL*a4w)
COVAL(A4Z:jj:,H1,ONLYMS, HINWBX)
COVAL(A4Z:jj:,KE,ONLYMS, DUMKE)
COVAL(A4Z:jj:,EP,ONLYMS, DUMEP)

ENDDO

   **** OFA INLET PATCHES ****

DO jj=1,2

CASE jj OF
WHEN 1
+ YAHT = 26
+ INLA1 = ZBURN * YOFAL1
+ INLA2 = XBURN * YOFAL1
WHEN 2
+ YAHT = 27
+ INLA1 = ZBURN * YOFAL2
+ INLA2 = XBURN * YOFAL2
ENDCASE

a1x:jj:=a1OFA:jj:*XCOMP
a2x:jj:=a2OFA:jj:*XCOMP
a3x:jj:=a3OFA:jj:*XCOMP
a4x:jj:=a4OFA:jj:*XCOMP

a1z:jj:=a1OFA:jj:*ZCOMP
a2z:jj:=a2OFA:jj:*ZCOMP
a3z:jj:=a3OFA:jj:*ZCOMP
a4z:jj:=a4OFA:jj:*ZCOMP

patch(O1X:jj:,WEST,2,2,#YAHT,#YAHT,4,5,1,1)
DUMFLX=A1x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(O1X:jj:,P1,FIXFLU, DUMFLX)
COVAL(O1X:jj:,U1,onlyms, DUMVEL*a1u)
COVAL(O1X:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(O1X:jj:,W1,ONLYMS, DUMVEL*a1w)
COVAL(O1X:jj:,H1,ONLYMS, HINWBX)
COVAL(O1X:jj:,KE,ONLYMS, DUMKE)
COVAL(O1X:jj:,EP,ONLYMS, DUMEP)

patch(O2X:jj:,WEST,2,2,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=A2x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(O2X:jj:,P1,FIXFLU, DUMFLX)
COVAL(O2X:jj:,U1,onlyms, DUMVEL*a2u)
COVAL(O2X:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(O2X:jj:,W1,ONLYMS,-DUMVEL*a2w)
COVAL(O2X:jj:,H1,ONLYMS, HINWBX)
COVAL(O2X:jj:,KE,ONLYMS, DUMKE)
COVAL(O2X:jj:,EP,ONLYMS, DUMEP)

patch(O3X:jj:,east,NX-1,NX-1,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=A3x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(O3X:jj:,P1,FIXFLU, DUMFLX)
COVAL(O3X:jj:,U1,onlyms,-DUMVEL*a3u)
COVAL(O3X:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(O3X:jj:,W1,ONLYMS,-DUMVEL*a3w)
COVAL(O3X:jj:,H1,ONLYMS, HINWBX)
COVAL(O3X:jj:,KE,ONLYMS, DUMKE)
COVAL(O3X:jj:,EP,ONLYMS, DUMEP)

patch(O4X:jj:,EAST,NX-1,NX-1,#YAHT,#YAHT,4,5,1,1)
DUMFLX=A4x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(O4X:jj:,P1,FIXFLU, DUMFLX)
COVAL(O4X:jj:,U1,onlyms,-DUMVEL*a4u)
COVAL(O4X:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(O4X:jj:,W1,ONLYMS, DUMVEL*a4w)
COVAL(O4X:jj:,H1,ONLYMS, HINWBX)
COVAL(O4X:jj:,KE,ONLYMS, DUMKE)
COVAL(O4X:jj:,EP,ONLYMS, DUMEP)

patch(O1Z:jj:,LOW,2,4,#YAHT,#YAHT,4,4,1,1)
DUMFLX=A1z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(O1Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(O1Z:jj:,U1,ONLYMS, DUMVEL*a1u)
COVAL(O1Z:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(O1Z:jj:,W1,onlyms, DUMVEL*a1w)
COVAL(O1Z:jj:,H1,ONLYMS, HINWBX)
COVAL(O1Z:jj:,KE,ONLYMS, DUMKE)
COVAL(O1Z:jj:,EP,ONLYMS, DUMEP)

patch(O2Z:jj:,HIGH,2,4,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=A2z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(O2Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(O2Z:jj:,U1,ONLYMS, DUMVEL*a2u)
COVAL(O2Z:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(O2Z:jj:,W1,onlyms,-DUMVEL*a2w)
COVAL(O2Z:jj:,H1,ONLYMS, HINWBX)
COVAL(O2Z:jj:,KE,ONLYMS, DUMKE)
COVAL(O2Z:jj:,EP,ONLYMS, DUMEP)

patch(O3Z:jj:,HIGH,nx-3,nx-1,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=A3z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(O3Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(O3Z:jj:,U1,ONLYMS,-DUMVEL*a3u)
COVAL(O3Z:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(O3Z:jj:,W1,onlyms,-DUMVEL*a3w)
COVAL(O3Z:jj:,H1,ONLYMS, HINWBX)
COVAL(O3Z:jj:,KE,ONLYMS, DUMKE)
COVAL(O3Z:jj:,EP,ONLYMS, DUMEP)

patch(O4Z:jj:,LOW,nx-3,nx-1,#YAHT,#YAHT,4,4,1,1)
DUMFLX=A4z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(O4Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(O4Z:jj:,U1,ONLYMS,-DUMVEL*a4u)
COVAL(O4Z:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(O4Z:jj:,W1,onlyms, DUMVEL*a4w)
COVAL(O4Z:jj:,H1,ONLYMS, HINWBX)
COVAL(O4Z:jj:,KE,ONLYMS, DUMKE)
COVAL(O4Z:jj:,EP,ONLYMS, DUMEP)

ENDDO
   
   **** AIR LEAKAGE THROUGH LOWER AND UPPER OIL GUNS ****

DO jj=1,3,2

CASE jj OF
WHEN 1
+ YAHT = 8
+ INLA1 = ZBURN * YOILLN
+ INLA2 = XBURN * YOILLN
WHEN 3
+ YAHT = 20
+ INLA1 = ZBURN * YOILLN
+ INLA2 = XBURN * YOILLN
ENDCASE

a1x:jj:=a1O:jj:*XCOMP
a2x:jj:=a2O:jj:*XCOMP
a3x:jj:=a3O:jj:*XCOMP
a4x:jj:=a4O:jj:*XCOMP

a1z:jj:=a1O:jj:*ZCOMP
a2z:jj:=a2O:jj:*ZCOMP
a3z:jj:=a3O:jj:*ZCOMP
a4z:jj:=a4O:jj:*ZCOMP

patch(OIL1X:JJ:,WEST,2,2,#YAHT,#YAHT,4,5,1,1)
DUMFLX=A1x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(OIL1X:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL1X:jj:,U1,onlyms, DUMVEL*B1u)
COVAL(OIL1X:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(OIL1X:jj:,W1,ONLYMS, DUMVEL*B1w)
COVAL(OIL1X:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL1X:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL1X:jj:,EP,ONLYMS, DUMEP)

patch(OIL2X:JJ:,WEST,2,2,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=A2x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(OIL2X:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL2X:jj:,U1,onlyms, DUMVEL*B2u)
COVAL(OIL2X:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(OIL2X:jj:,W1,ONLYMS,-DUMVEL*B2w)
COVAL(OIL2X:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL2X:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL2X:jj:,EP,ONLYMS, DUMEP)

patch(OIL3X:JJ:,east,NX-1,NX-1,#YAHT,#YAHT,nz-4,nz-3,1,1)
DUMFLX=A3x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(OIL3X:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL3X:jj:,U1,onlyms,-DUMVEL*B3u)
COVAL(OIL3X:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(OIL3X:jj:,W1,ONLYMS,-DUMVEL*B3w)
COVAL(OIL3X:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL3X:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL3X:jj:,EP,ONLYMS, DUMEP)

patch(OIL4X:JJ:,EAST,NX-1,NX-1,#YAHT,#YAHT,4,5,1,1)
DUMFLX=A4x:jj:/INLA1
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(OIL4X:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL4X:jj:,U1,onlyms,-DUMVEL*B4u)
COVAL(OIL4X:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(OIL4X:jj:,W1,ONLYMS, DUMVEL*B4w)
COVAL(OIL4X:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL4X:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL4X:jj:,EP,ONLYMS, DUMEP)

patch(OIL1Z:JJ:,LOW,2,4,#YAHT,#YAHT,4,4,1,1)
DUMFLX=A1z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(OIL1Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL1Z:jj:,U1,ONLYMS, DUMVEL*B1u)
COVAL(OIL1Z:jj:,V1,ONLYMS, DUMVEL*VWBX1)
COVAL(OIL1Z:jj:,W1,onlyms, DUMVEL*B1w)
COVAL(OIL1Z:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL1Z:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL1Z:jj:,EP,ONLYMS, DUMEP)

patch(OIL2Z:JJ:,HIGH,2,4,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=A2z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(OIL2Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL2Z:jj:,U1,ONLYMS, DUMVEL*B2u)
COVAL(OIL2Z:jj:,V1,ONLYMS, DUMVEL*VWBX2)
COVAL(OIL2Z:jj:,W1,onlyms,-DUMVEL*B2w)
COVAL(OIL2Z:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL2Z:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL2Z:jj:,EP,ONLYMS, DUMEP)

patch(OIL3Z:JJ:,high,nx-3,nx-1,#YAHT,#YAHT,NZ-3,NZ-3,1,1)
DUMFLX=A3z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(OIL3Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL3Z:jj:,U1,ONLYMS,-DUMVEL*B3u)
COVAL(OIL3Z:jj:,V1,ONLYMS, DUMVEL*VWBX3)
COVAL(OIL3Z:jj:,W1,onlyms,-DUMVEL*B3w)
COVAL(OIL3Z:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL3Z:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL3Z:jj:,EP,ONLYMS, DUMEP)

patch(OIL4Z:JJ:,LOW,nx-3,nx-1,#YAHT,#YAHT,4,4,1,1)
DUMFLX=A4z:jj:/INLA2
DUMVEL=DUMFLX/RHOWBX
DUMKE=KEAX**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(OIL4Z:jj:,P1,FIXFLU, DUMFLX)
COVAL(OIL4Z:jj:,U1,ONLYMS,-DUMVEL*B4u)
COVAL(OIL4Z:jj:,V1,ONLYMS, DUMVEL*VWBX4)
COVAL(OIL4Z:jj:,W1,onlyms, DUMVEL*B4w)
COVAL(OIL4Z:jj:,H1,ONLYMS, HINWBX)
COVAL(OIL4Z:jj:,KE,ONLYMS, DUMKE)
COVAL(OIL4Z:jj:,EP,ONLYMS, DUMEP)

ENDDO

IF (WWOOD) THEN

INLA1 = ZBURN * YOILLN
INLA2 = XBURN * YOILLN

patch(WPA1X1,WEST,2,2,#14,#14,4,5,1,1)
DUMFLX=twpa*xcomp/INLA1
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(WPA1X1,P1,FIXFLU,DUMFLX)
coval(wPA1x1,u1,onlyms, dumvel*b1u)
coval(wPA1x1,v1,onlyms, dumvel*vwbx1)
coval(wPA1x1,w1,onlyms, dumvel*b1w)
COVAL(WPA1X1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA1X1,KE,ONLYMS,DUMKE)
COVAL(WPA1X1,EP,ONLYMS,DUMEP)

patch(WPA2X1,WEST,2,2,#14,#14,nz-4,nz-3,1,1)
DUMFLX=twpa*xcomp/INLA1
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(WPA2X1,P1,FIXFLU,DUMFLX)
coval(wPA2x1,u1,onlyms, dumvel*b2u)
coval(wPA2x1,v1,onlyms, dumvel*vwbx2)
coval(wPA2x1,w1,onlyms,-dumvel*b2w)
COVAL(WPA2X1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA2X1,KE,ONLYMS,DUMKE)
COVAL(WPA2X1,EP,ONLYMS,DUMEP)

patch(WPA3X1,east,NX-1,NX-1,#14,#14,nz-4,nz-3,1,1)
DUMFLX=twpa*xcomp/INLA1
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(WPA3X1,P1,FIXFLU,DUMFLX)
coval(wPA3x1,u1,onlyms,-dumvel*b3u)
coval(wPA3x1,v1,onlyms, dumvel*vwbx3)
coval(wPA3x1,w1,onlyms,-dumvel*b3w)
COVAL(WPA3X1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA3X1,KE,ONLYMS,DUMKE)
COVAL(WPA3X1,EP,ONLYMS,DUMEP)

patch(WPA4X1,EAST,NX-1,NX-1,#14,#14,4,5,1,1)
DUMFLX=twpa*xcomp/INLA1
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC1*DUMKE**1.5
COVAL(WPA4X1,P1,FIXFLU,DUMFLX)
coval(wPA4x1,u1,onlyms,-dumvel*b4u)
coval(wPA4x1,v1,onlyms, dumvel*vwbx4)
coval(wPA4x1,w1,onlyms, dumvel*b4w)
COVAL(WPA4X1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA4X1,KE,ONLYMS,DUMKE)
COVAL(WPA4X1,EP,ONLYMS,DUMEP)

patch(WPA1Z1,LOW,2,4,#14,#14,4,4,1,1)
DUMFLX=twpa*zcomp/INLA2
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(WPA1Z1,P1,FIXFLU, DUMFLX)
coval(wPA1z1,u1,onlyms, dumvel*b1u)
coval(wPA1z1,v1,onlyms, dumvel*vwbx1)
coval(wPA1z1,w1,onlyms, dumvel*b1w)
COVAL(WPA1Z1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA1Z1,KE,ONLYMS,DUMKE)
COVAL(WPA1Z1,EP,ONLYMS,DUMEP)

patch(WPA2Z1,HIGH,2,4,#14,#14,NZ-3,NZ-3,1,1)
DUMFLX=twpa*zcomp/INLA2
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(WPA2Z1,P1,FIXFLU,DUMFLX)
coval(wPA2z1,u1,onlyms, dumvel*b1u)
coval(wPA2z1,v1,onlyms, dumvel*vwbx2)
coval(wPA2z1,w1,onlyms,-dumvel*b1w)
COVAL(WPA2Z1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA2Z1,KE,ONLYMS,DUMKE)
COVAL(WPA2Z1,EP,ONLYMS,DUMEP)

patch(WPA3Z1,high,nx-3,nx-1,#14,#14,NZ-3,NZ-3,1,1)
DUMFLX=twpa*zcomp/INLA2
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(WPA3Z1,P1,FIXFLU,DUMFLX)
coval(wPA3z1,u1,onlyms,-dumvel*b1u)
coval(wPA3z1,v1,onlyms, dumvel*vwbx3)
coval(wPA3z1,w1,onlyms,-dumvel*b1w)
COVAL(WPA3Z1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA3Z1,KE,ONLYMS,DUMKE)
COVAL(WPA3Z1,EP,ONLYMS,DUMEP)

patch(WPA4Z1,LOW,nx-3,nx-1,#14,#14,4,4,1,1)
DUMFLX=twpa*zcomp/INLA2
DUMVEL=DUMFLX/RHOWa
DUMKE=KEWD**2.*DUMVEL**2.
DUMEP=EPFAC2*DUMKE**1.5
COVAL(WPA4Z1,P1,FIXFLU,DUMFLX)
coval(wPA4z1,u1,onlyms,-dumvel*b1u)
coval(wPA4z1,v1,onlyms, dumvel*vwbx4)
coval(wPA4z1,w1,onlyms, dumvel*b1w)
COVAL(WPA4Z1,H1,ONLYMS,TWOOD*CP)
COVAL(WPA4Z1,KE,ONLYMS,DUMKE)
COVAL(WPA4Z1,EP,ONLYMS,DUMEP)

patch(W1X1,WEST,2,2,#14,#14,4,5,1,1)
DUMFLX=tw*xcomp/INLA1
COVAL(W1X1,P1,FIXFLU,DUMFLX)
COVAL(W1X1,FWD,ONLYMS,1.0)
COVAL(W1X1,WOOD,ONLYMS,1.)
COVAL(W1X1,H1,ONLYMS,HINWD)

patch(W2X1,WEST,2,2,#14,#14,nz-4,nz-3,1,1)
DUMFLX=tw*xcomp/INLA1
COVAL(W2X1,P1,FIXFLU,DUMFLX)
COVAL(W2X1,FWD,ONLYMS,1.0)
COVAL(W2X1,WOOD,ONLYMS,1.)
COVAL(W2X1,H1,ONLYMS,HINWD)

patch(W3X1,east,NX-1,NX-1,#14,#14,nz-4,nz-3,1,1)
DUMFLX=tw*xcomp/INLA1
COVAL(W3X1,P1,FIXFLU,DUMFLX)
COVAL(W3X1,FWD,ONLYMS,1.0)
COVAL(W3X1,WOOD,ONLYMS,1.)
COVAL(W3X1,H1,ONLYMS,HINWD)

patch(W4X1,EAST,NX-1,NX-1,#14,#14,4,5,1,1)
DUMFLX=tw*xcomp/INLA1
COVAL(W4X1,P1,FIXFLU,DUMFLX)
COVAL(W4X1,FWD,ONLYMS,1.0)
COVAL(W4X1,WOOD,ONLYMS,1.)
COVAL(W4X1,H1,ONLYMS,HINWD)

patch(W1Z1,LOW,2,4,#14,#14,4,4,1,1)
DUMFLX=tw*zcomp/INLA2
COVAL(W1Z1,P1,FIXFLU, DUMFLX)
COVAL(W1Z1,FWD,ONLYMS,1.0)
COVAL(W1Z1,WOOD,ONLYMS,1.)
COVAL(W1Z1,H1,ONLYMS,HINWD)

patch(W2Z1,HIGH,2,4,#14,#14,NZ-3,NZ-3,1,1)
DUMFLX=tw*zcomp/INLA2
COVAL(W2Z1,P1,FIXFLU,DUMFLX)
COVAL(W2Z1,FWD,ONLYMS,1.0)
COVAL(W2Z1,WOOD,ONLYMS,1.)
COVAL(W2Z1,H1,ONLYMS,HINWD)

patch(W3Z1,high,nx-3,nx-1,#14,#14,NZ-3,NZ-3,1,1)
DUMFLX=tw*zcomp/INLA2
COVAL(W3Z1,P1,FIXFLU,DUMFLX)
COVAL(W3Z1,FWD,ONLYMS,1.0)
COVAL(W3Z1,WOOD,ONLYMS,1.)
COVAL(W3Z1,H1,ONLYMS,HINWD)

patch(W4Z1,LOW,nx-3,nx-1,#14,#14,4,4,1,1)
DUMFLX=tw*zcomp/INLA2
COVAL(W4Z1,P1,FIXFLU,DUMFLX)
COVAL(W4Z1,FWD,ONLYMS,1.0)
COVAL(W4Z1,WOOD,ONLYMS,1.)
COVAL(W4Z1,H1,ONLYMS,HINWD)

SPEDAT(SET,WOODBURN,VLINWD,R,0.5)
SPEDAT(SET,WOODBURN,HINVL,R,0.02)
SPEDAT(SET,WOODBURN,MLWTVL,R,100)
SPEDAT(SET,WOODBURN,MOINWD,R,0.2)
ENDIF

    **** Outlet Patch ****
PATCH(outlet,NORTH,#4,#6,ny,ny,1,nz,1,#NREGT)
COVAL(outlet,P1,0.1,ZERO)
COVAL(outlet,P2,0.1*rho2,ZERO)
COVAL(outlet,U1,ONLYMS,SAME)
COVAL(outlet,V1,ONLYMS,SAME)
COVAL(outlet,W1,ONLYMS,SAME)
COVAL(outlet,U2,ONLYMS,SAME)
COVAL(outlet,V2,ONLYMS,SAME)
COVAL(outlet,W2,ONLYMS,SAME)
COVAL(outlet,H1,ONLYMS,SAME)
COVAL(outlet,H2,ONLYMS,SAME)
COVAL(outlet,KE,ONLYMS,SAME)
COVAL(outlet,EP,ONLYMS,SAME)
COVAL(outlet,MIXF,ONLYMS,SAME)

IF (WWOOD) THEN
COVAL(outlet,WOOD,ONLYMS,SAME)
COVAL(outlet,CHAR,ONLYMS,SAME)
COVAL(outlet,FWD,ONLYMS,SAME)
ENDIF

  * GRAVITY boundary condition, name BUOYANCY

  *  Set gravity resolutes

BUOYA=0.;BUOYB=-9.8100E+00;BUOYC=0.
PATCH(BUOYANCY,PHASEM,#1,#NREGX,#1,#NREGY,#1,#NREGZ,#1,#NREGT)
COVAL(BUOYANCY,V1,FIXFLU,GRND2)

EGWF    =    T

 ************************************************************
  Group 14. Downstream Pressure For PARAB
 ************************************************************
  Group 15. Terminate Sweeps
LSWEEP  = 3000; TSTSWP  = -1
SELREF  = T;    RESFAC  = 1.000E-02
LITER(P1)=30
LITER(U1)=2;LITER(U2)=2
LITER(V1)=2;LITER(V2)=2
LITER(W1)=2;LITER(W2)=2
 ************************************************************
  Group 16. Terminate Iterations
 ************************************************************
  Group 17. Relaxation
relax( p1 ,linrlx,1.0)
relax(RHO1,linrlx,0.1*2)
relax(MDOT,linrlx,0.1*2)
relax( u1 ,falsdt,1.0e-4*10)
relax( v1 ,falsdt,1.0e-4*10)
relax( w1 ,falsdt,1.0e-4*10)
relax( h1 ,falsdt,0.01*10)
relax( ke ,falsdt,0.01*10)
relax( ep ,falsdt,0.01*10)
relax( u2 ,falsdt,1.0e-4*10)
relax( v2 ,falsdt,1.0e-4*10)
relax( w2 ,falsdt,1.0e-4*10)
relax( h2 ,falsdt,0.01*10)
relax(MIXF,falsdt,0.01*10) 
relax( FUE,linrlx,0.1*2)
relax( GAS,linrlx,0.1*2)
relax(SHAD,linrlx,0.1*2)  
IF (WWOOD) THEN
+ relax(FWD ,falsdt,0.01*10)
+ relax(WOOD,falsdt,0.01*10)
+ relax(CHAR,falsdt,0.01*10)
ENDIF
 ************************************************************
  Group 18. Limits
VARMIN( H1 ) = 1.0E4 ; VARMAX( H1 ) = 1.0E8
VARMIN( H2 ) = 1.0E4 ; VARMAX( H2 ) = 1.0E8
VARMIN(TMP1) = 200.0 ; VARMAX(TMP1) = 4.0e3
VARMIN(TMP2) = 500.0 ; VARMAX(TMP2) = 4.0e3

VARMIN(RHO1) = 0.05  ; VARMAX(RHO1) = 100.0 
VARMIN(MIXF) = 0.0   ; VARMAX(MIXF) = CMDTC
VARMIN(SHAD) = 1.E-10; VARMAX(SHAD) = 1.0
VARMIN(FUE)  = 1.E-10; VARMAX(FUE)  = 1.0 
VARMIN(GAS)  = 1.E-10; VARMAX(GAS)  = 1.0 
 ************************************************************
  Group 19. EARTH Calls To GROUND Station
USEGRD=T; USEGRX=T
NAMGRD=NONE
IF (WWOOD) THEN
    
    **************************************************
    ****  Notes on combustion laws and constants  ****
    **************************************************
 
    Each of (i) pyrolysis and (ii) char generation can be modelled either
    as a power law or as an Arrhenius law.  However, the same law must be
    used for both.

    Note also that currently only a single set of constants CHSOA ... CHSOE
    are provided to cover both pyrolysis and char generation. 

    In the rate formulae below, "wood" means wood mass fraction, ditto for
    "char", "co2", etc.

    *********************
    ****  Pyrolysis  ****
    *********************
  
PATCH(CHSOW-,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
   
    Pyrolysis (wood => volatiles) can be set either with a power law or an
    Arrhenius relation (latter currently active).  
   
    Power law.  Rate = -chsoa.chsoe.wood.T**chsod
    =========
    COVAL(CHSOW-,WOOD,GRND5,0.0)
    Constants for power law.
    CHSOA=1.
    CHSOD=5.
    CHSOE=10.0*(1.E-3)**CHSOD
    CHSOB=0;CHSOC=0.0
 
    Arrhenius.  Rate = -chsoa.wood.e**(chsod/T)
    =========
COVAL(CHSOW-,WOOD,GRND6,0.0)
 
    Arrhenius rate constants.
real(acte,gascon)
chsoa=woodburn
    Activation energy ACTE in kj/mol (hence GASCON div by 1000)
acte=2.50e+2
gascon=8.314
chsod=-acte/gascon
    Multipliers for redundant terms in FN's in GXSOR (do not change these)
chsob=0.; chsoc=0.0; chsoe=0.0
 
    **************************
    ****  Char creation   ****
    **************************

PATCH(CHSOC+,PHASEM,1,NX,1,NY,1,NZ,1,1)
                              
    Char creation (wood => char) can be set either with a power law or an
    Arrhenius relation (latter currently active). (This is similar to the 
    pyrolysis formation above, and the same choice must be used for both.)
   
    Power law.  Rate = -chsoa.chsoe.wood.T**chsod
    =========
    COVAL(CHSOC+,CHAR,FIXFLU,GRND1)
   
    For constants see power law section of Pyrolysis (above).
   
    Arrhenius.  Rate = -chsoa.wood.*e**(chsod/T)
    =========
COVAL(CHSOC+,CHAR,FIXFLU,GRND2)

SPEDAT(SET,CHSOC+,REACT4,C,WOOD)
SPEDAT(SET,CHSOC+,CONSTA,R,CHSOA*(1.0-VLINWD))
SPEDAT(SET,CHSOC+,CONSTB,R,0.0)

    For constants see Arrhenius section of Pyrolysis (above).

    *************************
    ****  Char burning   ****
    *************************
 
PATCH(CHSOC-,PHASEM,1,NX,1,NY,1,NZ,1,1)                                         

    Char burning is assumed to obey the following relation.
    rate of burn =  - char . ( C.co2 + D.h2o + E.o2)

    The constants C, D, E are currently assumed to be equal and to have
    the value CHARBURN.
 
    Covals and spedats for CHSO1 patch.
 
COVAL(CHSOC-,CHAR,GRND1,0.0)
  
SPEDAT(SET,CHSOC-,REACT1,C,YCO2)
SPEDAT(SET,CHSOC-,REACT2,C,YH2O)
SPEDAT(SET,CHSOC-,REACT3,C,YO2)
SPEDAT(SET,CHSOC-,CONSTC,R,CHARBURN)
SPEDAT(SET,CHSOC-,CONSTD,R,CHARBURN)
SPEDAT(SET,CHSOC-,CONSTE,R,CHARBURN)
 
    **********************************************
    ****  End of combustion setting section   ****
    **********************************************
ENDIF

   ************************************************************
      GROUP 19. Data communicated by SATELLITE to GROUND
   ************************************************************
      Group 20. Preliminary Printout
ECHO    =    T
   ************************************************************
  Group 21. Print-out of Variables
INIFLD=F
output( p1 ,Y,n,n,Y,Y,Y)
OUTPUT( U1 ,Y,n,n,Y,Y,Y)
OUTPUT( V1 ,Y,n,n,Y,Y,Y)
OUTPUT( W1 ,Y,n,n,Y,Y,Y)
OUTPUT( H1 ,Y,n,n,Y,Y,Y)
OUTPUT( V2 ,Y,n,n,Y,Y,Y)
OUTPUT( H2 ,Y,n,n,Y,Y,Y)
OUTPUT(TMP1,Y,n,n,Y,Y,Y)
OUTPUT(SHAD,N,N,N,N,N,N)
OUTPUT( U2 ,N,N,N,N,N,N)
OUTPUT( W2 ,N,N,N,N,N,N)

 ************************************************************
  Group 22. Monitor Print-Out
IXMON=NX*2/3; IYMON=NY-2; IZMON=NZ/2
XZPR=T
 ************************************************************
  Group 24. Dumps For Restarts
 ************************************************************
     restrt(p1,gas,fue,shad,u1,v1,w1,u2,v2,w2,h1,h2,ke,ep,mixf)
solutn(h1,y,y,y,p,p,p)
solutn(h2,y,y,y,p,p,p)
STOP