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