CASE (MIRROR_BOUNDARY)                         WALL(IW)%ILW(N,IBND) = WALL(IW)%ILW(DLM(N,ABS(IOR)),IBND)                        IL(II,JJ,KK) = WALL(IW)%ILW(N,IBND)                     CASE (INTERPOLATED_BOUNDARY)                         IL(II,JJ,KK) = WALL(IW)%ILW(N,IBND)                     CASE DEFAULT ! solid wall                        WALL(IW)%ILW(N,IBND) = OUTRAD_W(IW) + RPI*(1._EB-E_WALL(IW))* INRAD_W(IW)                  END SELECT               ELSEIF (CYLINDRICAL) THEN                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) CYCLE WALL_LOOP1                  IL(II,JJ,KK) = WALL(IW)%ILW(N,IBND)               ENDIF            ENDDO WALL_LOOP1 ! Determine sweep direction in physical space             ISTART = 1            JSTART = 1            KSTART = 1            IEND   = IBAR            JEND   = JBAR            KEND   = KBAR            ISTEP  = 1            JSTEP  = 1            KSTEP  = 1            IF (DLX(N) < 0._EB) THEN               ISTART = IBAR               IEND   = 1               ISTEP  = -1            ENDIF            IF (DLY(N) < 0._EB) THEN               JSTART = JBAR               JEND   = 1               JSTEP  = -1            ENDIF            IF (DLZ(N) < 0._EB) THEN               KSTART = KBAR               KEND   = 1               KSTEP  = -1            ENDIF             GEOMETRY: IF (CYLINDRICAL) THEN  ! Sweep in axisymmetric geometry               J = 1               CKLOOP: DO K=KSTART,KEND,KSTEP                  CILOOP: DO I=ISTART,IEND,ISTEP                     IC = CELL_INDEX(I,J,K)                     IF (SOLID(IC)) CYCLE CILOOP                     ILXU = IL(I-ISTEP,J,K)                     ILYU = IL(I,J-JSTEP,K)                     ILZU = IL(I,J,K-KSTEP)                     IF (DLX(N)>=0._EB) THEN                        RU  = R(I-1)                        RD  = R(I)                     ELSE                        RU  = R(I)                        RD  = R(I-1)                     ENDIF                     RP  = SQRT(0.5_EB*(RU**2+RD**2))                     VC  = DX(I)  * RP*DPHI0 * DZ(K)                     AXU =          RU       * DZ(K) * ABS(DLX(N))                     AXD =          RD       * DZ(K) * ABS(DLX(N))                     AYU = DX(I)             * DZ(K) * ABS(DLB(N))                     AYD = DX(I)             * DZ(K) * ABS(DLY(N))                     AZ  = DX(I)  * RP*DPHI0         * ABS(DLZ(N))                     IF (MODULO(N,NRP(1))==1) AYD = 0._EB  ! Zero out the terms involving symmetric overhang                     IF (MODULO(N,NRP(1))==0) AYU = 0._EB                     IF (IC/=0) THEN                        IW = WALL_INDEX(IC,-ISTEP)                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                           ILXU = WALL(IW)%ILW(N,IBND)                           AYU = 0.5*AYU                           AZ = 0.5*AZ                        ENDIF                        IW = WALL_INDEX(IC,-JSTEP*2)                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                           ILYU = WALL(IW)%ILW(N,IBND)                           AXU = 0.5*AXU                           AZ = 0.5*AZ                        ENDIF                        IW = WALL_INDEX(IC,-KSTEP*3)                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                           ILZU = WALL(IW)%ILW(N,IBND)                           AXU = 0.5*AXU                           AYU = 0.5*AYU                        ENDIF                     ENDIF                     RAP = 1._EB/(AXD+AYD+AZ+EXTCOE(I,J,K)*VC*RSA(N))                     IL(I,J,K) = MAX(0._EB, RAP * (AXU*ILXU + AYU*ILYU + AZ*ILZU +  &                                 VC*RSA(N)*RFPI*( KFST4(I,J,K)+KFST4W(I,J,K) +RSA_RAT*SCAEFF(I,J,K)*UIIOLD(I,J,K) ) ) )                  ENDDO CILOOP               ENDDO CKLOOP            ELSEIF (TWO_D) THEN GEOMETRY  ! Sweep in 2D cartesian geometry               J = 1               K2LOOP: DO K=KSTART,KEND,KSTEP                  I2LOOP: DO I=ISTART,IEND,ISTEP                     IC = CELL_INDEX(I,J,K)                     IF (SOLID(IC)) CYCLE I2LOOP                     ILXU  = IL(I-ISTEP,J,K)                     ILZU  = IL(I,J,K-KSTEP)                     VC  = DX(I) * DZ(K)                     AX  =         DZ(K) * ABS(DLX(N))                     AZ  = DX(I)         * ABS(DLZ(N))                     IF (IC/=0) THEN                        IW = WALL_INDEX(IC,-ISTEP)                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                           ILXU = WALL(IW)%ILW(N,IBND)                           AZ = 0.5*AZ                        ENDIF                        IW = WALL_INDEX(IC,-KSTEP*3)                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                           ILZU = WALL(IW)%ILW(N,IBND)                           AX = 0.5*AX                        ENDIF                     ENDIF                     RAP = 1._EB/(AX+AZ+EXTCOE(I,J,K)*VC*RSA(N))                     IL(I,J,K) = MAX(0._EB, RAP * (AX*ILXU + AZ*ILZU + &                                 VC*RSA(N)*RFPI*(KFST4(I,J,K)+KFST4W(I,J,K) +  RSA_RAT*SCAEFF(I,J,K)*UIIOLD(I,J,K) ) ) )                   ENDDO I2LOOP               ENDDO K2LOOP            ELSE GEOMETRY  ! Sweep in 3D cartesian geometry               KLOOP: DO K=KSTART,KEND,KSTEP                  JLOOP: DO J=JSTART,JEND,JSTEP                     ILOOP: DO I=ISTART,IEND,ISTEP                        IC = CELL_INDEX(I,J,K)                        IF (SOLID(IC)) CYCLE ILOOP                        ILXU  = IL(I-ISTEP,J,K)                        ILYU  = IL(I,J-JSTEP,K)                        ILZU  = IL(I,J,K-KSTEP)                        VC  = DX(I) * DY(J) * DZ(K)                        AX  =         DY(J) * DZ(K) * ABS(DLX(N))                        AY  = DX(I)         * DZ(K) * ABS(DLY(N))                        AZ  = DX(I) * DY(J)         * ABS(DLZ(N))                        IF (IC/=0) THEN                           IW = WALL_INDEX(IC,-ISTEP)                           IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                              AY = 0.5*AY                              AZ = 0.5*AZ                              ILXU = WALL(IW)%ILW(N,IBND)                           ENDIF                           IW = WALL_INDEX(IC,-JSTEP*2)                           IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                              AX = 0.5*AX                              AZ = 0.5*AZ                              ILYU = WALL(IW)%ILW(N,IBND)                           ENDIF                           IW = WALL_INDEX(IC,-KSTEP*3)                           IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN                              AX = 0.5*AX                              AY = 0.5*AY                              ILZU = WALL(IW)%ILW(N,IBND)                           ENDIF                        ENDIF                        RAP = 1._EB/(AX+AY+AZ+EXTCOE(I,J,K)*VC*RSA(N))                        IL(I,J,K) = MAX(0._EB, RAP * ( AX*ILXU + AY*ILYU + AZ*ILZU + &                                    VC*RSA(N)*RFPI*( KFST4(I,J,K)+KFST4W(I,J,K) + RSA_RAT*SCAEFF(I,J,K)*UIIOLD(I,J,K) ) ) )                     ENDDO ILOOP                  ENDDO JLOOP               ENDDO KLOOP             ENDIF GEOMETRY ! Boundary values: Incoming radiation             WALL_LOOP2: DO IW=1,NWC               IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY)   CYCLE WALL_LOOP2                    IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY)   CYCLE WALL_LOOP2                 IF (BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_LOOP2                 IOR = IJKW(4,IW)               IF (TWO_D .AND. .NOT.CYLINDRICAL  .AND. ABS(IOR)==2) CYCLE WALL_LOOP2  ! 2-D non cylindrical               IF (DLN(IOR,N)>=0._EB) CYCLE WALL_LOOP2     ! outgoing               WAXIDLN = - W_AXI*DLN(IOR,N)               IIG = IJKW(6,IW)               JJG = IJKW(7,IW)               KKG = IJKW(8,IW)               INRAD_W(IW) = INRAD_W(IW) - WAXIDLN * WALL(IW)%ILW(N,IBND) ! update incoming radiation,step 1               WALL(IW)%ILW(N,IBND) = IL(IIG,JJG,KKG)               INRAD_W(IW) = INRAD_W(IW) + WAXIDLN * WALL(IW)%ILW(N,IBND) ! update incoming radiation,step 2            ENDDO WALL_LOOP2 ! Copy the Y-downwind intensities to Y-upwind in cylindrical case             IF (CYLINDRICAL) THEN               J=1               DO K=1,KBAR                  DO I=1,IBAR                     IWUP   = WALL_INDEX(CELL_INDEX(I,J,K),-2)                     IWDOWN = WALL_INDEX(CELL_INDEX(I,J,K), 2)                     WALL(IWUP)%ILW(MAX(1,N-1),IBND) = WALL(IWDOWN)%ILW(N,IBND)                  ENDDO               ENDDO            ENDIF ! Calculate integrated intensity UIID             IF (WIDE_BAND_MODEL) THEN               UIID(:,:,:,IBND) = UIID(:,:,:,IBND) + W_AXI*RSA(N)*IL            ELSE               UIID(:,:,:,ANGLE_INC_COUNTER) = UIID(:,:,:,ANGLE_INC_COUNTER) + W_AXI*RSA(N)*IL            ENDIF ! Interpolate boundary intensities onto other meshes             INTERPOLATION_LOOP: DO NOM=1,NMESHES               IF (NM==NOM) CYCLE INTERPOLATION_LOOP               IF (EVACUATION_ONLY(NOM)) CYCLE INTERPOLATION_LOOP               IF (NIC(NOM,NM)==0) CYCLE INTERPOLATION_LOOP               M2=>OMESH(NOM)               OTHER_WALL_LOOP: DO IW=1,MESHES(NOM)%NEWC                  IF (M2%IJKW(9,IW)/=NM .OR. M2%BOUNDARY_TYPE(IW)/=INTERPOLATED_BOUNDARY) CYCLE OTHER_WALL_LOOP                  IOR = M2%IJKW(4,IW)                  IF (DLN(IOR,N)<=0._EB) CYCLE OTHER_WALL_LOOP                  M2%WALL(IW)%ILW(N,IBND)=IL(M2%IJKW(10,IW),M2%IJKW(11,IW),M2%IJKW(12,IW))               ENDDO OTHER_WALL_LOOP            ENDDO INTERPOLATION_LOOP          ENDDO ANGLE_LOOP      ENDDO UPDATE_LOOP       ! Compute incoming flux on walls       DO IW=1,NWC         IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) CYCLE          IBC = IJKW(5,IW)         QRADIN(IW)  = QRADIN(IW) + E_WALL(IW)*(INRAD_W(IW)+SURFACE(IBC)%EXTERNAL_FLUX/NUMBER_SPECTRAL_BANDS)      ENDDO    ENDIF INTENSITY_UPDATE    ! Save source term for the energy equation (QR = -DIV Q)   IF (WIDE_BAND_MODEL) THEN      QR = QR + KAPPA*UIID(:,:,:,IBND)-KFST4      IF (NLP>0 .AND. N_EVAP_INDICIES>0) QR_W = QR_W + KAPPAW*UIID(:,:,:,IBND) - KFST4W   ENDIFENDDO BAND_LOOPIF (UPDATE_INTENSITY) THEN   ! Sum up the parts of the intensity   UII = SUM(UIID, DIM = 4)ENDIF! Save source term for the energy equation (QR = -DIV Q). Done only in one-band (gray gas) case.IF (.NOT. WIDE_BAND_MODEL) THEN   QR  = KAPPA*UII - KFST4   IF (NLP>0 .AND. N_EVAP_INDICIES>0) QR_W = QR_W + KAPPAW*UII - KFST4WENDIFEND SUBROUTINE RADIATION_FVMEND SUBROUTINE COMPUTE_RADIATIONREAL(EB) FUNCTION ZZ2KAPPA(ZZ,ZZF,NIYY,TTD,BAND,NS) ! Calculate KAPPA as a function of ZZ. REAL(EB) :: ZZ,ZZF,IYYINTEGER  :: TTD,BAND,IY1,IY2,ZINT,NIYY,NSREAL(EB), POINTER, DIMENSION(:,:,:) :: KAPPAKAPPA => SPECIES(NS)%KAPPA IF (ZZ<=ZZF .OR. ZZF==1._EB) THEN   IYY = ZZ*REAL(NIYY/2,EB)/ZZFELSE   IYY = NIYY/2 + REAL(NIYY-NIYY/2,EB)*(ZZ-ZZF)/(1._EB-ZZF) ENDIFZINT   = MAX(0,MIN(NIYY,NINT(IYY)))IF(IYY==ZINT) THEN   ZZ2KAPPA=KAPPA(ZINT,TTD,BAND)ELSE   ZINT = MAX(0,MIN(NIYY,NINT(IYY)))   IY2  = MIN(NIYY,ZINT + 1)   IY1  = ZINT   ZZ2KAPPA=KAPPA(IY1,TTD,BAND)+(IYY-REAL(IY1,EB))* (KAPPA(IY2,TTD,BAND)-KAPPA(IY1,TTD,BAND))ENDIFEND FUNCTION ZZ2KAPPA  REAL(EB) FUNCTION IZ2ZZ(IZ,ZZF,NIZ) ! Function gives the mixture fraction as a function of index IZ REAL(EB) :: ZZFINTEGER  :: IZ, NIZ IF (IZ<=NIZ/2) THEN   IZ2ZZ = REAL(IZ,EB)*ZZF/REAL(NIZ/2,EB)ELSE   IZ2ZZ = ZZF+(1._EB-ZZF)*REAL(IZ-NIZ/2,EB)/REAL(NIZ-NIZ/2,EB)ENDIF END FUNCTION IZ2ZZ REAL(EB) FUNCTION BLACKBODY_FRACTION(L1,L2,TEMP) ! Calculates the fraction of black body radiation between wavelengths L1 and L2 (micron) in Temperature TEMP REAL(EB) :: L1,L2,TEMP,LT1,LT2,BBFLOW,BBFHIGHINTEGER  :: IYYLT1    =   L1 * TEMP/LTSTEPLT2    =   L2 * TEMP/LTSTEPIYY = MIN(NLAMBDAT,MAX(0,NINT(LT1)))BBFLOW = BBFRAC(IYY)IYY = MIN(NLAMBDAT,MAX(0,NINT(LT2)))BBFHIGH = BBFRAC(IYY)BLACKBODY_FRACTION = BBFHIGH - BBFLOWEND FUNCTION BLACKBODY_FRACTIONSUBROUTINE GET_KAPPA(Z1,Z2,Z3,YY_SUM,KAPPA,TYY,IBND)USE PHYSICAL_FUNCTIONS, ONLY : GET_F,GET_F_C! GET_KAPPA returns the radiative absorptionREAL(EB) :: Z_1,Z_2,Z_3,YY_SUM,F,C,ZZ,Z_F,KAPPA,KAPPA_N(N_SPECIES),Z1,Z2,Z3INTEGER :: IBND,TYYIF (YY_SUM >=1._EB) THEN   KAPPA = 0._EB   RETURNELSE   YY_SUM = MAX(0._EB,YY_SUM)   Z_1 = MAX(0._EB,Z1)/(1._EB - YY_SUM)   Z_2 = MAX(0._EB,Z2)/(1._EB - YY_SUM)            Z_3 = MAX(0._EB,Z3)/(1._EB - YY_SUM)           ENDIFselectmethod: IF (.NOT. CO_PRODUCTION) THEN   Z_F = REACTION(1)%Z_F   CALL GET_F(Z_1,Z_3,F,Z_F)   ZZ = MIN(1._EB,Z_1 + Z_3)   KAPPA_N(I_FUEL)   = ZZ2KAPPA(ZZ,REACTION(1)%Z_F,SPECIES(I_FUEL)%NKAP_MASS,TYY,IBND,I_FUEL)   KAPPA_N(I_PROG_F) = ZZ2KAPPA(ZZ,REACTION(2)%Z_F,SPECIES(I_PROG_F)%NKAP_MASS,TYY,IBND,I_PROG_F)   KAPPA = (1._EB - F) * KAPPA_N(I_FUEL) + F * KAPPA_N(I_PROG_F)    ELSE selectmethod   CALL GET_F_C(Z_1,Z_2,Z_3,F,C,Z_F)           ZZ = MIN(1._EB,Z_1 + Z_2 + Z_3)   KAPPA_N(I_FUEL)    = ZZ2KAPPA(ZZ,REACTION(1)%Z_F,SPECIES(I_FUEL)%NKAP_MASS,TYY,IBND,I_FUEL)   KAPPA_N(I_PROG_CO) = ZZ2KAPPA(ZZ,REACTION(2)%Z_F,SPECIES(I_PROG_CO)%NKAP_MASS,TYY,IBND,I_PROG_CO)      KAPPA_N(I_PROG_F) = ZZ2KAPPA(ZZ,REACTION(3)%Z_F,SPECIES(I_PROG_F)%NKAP_MASS,TYY,IBND,I_PROG_F)      KAPPA = (1._EB - F) * ( (1._EB - C) * KAPPA_N(I_FUEL) + C * KAPPA_N(I_PROG_CO)) + F *KAPPA_N(I_PROG_F) ENDIF selectmethodKAPPA = KAPPA/(1._EB - YY_SUM)END SUBROUTINE GET_KAPPA SUBROUTINE GET_REV_radi(MODULE_REV,MODULE_DATE)INTEGER,INTENT(INOUT) :: MODULE_REVCHARACTER(255),INTENT(INOUT) :: MODULE_DATEWRITE(MODULE_DATE,'(A)') radirev(INDEX(radirev,':')+1:LEN_TRIM(radirev)-2)READ (MODULE_DATE,'(I5)') MODULE_REVWRITE(MODULE_DATE,'(A)') radidateEND SUBROUTINE GET_REV_radi END MODULE RAD