c

c     file name           GXRADIAT.HTM         160504
C**** SUBROUTINE GXRADI provides coding for the composite-flux
C     and radiosity models for describing thermal radiation in a
C     participating media. The temperature may be solved either directly
C     (via TEM1) and indirectly (via H1) .
C
C     Subroutine GXRADI is called from the following Groups of
C     subroutine GREX3:-
C
C     * Group  1, Section 1:  in order to define various constants
C       in the radiation models and to allocate any storage.
C
C     * Group  9, Section 7:  in order to define the diffusivity for
C       the diffusion fluxes in the equations for the composite fluxes
C       and radiosity.
C
C     * Group 13, Section 13: in order to introduce the source terms
C       in the energy equation (TEM1 or H1) and the equations for the
C       composite fluxes and the radiosity.
C
C How the subroutine is activated:
C
C     In the Q1 file, the user sets the following PIL commands:
C        RADIAT(MODEL,ABSORB,SCAT,ENERGY)
C     where MODEL  = FLUX for the composite-flux model or
C                  = RADI for the radiosity      model;
C           ABSORB = absorption coefficient in m**-1 (=RADIA)
C           SCAT   = scattering coefficient in m**-1 (=RADIB)
C           ENERGY = H1   for energy eqn solved via enthalpy or
C                  = TEM1 for energy eqn solved via temperature.
C
C Limitations:
C
C    The radiation models are coded for ONEPHS=T only, and a gray
C    medium is presumed. The composite-flux model is not valid for
C    BFC=T, although the radiosity model can be used in such
C    applications.
C
C For further information see PHOENICS Encyclopedia entry:
C                                      'Radiation Heat Transfer'
C
      SUBROUTINE GXRADI(CARTES,TEMP1,DEN1,RHO1)
C----------------------------------------------------------------
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      COMMON/IDATA/NX,NY,NZ,LUPR1,IGFILL(116)
      COMMON/RDATA/TINY,RDFIL1(38),HUNIT,RDFIL2(45)
      COMMON/GENI/NXNY,IGFIL(48),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
     1 IPRPS,IRADX,IRADY,IRADZ,IVFOL
      LOGICAL CARTES,NEZ
      INTEGER TEMP1,DEN1,TEMVAR
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      common/dbe/dbear
      logical dbear
      SAVE CORAD,COTEM,FACRAD,VALTEM,IEMPO,ISRAD,NFLX,SIGMA,TEMVAR
C
      NAMSUB = 'GXRADI'
      IF(IGR.EQ.1) THEN ! Group 1
C.... SIGMA is the Stefan-Boltzmann constant in watts m**-2 K**-4
        SIGMA = 5.6697E-8
        TEMVAR= H1
        IF(SOLVE(ITEM1)) TEMVAR=ITEM1
        CALL SUB2(ISRAD,LBNAME('SRAD'),IEMPO,LBNAME('EMPO'))
        IF(ISRAD.GT.0) THEN  ! radiosity model
          CALL SUB3R(COTEM,4.*RADIA*SIGMA,VALTEM,4.*RADIA,
     1               CORAD,4.*RADIA)
        ELSE                 ! composite-flux (i.e. six-flux) model
          NFLX=0         ! count number of flux-pairs, nflx
          IF(SOLVE(IRADX)) NFLX = NFLX + 1
          IF(SOLVE(IRADY)) NFLX = NFLX + 1
          IF(SOLVE(IRADZ)) NFLX = NFLX + 1
          IF(NFLX.EQ.0) THEN
            WRITE(LUPR1,*)'GXRADI has been called because'
            WRITE(LUPR1,*)'a patch named RADI.. exists, but'
            WRITE(LUPR1,*)'no radiation flux is being solved'
            CALL SET_ERR(227,'Error. See result file.',1)
C            CALL EAROUT(1)
          ENDIF
          CORAD  = RADIA + FLOAT(NFLX-1)*RADIB/FLOAT(NFLX)
          COTEM  = FLOAT(2*NFLX)*RADIA*SIGMA
          IF(CORAD.GT.0.0) THEN
            FACRAD = RADIB/(FLOAT(NFLX)*CORAD)
          ELSE
            FACRAD=0.0
          ENDIF
        ENDIF
      ELSEIF(IGR.EQ.9) THEN  ! Group 9 for exchange coefficient
        IF(INDVAR.EQ.IRADX.OR.INDVAR.EQ.IRADY.OR.INDVAR.EQ.IRADZ
     1       .OR.INDVAR.EQ.ISRAD) THEN
C.... PIL variables RADIA and RADIB are the absorption coefficient and
c         scattering coefficients respectively
          CONST=RADIA+RADIB
          IF(INDVAR.EQ.ISRAD.OR.CARTES.OR.INDVAR.NE.IRADY) THEN
            RCONST=1.0/CONST
            IF(ISRAD.GT.0) RCONST=4.0*RCONST/3.0
            CALL FN1(LAMPR,RCONST)  ! fn1(y,a) sets  y = a
          ELSE
            CALL FN7(LAMPR,RV2D,0.0,1.0,CONST) ! fn7(y,x,a,b,c) sets
                                               ! y = (a + b*x)/(1.+c*x)
          ENDIF
          CALL FN27(LAMPR,DEN1)
        ENDIF
      ELSEIF(IGR.EQ.13) THEN  ! Group 13: Sources for composite-flux variables
                              !           and enthalpy or temperature
        IF(INDVAR.EQ.IRADX.OR.INDVAR.EQ.IRADY.OR.INDVAR.EQ.IRADZ
     1     .OR.INDVAR.EQ.ISRAD)                                  THEN
          IF(ISC.EQ.2) THEN      ! CO
            CALL FN1(CO,CORAD)
            CALL FNZERO_IN_SLD_OR_POR(CO)
          ELSEIF(ISC.EQ.13) THEN ! VAL
            ITEMP=ITWO(TEMP1,TEMVAR,TEMVAR.EQ.H1)
            IF(NEZ(TEMP0)) THEN  ! put absolute temperature in EASP1
              CALL FN2(EASP1,ITEMP,TEMP0,1.0)
              ITEMP=EASP1
            ENDIF
            CALL FN63(VAL,ITEMP,SIGMA) ! fn63(y,x,a) sets y = a*x**4
            IF(IEMPO.GT.0) THEN        ! store emissive power 3D
              CALL FN0(IEMPO,VAL)
            ENDIF
            IF(ISRAD.EQ.0.AND.NFLX.GT.1) THEN
              CALL FN25(VAL,RADIA/(CORAD+TINY)) ! fn25(y,a) sets y = a*y
              IF(INDVAR.EQ.IRADX) THEN
                CALL SUB2(MRJ,IRADY,MRK,IRADZ)
              ELSEIF(INDVAR.EQ.IRADY) THEN
                CALL SUB2(MRJ,IRADX,MRK,IRADZ)
              ELSE
                CALL SUB2(MRJ,IRADX,MRK,IRADY)
              ENDIF
              IF(MRJ.GT.0) CALL FN34(VAL,MRJ,FACRAD) ! fn34(y,x,a) sets
                                                     ! y = y + a * x
              IF(MRK.GT.0) CALL FN34(VAL,MRK,FACRAD)
              CALL L0F1(VAL,I,IADD,'GXRADI')
            ENDIF
            CALL FNZERO_IN_SLD_OR_POR(VAL)
          ENDIF
        ELSEIF(INDVAR.EQ.TEMVAR) THEN ! Thermal energy source from radiation
          IF(TEMVAR.EQ.H1) THEN       ! Enthalpy
            IF(NEZ(TEMP0)) THEN  ! Put absolute temperature in EASP1
              CALL FN2(EASP1,TEMP1,TEMP0,1.0)
              ITEMP=EASP1
            ELSE                      ! temp0 = 0
              ITEMP=TEMP1
            ENDIF
            IF(ISC.EQ.2) THEN         ! CO
              CALL FN22(ISPH1,1.E-10)
              CALL FN28(CO,ISPH1,4.*COTEM*HUNIT)
              CALL FN37(CO,ITEMP,3.0)
            ELSEIF(ISC.EQ.13) THEN    ! VAL
              IF(ISRAD.GT.0) THEN
                CALL FN2(VAL,ISRAD,0.0,VALTEM) ! fn2(y,x,a,b) sets y = a + b*x
              ELSE
                CALL FN1(VAL,0.0)
                IF(IRADX.GT.0) CALL FN34(VAL,IRADX,2.*RADIA)
                IF(IRADY.GT.0) CALL FN34(VAL,IRADY,2.*RADIA)
                IF(IRADZ.GT.0) CALL FN34(VAL,IRADZ,2.*RADIA)
              ENDIF
              CALL L0F4(VAL,H1,ITEMP,CO,I,IH1,ITM,ICO,IADD,'GXRADI')
              DO IIX = IXF,IXL
                I = I + IADD
                DO IIY = IYF,IYL
                  I = I + 1
                  TEM4= F(I+ITM)**4
                  F(I)= F(I+IH1) + HUNIT*(F(I)-COTEM*TEM4)/
     1                            (F(I+ICO) + TINY)
                ENDDO
              ENDDO
              CALL FNZERO_IN_SLD_OR_POR(VAL)
            ENDIF
          ELSE                        ! true temperature
            IF(NEZ(TEMP0)) THEN
              CALL FN2(EASP1,TEMVAR,TEMP0,1.0)
              ITEMP=EASP1
            ELSE
              ITEMP=TEMVAR
            ENDIF
            IF(ISC.EQ.2) THEN                 ! CO
              CALL FN61(CO,ITEMP,COTEM)
              CALL FNZERO_IN_SLD_OR_POR(CO)
            ELSEIF(ISC.EQ.13) THEN            ! VAL
              IF(ISRAD.GT.0) THEN
                CALL FN2(VAL,ISRAD,0.0,VALTEM)
              ELSE
                CALL FN1(VAL,0.0)
                IF(IRADX.GT.0) CALL FN34(VAL,IRADX,2.*RADIA)
                IF(IRADY.GT.0) CALL FN34(VAL,IRADY,2.*RADIA)
                IF(IRADZ.GT.0) CALL FN34(VAL,IRADZ,2.*RADIA)
              ENDIF
              CALL FN62(VAL,ITEMP,1.0/(COTEM+TINY))! FN62(Y,X,A) sets
                                                   !      Y = A * Y/X**3
              IF(NEZ(TEMP0)) CALL FN33(VAL,-TEMP0) ! FN33(Y,A) sets Y = Y + A
            ENDIF
            CALL FNZERO_IN_SLD_OR_POR(VAL)
          ENDIF
        ENDIF
      ENDIF
      NAMSUB = 'gxradi'
      END
c