c
C.... file name GXCHEMSO.F 061110
C**** SUBROUTINE GXCHSO is called from group 13 of GREX3, and is
C entered only when the PATCH name begins with the characters 'CHSO'.
C
C.... The argument TEMP is the temperature of the first-phase,
C TMP1, in the SATELLITE.
C
C.... The library cases 109(CO=GRND7), 492(CO=GRND9 and VAL=GRND9)
C make use of it.
C
SUBROUTINE GXCHSO(TEMP)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON/GENI/NXNY,IFIL5(54),IPRPS,IFIL1(4)
COMMON/RDA6/VARMIN(150) /RDA7/ VARMAX(150)/RDA10/CINT(150)
COMMON/IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
COMMON/IDATA/NX,NY,IFL1(23),FSWEEP,IFL2(64),LEN1,IFIL3(27),VIST,
1 IDFIL3
COMMON/RDATA/TINY,RD1(84)
COMMON/LRNTM2/LBWDIS,JFMU,JFONE,JFTWO,JREYN,JREYT,LBEPKE
INTEGER VAL,CO,WALDIS,PATGEO,TEMP,FSWEEP,VIST
COMMON/NAMFN/NAMFUN,NAMSUB
LOGICAL SLD,WOOD,FIRST
CHARACTER*6 NAMFUN,NAMSUB
CHARACTER*4 RENAM1,RENAM2,RENAM3,RENAM4
SAVE MFUEL,MIXF,LBSOR,LBREA1,LBREA2,LBREA3,LBREA4,A,B,C,D,E,
1 RENAM1,RENAM2,RENAM3,RENAM4,WOOD
DATA FIRST/.TRUE./
C
C CO=GRNDn with n= 1 2 3 4 5 6 7 8 9 10
C.... corresponds to ISC=1 2 3 4 5 6 7 8 9 10 11
C
C VAL=GRNDn with n= 1 2 3 4 5 6 7 8 9 10
C.... corresponds to ISC=12 13 14 15 16 17 18 19 20 21 22
C
C INDVAR=unburned fuel=LBNAME('FU')
C
C functions used in this subroutine are:
c.... fn8(y,x,a,b,c,d) y=a * (x + b)**c + d
c fn10(y,x1,x2,a,b1,b2) y=a + b1*x1 + b2*x2
c fn22(y,a) y=amax1(y,a)
c fn25(y,a) y=a * y
c fn32(y,x,a,b) y=y * exp(a/(b + x))
c fn37(y,x,a) y=y * x**a
c fn48(y,x1,x2,x3,x4,a,b) y=a * x1*x2 + b * x3*x4
c
c
NAMSUB='GXCHSO'
c ISC=6; CO=GRND5 or TEMPOWERA : rate prp. to power of temperature
c coefficient=CHSOE * mox * temperature ** CHSOD
c with mox=CHSOA + CHSOB * C1 + CHSOC * FU
c CHSOA=stoichiometric mixture fraction
c CHSOB=-1.0/CHSOA
c CHSOC=1.0/CHSOA - 1.0
c
IF(FIRST) THEN
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use start
FIRST=.FALSE. ; LBSOR=LBNAME('CHSO') ; LBWOOD=LBNAME('WOOD')
LBREA3=C1 ; LBMIXF=LBNAME('MIXF')
IF(STORE(LBMIXF)) THEN
LBREA3=LBMIXF
ELSE
LBREA3=C1
ENDIF
WOOD=LBWOOD.NE.0
IF(WOOD) THEN
c.... read data inserted by SPEDAT commands in Q1, in first sweep
CALL GETSDC('CHSOC-','REACT1',RENAM1)
CALL GETSDC('CHSOC-','REACT2',RENAM2)
CALL GETSDC('CHSOC-','REACT3',RENAM3)
CALL GETSDR('CHSOC-','CONSTC',C)
CALL GETSDR('CHSOC-','CONSTD',D)
CALL GETSDR('CHSOC-','CONSTE',E)
c
CALL GETSDC('CHSOC+','REACT4',RENAM4)
CALL GETSDR('CHSOC+','CONSTA',A)
CALL GETSDR('CHSOC+','CONSTB',B)
c.... set required location-in-common-block indices
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use start
LBREA1=LBNAME(RENAM1) ; LBREA2=LBNAME(RENAM2)
LBREA3=LBNAME(RENAM3) ; LBREA4=LBNAME(RENAM4)
write(14,*)'Data transmitted by SPEDAT to GXCHSO'
WRITE(14,*)'RENAM1= ',RENAM1,' RENAM2=',RENAM2,
1 ' RENAM3=',RENAM3,' RENAM4=',RENAM4
CALL WRIT2R('A ',A,'B ',B)
CALL WRIT3R('C ',C,'D ',D,'E ',E)
ENDIF
ENDIF
c
c......................... char generation and burning ................
IF(ISC.EQ.2.OR.ISC.EQ.13.OR.ISC.EQ.14) THEN
IF(ISC.EQ.2) THEN
C ISC=2, ie co=grnd1: e.g. char burning, at rate
c - char * ( c * co2 + d * h2o + e * o2)
L0CO=L0F(CO) ; L0REA1=L0F(LBREA1) ; L0REA2=L0F(LBREA2)
L0REA3=L0F(LBREA3)
DO I=1,NXNY
F(L0CO+I)=C*F(L0REA1+I) + D*F(L0REA2+I) + E*F(L0REA3+I)
ENDDO
ELSEIF(ISC.EQ.13) THEN
C ISC=13, ie val=grnd1: eg power-law production of char, at rate
c a * wood_mass_fraction * t ** b
L0REA4=L0F(LBREA4) ; L0VAL=L0F(VAL) ; L0TEMP=L0F(TEMP)
DO I=1,NXNY
F(L0VAL+I)=A * F(L0REA4+I) * F(L0TEMP +I) ** B
ENDDO
ELSEIF(ISC.EQ.14) THEN
C ISC=14, i.e. val=GRND2: e.g. arrhenius-law production of char,
c at rate
L0REA4=L0F(LBREA4) ; L0VAL=L0F(VAL) ; L0TEMP=L0F(TEMP)
CALL FN22(TEMP,VARMIN(TEMP))
DO I=1,NXNY
IF(SLD(I)) THEN
F(L0VAL+I)=0.0
ELSE
F(L0VAL+I)=A * F(L0REA4+I) *
1 EXP(-CHSOD/(F(L0TEMP +I) + CHSOE + 1.E-20))
ENDIF
ENDDO
ENDIF
ELSEIF(ISC.EQ.3) THEN
ELSEIF(ISC.EQ.4) THEN
ELSEIF(ISC.EQ.5) THEN
ELSEIF(ISC.EQ.6) THEN
C.... power-law-of temperature depletion rate
CALL FN10(CO,LBREA3,INDVAR,CHSOA,CHSOB,CHSOC)
CALL FN22(CO,0.0)
CALL FN37(CO,TEMP,CHSOD)
CALL FN25(CO,CHSOE)
C
c ISC=7; CO=GRND6 or ARRHEN : rate prp. to Arrhenius expression
C coefficient=mox * mfuel* exp(CHSOD/(CHSOE+temp.))
C with mox, CHSOA, CHSOB, CHSOC as above
C
ELSEIF(ISC.EQ.7) THEN
C Arrhenius depletion rate
CALL FN10(CO,LBREA3,INDVAR,CHSOA,CHSOB,CHSOC)
CALL FN22(CO,0.0)
CALL FN22(TEMP,VARMIN(TEMP))
CALL FN32(CO,TEMP,CHSOD,CHSOE)
C
c ISC=8; CO=GRND7 or POLYNOM : rate prp. to power of reactedness
C coefficient=CHSOA * concentration ** CHSOB
C
ELSEIF(ISC.EQ.8) THEN
CALL FN8(CO,INDVAR,CHSOA,0.0,CHSOB,0.0)
C
C ISC=10 or 21; CO and VAL=GRND9 or EBU (Eddy-break-up model)
C coefficient=CHSOB * (EP/KE or EPKE) for IEBU=CHSOE=0
C =CHSOB * GEN1 ** 0.5 for IEBU=CHSOE=1
C =CHSOB * ABS(VEL) / LEN1 for IEBU=CHSOE=2
C =CHSOB * VIST / LEN1 ** 2 for IEBU=CHSOE=3
C except coefficient=0 for solid, or for
c CHSOD > 0.0 and (fu-fubrnt)/(mixture fraction - fubrnt) < chsod
C or for
C CHSOD < 0.0 coefficient is multiplied by
c ((fu-fubrnt)/(mixture fraction - fubrnt))**(-chsod)
c
ELSEIF(ISC.EQ.10 .OR. ISC.EQ.21) THEN
IF(ISWEEP.LE.FSWEEP) THEN
MFUEL=LBNAME('FUEL') ; MIXF=LBNAME('MIXF')
ENDIF
L0MIXF=L0F(MIXF) ; L0MFU=L0F(MFUEL)
IF(ISC.EQ.10) THEN
L0CO=L0F(CO)
ELSEIF(ISC.EQ.21) THEN
L0VAL=L0F(VAL)
ENDIF
IEBU=CHSOE + 0.1
IF(IEBU.EQ.0) THEN
IF(LBEPKE.NE.0) THEN
L0EPKE=L0F(LBEPKE)
ELSE
L0KE=L0F(KE) ; L0EP=L0F(EP)
ENDIF
ELSEIF(IEBU.EQ.1) THEN
L0GEN1=L0F(GEN1)
ELSEIF(IEBU.EQ.2) THEN
L0LEN1=L0F(LEN1)
CALL FNVLSQ(CO,1)
ELSEIF(IEBU.EQ.3) THEN
L0LEN1=L0F(LEN1) ; L0VIST=L0F(VIST)
ENDIF
REC1=1.0/(1.0 - CHSOA)
REC2=CHSOA/(1.0-CHSOA)
DO 20 IX=IXF,IXL
IPLUS=(IX-1)*NY
DO 30 IY=IYF,IYL
I=IY + IPLUS ; FU=F(L0MFU+I)
FUBRNT=AMIN1(FU , AMAX1(0.0, (F(L0MIXF+I)-CHSOA)*REC1) )
IF(ISC.EQ.10) THEN
C.... Set coefficient
IF(SLD(I)) THEN
F(L0CO+I)=0.0
GO TO 32
ELSEIF(CHSOD.GT.0.0) THEN
IF( (FU-FUBRNT) .GT. CHSOD*(F(L0MIXF+I)-FUBRNT) ) THEN
F(L0CO+I)=0.0
GO TO 32
ENDIF
c.... conditionally multiply CHSOD by function of FU, FUBRNT & CHSOD
ELSEIF(CHSOD.LT.0.0) THEN
COEF=CHSOB*((FU-FUBRNT)/(F(L0MIXF+I)-FUBRNT))**(-CHSOD)
ELSE
COEF=CHSOB
ENDIF
C.... multiply by IEBU-dependent rate-controlling factor
IF(IEBU.EQ.0) THEN
IF(LBEPKE.EQ.0) THEN
F(L0CO+I)=COEF * F(L0EP+I)/(F(L0KE+I)+TINY)
ELSE
F(L0CO+I)=COEF * F(L0EPKE+I)
ENDIF
ELSEIF(IEBU.EQ.1)THEN
F(L0CO+I)=COEF * SQRT(F(L0GEN1+I))
ELSEIF(IEBU.EQ.2) THEN
VELSQ=F(L0CO+I)
F(L0CO+I)=COEF * SQRT(VELSQ)/(F(L0LEN1+I) + TINY)
ELSEIF(IEBU.EQ.3) THEN
F(L0CO+I)=COEF * F(L0VIST+I)/(F(LEN1+I)**2 + TINY)
ENDIF
32 CONTINUE
ELSEIF(ISC.EQ.21) THEN
C.... Set value
F(L0VAL+I)=FUBRNT
ENDIF
30 CONTINUE
20 CONTINUE
ENDIF
IF(LBSOR.NE.0) THEN ! conditionally store reaction rate in 3D CHSO store
IF(ISC.GT.11) CALL FN48(LBSOR,CO,VAL,CO,INDVAR,1.0,-1.0)
ENDIF
NAMSUB='gxchso'
END