c
 File Name .....   gxfire.htm ...  031114
c                  called only when that condition is true.
C.... File name ................ gxfire.htm ..................... 220600
C-----------------------------------------------------------------------
C.... This subroutine introduces the heat and mass sources associated
C     with FIRE objects.
C     The Mass sources are:
C        1) No mass source
C        2) Heat related mass source.
C        3) Fixed total mass source.
C        4) POOL fire mass source. (transient only)
C        5) Piece-wise Linear fire mass source. (transient only)
C
C     The Heat Sources are:
C        1) Mass-related heat source.
C        1) Fixed Temperature.
C        2) Fixed total heat flux.
C        4) Heat source Linear with Temperature.
C        5) Heat source power of time. (transient only)
C        6) Piece-wise Linear fire heat source. (transient only)
C
C     The POOL fire mass source is calculated from:
C         Area(t) = a + b*t^c
C                   where a,b,c are constants & t is time.
C         Mass(t) = Area(t)*(1.0 - exp(-B*Area(t)^0.5))
C                   where B is a constant
C
C     The mass-related heat source is calculated from
C         Q = Mass * (Efficiency*Heat of Combustion + Cp*Tpre)
C             where Tpre is the pre-combustion temperature, and
C             Mass is calculated using one of the above expressions.
C             Mass * (Efficiency*Heat of Combustion) is added by the FIRnA patch
C             and Mass * (Cp*Tpre) is added by FIRn
C
C     The Linear with Temperature heat source is calculated from
C         Q = a + b*(TEM1+TEMP0) when Tmin <= TEM1 <= Tmax
C         When TEM1 < Tmin, Tmin is used in the expression
C         When TEM1 > Tmax, Tmax is used in the expression
C     Note that energy sources are multiplied by (1-radiative factor) through
C     adjustment of the constants in Satellite. When mass sources are derived
C     from a heat source, the radiative factor is ignored.
C     The efficiency, EFFI, below, is set in Satellite to be (1-Rf)/(1+Rox)
C-----------------------------------------------------------------------
      SUBROUTINE GXFIRE
      INCLUDE 'farray'
      INCLUDE 'patcmn'
      INCLUDE 'satear'
      INCLUDE 'satgrd'
      INCLUDE 'grdloc'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      COMMON /INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
     1        L0TTFL,L0TTLS,L0MAXC,L0MAXV,L0MINV,L0RATE,L0MAXI,L0NETT,
     1        L0TTFLM,L0TTFLX,L0TTFLY,L0TTFLZ
      COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
     1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
      COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
      COMMON /CMOBTP/OBJTYP(20) /CMNOTP/NMOTYP,NMCHAM
      CHARACTER*14  OBJTYP, CTEMP, BUFF
      REAL TIMVAL(11),MASSVAL(11), QTIM(11),QDOT(11)
      LOGICAL NEZ,QLE,QGE,QGT,MASKPT,EQZ,MASSOR
C
      VALMASS=0.0
      IF(IGR.EQ.13) THEN
C--- Group 9
        IF(ISC.EQ.9) THEN
C------------------- SECTION  9 ------------- coefficient = GRND8
C... initialise CO
          CALL FN1(CO,ONLYMS)
        ELSEIF(ISC.EQ.20) THEN
C------------------- SECTION 20 ------------------- value = GRND8
c...............Note: FN28(Y,X,A)   --> Y= A/X
c                     FN2 (Y,X,A,B) --> Y= A+B*X
c                     FN1 (Y,A)     --> Y= A
C... check if patch is object-related
          IF(MASKPT(IPAT)) THEN
            CTEMP='^'//NPATCH
          ELSE
            CTEMP=NPATCH
          ENDIF
C... get time limits
          CALL GETSDR('TSRT',CTEMP,TSTART)
          CALL GETSDR('TDUR',CTEMP,TDUR)
C... get mass source type
          IMAS=0; CALL GETSDI('IMAS',CTEMP,IMAS)
          MASSOR=IMAS.GT.0
          IMAS=IABS(IMAS)
C... get heat source type
          IQHT=0; CALL GETSDI('IQHT',CTEMP,IQHT)
C... get scalar source type
          ISCAL=0; CALL GETSDI('ISCAL',CTEMP,ISCAL)
C... get Mass source parameters
          IF(IMAS.EQ.2.OR.IQHT.EQ.1.OR.ISCAL.EQ.3) THEN
C... Heat-related source
            CALL GETSDR('EFFI',CTEMP,EFFI)
            CALL GETSDR('HCMB',CTEMP,HCMB)
          ENDIF
          IF(IMAS.EQ.3) THEN
C... Fixed mass source
            FMAS=0.0; CALL GETSDR('FMAS',CTEMP,FMAS)
          ELSEIF(IMAS.EQ.4) THEN
C... POOL fire
            CALL GETSDR('BETA',CTEMP,BETA)
            CALL GETSDR('COFA',CTEMP,COFA)
            CALL GETSDR('COFB',CTEMP,COFB)
            CALL GETSDR('COFC',CTEMP,COFC)
            RAREA=0.0
            CALL GETSDR('AREA',CTEMP,RAREA)
          ELSEIF(IMAS.EQ.5) THEN
C... Piecewise Linear fire
            CALL GETSDI('NSEG',CTEMP,NSEG)
            DO I=1,NSEG+1
              IF(I-1.LT.10) WRITE(BUFF,'(''TIM_'',I1)') I-1
              IF(I-1.GE.10) WRITE(BUFF,'(''TIM_'',I2)') I-1
              CALL GETSDR(BUFF,CTEMP,TIMVAL(I))
              IF(I-1.LT.10) WRITE(BUFF,'(''MAS_'',I1)') I-1
              IF(I-1.GE.10) WRITE(BUFF,'(''MAS_'',I2)') I-1
              CALL GETSDR(BUFF,CTEMP,MASSVAL(I))
            ENDDO
          ENDIF
C... get Heat source parameters
          IF(IQHT.EQ.1) THEN
C... Mass-related source
            CALL GETSDR('EFFI',CTEMP,EFFI)
            CALL GETSDR('HCMB',CTEMP,HCMB)
          ELSEIF(IQHT.EQ.2) THEN
C... fixed temperature
            CALL GETSDR('FTEM',CTEMP,FTEM)
          ELSEIF(IQHT.EQ.3) THEN
C... Fixed source
            CALL GETSDR('FFLU',CTEMP,FFLU)
          ELSEIF(IQHT.EQ.4) THEN
C... Linear with Temperature
            CALL GETSDR('COEFA',CTEMP,COEFA)
            CALL GETSDR('COEFB',CTEMP,COEFB)
            CALL GETSDR('TMIN',CTEMP,TMIN)
            CALL GETSDR('TMAX',CTEMP,TMAX)
C... Power of time
          ELSEIF(IQHT.EQ.5) THEN
            CALL GETSDR('COEFA',CTEMP,COEFA)
            CALL GETSDR('COEFB',CTEMP,COEFB)
            CALL GETSDR('T0',CTEMP,TMIN)
            CALL GETSDR('QMAX',CTEMP,TMAX)
          ELSEIF(IQHT.EQ.6) THEN
C... Piecewise Linear fire
            CALL GETSDI('NSEGQ',CTEMP,NSEGQ)
            DO I=1,NSEGQ+1
              IF(I-1.LT.10) WRITE(BUFF,'(''QTIM_'',I1)') I-1
              IF(I-1.GE.10) WRITE(BUFF,'(''QTIM_'',I2)') I-1
              CALL GETSDR(BUFF,CTEMP,QTIM(I))
              IF(I-1.LT.10) WRITE(BUFF,'(''QDOT_'',I1)') I-1
              IF(I-1.GE.10) WRITE(BUFF,'(''QDOT_'',I2)') I-1
              CALL GETSDR(BUFF,CTEMP,QDOT(I))
            ENDDO
          ENDIF
C... get pre-combustion temperature
          CALL GETSDR('PRTM',CTEMP,PRETEM)
C... Set start and end times
          TFIN=TDUR
          TIMABS=TIM
C... For fixed-flux heat or mass source, source is average over timestep
          IF(.NOT.STEADY.AND.(IQHT.GT.4.OR.IMAS.GT.3))
     1                                              TIMABS=TIMABS-0.5*DT
          DEL=0.01*DT
          IF((.NOT.STEADY.AND.(QGE(TIMABS,TSTART-DEL).AND.
     1                         QLE(TIMABS,TFIN+DEL))).OR.STEADY)
     1                                                              THEN
C... check mass source type
            IF(IMAS.EQ.1) THEN
C... No mass source
              VALMASS=0.0
            ELSEIF(IMAS.EQ.2) THEN
C... Heat-related fire
              VALHEAT=FIRE_HEAT(IQHT,EFFIM,HCMB,FTEM,FFLU,VALMASS,
     1                     COEFA,COEFB,TMIN,TMAX,NSEGQ,QTIM,QDOT,TIMABS)
              VALMASS=VALHEAT/(EFFI*HCMB+TINY)
            ELSEIF(IMAS.EQ.3) THEN
C... Fixed-flow-rate fire
              VALMASS=FMAS
            ELSEIF(IMAS.EQ.4) THEN
C... POOL fire - get current area and mass source
              AREA= COFA+ COFB*(TIMABS**COFC)
              AREA=AMAX1(0.,AREA)
              VALMASS= AREA*(1.-EXP(-BETA*(AREA**0.5)))
              IF(NEZ(RAREA)) VALMASS=VALMASS/RAREA
            ELSEIF(IMAS.EQ.5) THEN
C... Piecewise Linear fire - get current mass source
              IF(QLE(TIMABS,TIMVAL(1))) THE  N        ! Before start
                VALMASS=MASSVAL(1)
              ELSEIF(QGE(TIMABS,TIMVAL(NSEG+1))) THEN ! After finish
                VALMASS=MASSVAL(NSEG+1)
              ELSE                                    ! in range
                I=2
                DO WHILE (QGT(TIMABS,TIMVAL(I)))
                  I=I+1
                ENDDO
                VALMASS=MASSVAL(I-1)+(MASSVAL(I)-MASSVAL(I-1))*
     1                      (TIMABS-TIMVAL(I-1))/(TIMVAL(I)-TIMVAL(I-1))
              ENDIF
            ELSEIF(IMAS.EQ.6) THEN
C... From table file using InForm
              LL=LENGZZ(NAMPAT(IPAT))
              IF(NAMPAT(IPAT)(LL:LL).EQ.'A') THEN
                DO IP=1,NUMPAT
                  IF(NAMPAT(IP).EQ.NAMPAT(IPAT)(1:LL-1)) THEN
                    CALL GETSO(IP,9,VALMASS) ! get source from mass_source patch
                    EXIT
                  ENDIF
                ENDDO
              ELSE
                CALL GETSO(IPAT,9,VALMASS) ! get source from mass_source patch
              ENDIF
            ENDIF
            IF(INDVAR.EQ.P1) THEN
              IF(MASSOR.AND.NEZ(VALMASS)) THEN
C... Mass not zero, so set CO=FIXFLU, VAL=mass/FIXFLU
                CALL FN1(CO,FIXFLU)
                VALMASS=VALMASS/FIXFLU
                CALL FN1(VAL,VALMASS)
              ELSE
C... Mass is zero, so set VAL=0
c                CALL FN1(VAL,VALMASS)
                CALL FN1(VAL,0.0)
              ENDIF
            ELSEIF(INDVAR.EQ.ITEM1) THEN
C... Get heat source
              VALHEAT=FIRE_HEAT(IQHT,EFFI,HCMB,FTEM,FFLU,VALMASS,
     1                     COEFA,COEFB,TMIN,TMAX,NSEGQ,QTIM,QDOT,TIMABS)
C... Apply heat source
              IF(IQHT.NE.2) THEN
C... All heat sources except fixed temperature
                IF(MASSOR.AND.NEZ(VALMASS)) THEN
C... Mass source not zero. Set VAL = Tpre. Rest of heat source in separate patch
                  CALL FN1(VAL,PRETEM)
                ELSE
                  IF((VARMAX(INDVAR).LT.1.E10.AND.VARMIN(INDVAR).LT.
     1                                          -1.E10).OR.CONWIZ) THEN
C.... No mass source, but max increment set. Use source linearisation
                    IF(CONWIZ) THEN
                      CPHI=ABS(VALHEAT)/F(L0MAXI+ISL(INDVAR)) ! set cphi=flux/maxi
                    ELSE
                      CPHI=ABS(VALHEAT)/VARMAX(INDVAR)        ! set cphi=flux/varmax
                    ENDIF
                    CALL FN1(CO,CPHI)
                    IF(VALHEAT.GE.0.0) THEN
                      IF(CONWIZ) THEN
                        VARM=F(L0MAXI+ISL(INDVAR))
                      ELSE
                        VARM=VARMAX(INDVAR)
                      ENDIF
                    ELSE
                      IF(CONWIZ) THEN
                        VARM=-F(L0MAXI+ISL(INDVAR))
                      ELSE
                        VARM=-VARMAX(INDVAR)
                      ENDIF
                    ENDIF
                    CALL FN2(VAL,INDVAR,VARM,1.0) ! & val = var+varmax
                  ELSE
C... No mass source, no max increment, set CO=FIXFLU, VAL=Q/FIXFLU
                    CALL FN1(CO,FIXFLU)
                    VALHEAT=VALHEAT/FIXFLU
                    CALL FN1(VAL,VALHEAT)
                  ENDIF
                ENDIF
              ELSE
C... Fixed temperature
                CALL FN1(VAL,VALHEAT)
                CALL FN1(CO,FIXVAL)
              ENDIF
            ELSE
C... Scalar source
              IF(ISCAL.GT.1) THEN
                LL=LENGZZ(NAME(INDVAR))
                CALL GETSDR(NAME(INDVAR)(1:LL)//'_S',CTEMP,SCAL)
                IF(ISCAL.EQ.4) THEN
C... Scalar source - fixed value
                  CALL FN1(CO,FIXVAL)
                ELSEIF(ISCAL.EQ.3) THEN
C... Scalar source - heat-related
                  VALHEAT=FIRE_HEAT(IQHT,EFFI,HCMB,FTEM,FFLU,VALMASS,
     1                     COEFA,COEFB,TMIN,TMAX,NSEGQ,QTIM,QDOT,TIMABS)
                  VALMASQ=VALHEAT/(EFFI*HCMB+TINY)
                  IF(massor.and.IMAS.NE.1) VALMASQ=VALMASQ-VALMASS
C... Set CO to heat-related mass source
                  CALL FN1(CO,VALMASQ)
                ELSEIF(ISCAL.EQ.2) THEN
C... Scalar source - mass-related. Do nothing to CO as already 0.0
                ENDIF
C... Set VAL to scalar value
                CALL FN1(VAL,SCAL)
              ELSE
                CALL FN1(VAL,0.)
              ENDIF
            ENDIF
          ELSE
C... Outside active time range, set VAL=0
            CALL FN1(VAL,0.)
          ENDIF
        ENDIF
      ENDIF
      END
C-----------------------------------------------------------------------
      REAL FUNCTION FIRE_HEAT(IQHT,EFFI,HCMB,FTEM,FFLU,VALMASS,
     1                     COEFA,COEFB,TMIN,TMAX,NSEGQ,QTIM,QDOT,TIMABS)
      INCLUDE 'farray'
      INCLUDE 'patcmn'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE 'satgrd'
      INCLUDE 'grdloc'
      COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
     1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
      REAL QTIM(11),QDOT(11)
      LOGICAL QLE,QGE,QGT,maskpt,lvrcel
C
      VALHEAT=0.0
C... check heat source type
      IF(IQHT.EQ.1) THEN
C... Heat source function of mass source. Set VAL= Ef*Hcmb
        VALHEAT=EFFI*HCMB*VALMASS
      ELSEIF(IQHT.EQ.2) THEN
C... Fixed temperature
        VALHEAT=FTEM
      ELSEIF(IQHT.EQ.3) THEN
C... Fixed heat source
        VALHEAT=FFLU
      ELSEIF(IQHT.EQ.4) THEN
C... Heat source linear in Temperature
        L0VAL=L0F(VAL)
        L0TEM=L0F(ITEM1)
        L0CP1=L0F(ISPH1)
        L0DEN=L0F(DEN1)
        L0VOL=L0F(VOL)
        CALL SUB2R(SUM1,0.0,SUM2,0.0)
        IPV=0
        DO IX=IXF,IXL
          DO IY=IYF,IYL
            I=IY+(IX-1)*NY
            IF(MASKPT(IPAT)) THEN
              IPV=IPV+1
              IF(.NOT.LVRCEL(IPV)) CYCLE ! don't include cells not inside facets
            ENDIF
            SUM1=SUM1+F(L0DEN+I)*F(L0CP1+I)*F(L0VOL+I)*(F(L0TEM+I)+
     1                                                            TEMP0)
            SUM2=SUM2+F(L0DEN+I)*F(L0CP1+I)*F(L0VOL+I)
          ENDDO
        ENDDO
        IF(SUM2.GT.0.0) THEN ! don't include cells not inside facets
          TBAR=SUM1/SUM2-TEMP0
          IF(QLE(TBAR,TMIN)) THEN
            TTOT=TMIN+TEMP0
          ELSEIF(QGE(TBAR,TMAX)) THEN
            TTOT=TMAX+TEMP0
          ELSE
            TTOT=TBAR+TEMP0
          ENDIF
          VALHEAT=COEFA+COEFB*TTOT
        ELSE
          VALHEAT=0.0
        ENDIF
      ELSEIF(IQHT.EQ.5) THEN
C... Heat source power of time
        VALHEAT=AMIN1(TMAX, COEFA*(TIMABS-TMIN)**COEFB)
      ELSEIF(IQHT.EQ.6) THEN
C... Piecewise Linear fire - get current heat source
        IF(QLE(TIMABS,QTIM(1))) THEN           ! Before start
          VALHEAT=QDOT(1)
        ELSEIF(QGE(TIMABS,QTIM(NSEGQ+1))) THEN ! After finish
          VALHEAT=QDOT(NSEGQ+1)
        ELSE                                   ! in range
          I=2
          DO WHILE (QGT(TIMABS,QTIM(I)))
            I=I+1
          ENDDO
          VALHEAT=QDOT(I-1)+(QDOT(I)-QDOT(I-1))*
     1                       (TIMABS-QTIM(I-1))/(QTIM(I)-QTIM(I-1))
        ENDIF
      ELSEIF(IQHT.EQ.7) THEN
C... From table file using InForm
        LL=LENGZZ(NAMPAT(IPAT))
        DO IP=1,NUMPAT
          IF(NAMPAT(IP).EQ.NAMPAT(IPAT)(1:LL)//'A') THEN
            CALL GETSO(IP,ITEM1,VALHEAT) ! get source from heat-source patch
            EXIT
          ENDIF
        ENDDO
      ENDIF
      FIRE_HEAT=VALHEAT
      END
c