c

C.... file name GXCHEMSO.F              02.01.00
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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/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
        FIRST=.FALSE.
        LBSOR = LBNAME('CHSO')
        LBWOOD= LBNAME('WOOD')
        LBREA3= C1
C                 to be used for SCRS
        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
          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:  eg  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 1 I=1,NXNY
    1     F(L0CO+I) = C*F(L0REA1+I) + D*F(L0REA2+I) + E*F(L0REA3+I)
        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 11 I=1,NXNY
   11     F(L0VAL+I) = A * F(L0REA4+I) * F(L0TEMP +I) ** B
        ELSEIF(ISC.EQ.14) THEN
C     ISC=14, ie val = GRND2: eg  arrhenius-law production of char,
c                                                              at rate
          L0REA4=L0F(LBREA4)
          L0VAL=L0F(VAL)
          L0TEMP=L0F(TEMP)
          CALL FN22(TEMP,VARMIN(TEMP))
          DO 12 I=1,NXNY
c lj 24.02.97 protection against solid regions
          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
   12     CONTINUE
        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) CALL SUB2(MFUEL,LBNAME('FUEL'),
     1                                  MIXF,LBNAME('MIXF'))
        CALL SUB2(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
            CALL SUB2(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
c.... conditionally store reaction rate in 3D CHSO store
      IF(LBSOR.NE.0) THEN
        IF(ISC.GT.11) CALL FN48(LBSOR,CO,VAL,CO,INDVAR,1.0,-1.0)
      ENDIF
      NAMSUB = 'gxchso'
      END