MW_RADCAL = 32            ENDIF         END IF         BAND_LOOP_FR1: DO IBND = 1,NSB            IF (NSB>1) THEN               OMMIN = REAL(NINT(1.E4_EB/WL_HIGH(IBND)),EB)               OMMAX = REAL(NINT(1.E4_EB/WL_LOW(IBND)),EB)            ELSE               OMMIN = 50._EB               OMMAX = 10000._EB            ENDIF            CALL INIT_RADCAL            T_LOOP_MF2: DO K = 0,SS%NKAP_TEMP               RCT(1) = 300._EB + K*(2400._EB-300._EB)/SS%NKAP_TEMP               Z_LOOP_MF2: DO J=0,SS%NKAP_MASS                  YY = REAL(J)*SS%MAXMASS/REAL(SS%NKAP_MASS)                  MTOT = YY/MW_RADCAL+(1._EB-YY)/MW_N2                  SPECIE = 0._EB                  SPECIE(I_RADCAL) = YY                  P(:,1) = 0._EB                  SVF = 0._EB                  P(I_RADCAL,1) = YY/MW_RADCAL/MTOT                  P(6,1) = (1-YY)/MW_N2/MTOT                  CALL RADCAL(AMEAN,AP0)                  IF (NSB==1 .AND. PATH_LENGTH > 0._EB) THEN                     SS%KAPPA(J,K,IBND) = MIN(AMEAN,AP0)                  ELSE                     IF (NSB==1) THEN                        BBF = 1._EB                     ELSE                        BBF = BLACKBODY_FRACTION(WL_LOW(1),WL_HIGH(1),RCT(1))                     ENDIF                     SS%KAPPA(J,K,IBND) = AP0/BBF                  ENDIF               ENDDO Z_LOOP_MF2            ENDDO T_LOOP_MF2         ENDDO BAND_LOOP_FR1      CASE (AEROSOL_SPECIES) GAS_TYPE         SS%NKAP_TEMP = 21         SS%NKAP_MASS = 200           SS%MAXMASS=0.2_EB            ALLOCATE(SS%KAPPA(0:SS%NKAP_MASS,0:SS%NKAP_TEMP,1:NSB),STAT=IZERO)         CALL ChkMemErr('RADI','KAPPA',IZERO)         BAND_LOOP_FR2: DO IBND = 1,NSB            IF (NSB>1) THEN               OMMIN = REAL(NINT(1.E4_EB/WL_HIGH(IBND)),EB)               OMMAX = REAL(NINT(1.E4_EB/WL_LOW(IBND)),EB)            ELSE               OMMIN = 50._EB               OMMAX = 10000._EB            ENDIF            CALL INIT_RADCAL            I_RADCAL=5            T_LOOP_MF3: DO K = 0,SS%NKAP_TEMP               RCT(1) = 300._EB + K*(2400._EB-300._EB)/SS%NKAP_TEMP               Z_LOOP_MF3: DO J=0,SS%NKAP_MASS                  YY = REAL(J)*SS%MAXMASS/REAL(SS%NKAP_MASS)                  RCRHO = 29._EB * P_INF/(R0*RCT(1)) !Good enough                  SPECIE = 0._EB                                 P(:,1) = 0._EB                                 P(6,1) = 1._EB                  SPECIE(I_RADCAL) = YY*RCRHO/RHO_SOOT                                    SVF(1) = YY*RCRHO/RHO_SOOT                                                      CALL RADCAL(AMEAN,AP0)                  IF (NSB==1 .AND. PATH_LENGTH > 0._EB) THEN                     SS%KAPPA(J,K,IBND) = MIN(AMEAN,AP0)                  ELSE                     IF (NSB==1) THEN                        BBF = 1._EB                     ELSE                        BBF = BLACKBODY_FRACTION(WL_LOW(1),WL_HIGH(1),RCT(1))                     ENDIF                     SS%KAPPA(J,K,IBND) = AP0/BBF                  ENDIF               ENDDO Z_LOOP_MF3            ENDDO T_LOOP_MF3         ENDDO BAND_LOOP_FR2   END SELECT GAS_TYPEENDDO SPECIES_LOOP CALL RCDEALLOC  ! Deallocate RadCal arrays  !-----------------------------------------------------!!     Tables for droplet absorption coefficients! !-----------------------------------------------------DROPLETS: IF (N_EVAP_INDICIES>0) THEN   GET_PC_RADI: DO J=1,N_PART      PC => PARTICLE_CLASS(J)      IF ((.NOT. PC%FUEL) .AND. (.NOT. PC%WATER)) CYCLE GET_PC_RADI      IF (PC%FUEL) THEN  ! Tuntomo, Tien and Park, Comb. Sci. and Tech., 84, 133-140, 1992         PC%WQABS(:,1) = (/0,10,16,28,52,98,191,386,792,1629,3272,6163, &            10389,15588,20807,23011,22123,22342,22200,22241,21856, &            22795,23633,24427,25285,26207,27006,27728,28364,28866, &            29260/)!,29552,29748,30000/)         PC%KWR(0) = 0._EB         DO I=1,NDG            PC%KWR(I)=EXP(I/2.5_EB-4._EB)         ENDDO         PC%WQABS=PC%WQABS/10000._EB         PC%KWR=PC%KWR/1000000._EB      ENDIF      IF (PC%WATER) CALL MEAN_CROSS_SECTIONS(J)   ENDDO GET_PC_RADIENDIF DROPLETS ! A few miscellaneous constantsFOUR_SIGMA = 4._EB*SIGMARPI_SIGMA  = RPI*SIGMA ! In axially symmetric case, each angle represents two symmetric angles. So weight the intensities by two.W_AXI = 1._EBIF (CYLINDRICAL) W_AXI = 2._EB END SUBROUTINE INIT_RADIATION SUBROUTINE COMPUTE_RADIATION(NM)! Call radiation routine or simply specify the radiative lossUSE MESH_POINTERSUSE COMP_FUNCTIONS, ONLY : SECOND  REAL(EB) :: TNOWINTEGER, INTENT(IN) :: NM IF (EVACUATION_ONLY(NM)) RETURNTNOW=SECOND() CALL POINT_TO_MESH(NM)IF (RADIATION) THEN   CALL RADIATION_FVM(NM)ELSE   IF (N_REACTIONS>0) QR = -RADIATIVE_FRACTION*QENDIFTUSED(9,NM)=TUSED(9,NM)+SECOND()-TNOWCONTAINS  SUBROUTINE RADIATION_FVM(NM)USE MIEVUSE MATH_FUNCTIONS, ONLY : INTERPOLATE1DREAL(EB) :: ZZ, RAP, AX, AXU, AXD, AY, AYU, AYD, AZ, VC, RU, RD, RP, &            ILXU, ILYU, ILZU, QVAL, BBF, BBFA, NCSDROP, RSA_RAT, WAXIDLN , KAPPA_1, Z_2INTEGER  :: N,IIG,JJG,KKG,I,J,K,IW,II,JJ,KK,IOR,IC,IWUP,IWDOWN, &            ISTART, IEND, ISTEP, JSTART, JEND, JSTEP, &            KSTART, KEND, KSTEP, NSTART, NEND, NSTEP, &            I_UIID, N_UPDATES, IBND, TYY, NOM, NS, IBC,EVAP_INDEXLOGICAL :: UPDATE_INTENSITYREAL(EB), POINTER, DIMENSION(:,:,:) :: KFST4, IL, UIIOLD, KAPPAW, KFST4W, EXTCOE, SCAEFFREAL(EB), POINTER, DIMENSION(:)     :: OUTRAD_W, INRAD_WINTEGER, INTENT(IN) :: NMTYPE (OMESH_TYPE), POINTER :: M2TYPE(SURFACE_TYPE), POINTER :: SFTYPE(PARTICLE_CLASS_TYPE), POINTER :: PCKFST4    => WORK1IL       => WORK2UIIOLD   => WORK3EXTCOE   => WORK4KAPPAW   => WORK5SCAEFF   => WORK6KFST4W   => WORK7OUTRAD_W => WALL_WORK1INRAD_W  => WALL_WORK2 ! Ratio of solid angle, used in scattering RSA_RAT = 1._EB/(1._EB-1._EB/NRA) ! Check if it time to update radiation intensity field RAD_CALL_COUNTER  = RAD_CALL_COUNTER+1IF ( MOD(RAD_CALL_COUNTER,TIME_STEP_INCREMENT)==1 .OR. TIME_STEP_INCREMENT==1) THEN   UPDATE_INTENSITY = .TRUE.ELSE   UPDATE_INTENSITY = .FALSE.ENDIFIF (WIDE_BAND_MODEL) QR = 0._EBIF (UPDATE_INTENSITY) QRADIN = 0._EB ! Loop over spectral bands BAND_LOOP: DO IBND = 1,NUMBER_SPECTRAL_BANDS    KAPPAW = 0._EB   KFST4W = 0._EB   SCAEFF = 0._EB ! Calculate fraction on ambient black body radiation    IF (NUMBER_SPECTRAL_BANDS==1) THEN      BBFA = 1._EB   ELSE      BBFA = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMPA)   ENDIF      ! Generate water absorption and scattering coefficients    IF_DROPLETS_INCLUDED: IF (NLP>0 .AND. N_EVAP_INDICIES>0) THEN      IF (NUMBER_SPECTRAL_BANDS==1) THEN         BBF = 1._EB      ELSE         BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),RADTMP)      ENDIF            DO K=1,KBAR         DO J=1,JBAR            ZLOOPM: DO I=1,IBAR               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ZLOOPM                  PC_LOOP: DO N = 1,N_PART                     PC => PARTICLE_CLASS(N)                     EVAP_INDEX = PC%EVAP_INDEX                     IF (EVAP_INDEX==0) CYCLE PC_LOOP                     IF (AVG_DROP_DEN(I,J,K,EVAP_INDEX)==0._EB) CYCLE PC_LOOP                     NCSDROP = THFO*AVG_DROP_DEN(I,J,K,EVAP_INDEX)/ (PC%DENSITY*AVG_DROP_RAD(I,J,K,EVAP_INDEX))                     ! Absorption and scattering efficiency                     CALL INTERPOLATE1D(PC%KWR,PC%WQABS(:,IBND), 0.95_EB*AVG_DROP_RAD(I,J,K,EVAP_INDEX),QVAL)                      KAPPAW(I,J,K) = KAPPAW(I,J,K) + NCSDROP*QVAL                     IF (PC%WATER) THEN                        CALL INTERPOLATE1D(PC%KWR,PC%WQSCA(:,IBND),0.95_EB*AVG_DROP_RAD(I,J,K,EVAP_INDEX),QVAL)                        SCAEFF(I,J,K) = SCAEFF(I,J,K) + NCSDROP*QVAL                     ENDIF                     KFST4W(I,J,K) = KFST4W(I,J,K)+ BBF*KAPPAW(I,J,K)*FOUR_SIGMA*AVG_DROP_TMP(I,J,K,EVAP_INDEX)**4                  END DO PC_LOOP            ENDDO ZLOOPM         ENDDO      ENDDO      QR_W = 0._EB   ENDIF IF_DROPLETS_INCLUDED ! Compute absorption coefficient KAPPA and source term KAPPA*4*SIGMA*TMP**4    BBF = 1._EB   KAPPA=KAPPA0   DO K=1,KBAR      DO J=1,JBAR         ZLOOP: DO I=1,IBAR            IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ZLOOP            SUM_SPECIES: DO NS=1,N_SPECIES               IF (SPECIES(NS)%MODE == MIXTURE_FRACTION_SPECIES) CYCLE SUM_SPECIES               IF (.NOT. SPECIES(NS)%ABSORBING)                  CYCLE SUM_SPECIES               TYY = NINT(SPECIES(NS)%NKAP_TEMP*(TMP(I,J,K)-300._EB)/2100._EB)               TYY = MAX(0,MIN(SPECIES(NS)%NKAP_TEMP,TYY))               ZZ = MIN(1._EB,MAX(0._EB,YY(I,J,K,NS)))               KAPPA(I,J,K) = KAPPA(I,J,K) + ZZ2KAPPA(ZZ,SPECIES(NS)%MAXMASS,2*SPECIES(NS)%NKAP_MASS,TYY,IBND,NS)            ENDDO SUM_SPECIES            IF (MIXTURE_FRACTION) THEN               TYY = NINT(SPECIES(I_FUEL)%NKAP_TEMP*(TMP(I,J,K)-300._EB)/2100._EB)               TYY = MAX(0,MIN(SPECIES(I_FUEL)%NKAP_TEMP,TYY))               IF (CO_PRODUCTION) THEN                  Z_2 = YY(I,J,K,I_PROG_CO)               ELSE                  Z_2 = 0._EB               ENDIF               CALL GET_KAPPA(YY(I,J,K,I_FUEL),Z_2,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),KAPPA_1,TYY,IBND)                             KAPPA(I,J,K) = KAPPA(I,J,K) + KAPPA_1            ENDIF            IF (WIDE_BAND_MODEL)  BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMP(I,J,K))            KFST4(I,J,K) = BBF*KAPPA(I,J,K)*FOUR_SIGMA*TMP(I,J,K)**4            IF (LES) THEN               IF (BBF*RADIATIVE_FRACTION*Q(I,J,K) > KFST4(I,J,K)) THEN                  KFST4(I,J,K) = BBF*RADIATIVE_FRACTION*Q(I,J,K)                  KAPPA(I,J,K) = 0._EB               ENDIF            ENDIF         ENDDO ZLOOP      ENDDO   ENDDO! Calculate extinction coefficient    EXTCOE = KAPPA + KAPPAW + SCAEFF*RSA_RAT! Update intensity field    INTENSITY_UPDATE: IF (UPDATE_INTENSITY) THEN       IF (WIDE_BAND_MODEL) THEN         UIIOLD = UIID(:,:,:,IBND)      ELSE         UIIOLD = UII      ENDIF      UII = 0._EB! Compute boundary condition intensity emissivity*sigma*Tw**4/pi !     or emissivity*QRADOUT/pi for wall with internal radiation       BBF = 1.0_EB      DO IW = 1,NWC         IF (WIDE_BAND_MODEL) BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMP_F(IW))         IBC = IJKW(5,IW)         SF  => SURFACE(IBC)         IF (.NOT. SF%INTERNAL_RADIATION) THEN            QRADOUT(IW) = E_WALL(IW)*SIGMA*TMP_F(IW)**4         ENDIF         OUTRAD_W(IW) = BBF*RPI*QRADOUT(IW)      ENDDO! Compute boundary condition term incoming radiation integral       DO IW = 1,NWC         IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) CYCLE         IOR = IJKW(4,IW)         INRAD_W(IW) = SUM(-W_AXI*DLN(IOR,:)* WALL(IW)%ILW(:,IBND),1, DLN(IOR,:)<0._EB)      ENDDO ! If updating intensities first time, sweep ALL angles       N_UPDATES = 1      IF (RAD_CALL_COUNTER==1) N_UPDATES = ANGLE_INCREMENT      UPDATE_LOOP: DO I_UIID = 1,N_UPDATES ! Update counters inside the radiation routine          ANGLE_INC_COUNTER = MOD(ANGLE_INC_COUNTER,ANGLE_INCREMENT) + 1         IF (WIDE_BAND_MODEL) THEN            UIID(:,:,:,IBND) = 0._EB         ELSE            UIID(:,:,:,ANGLE_INC_COUNTER) = 0._EB         ENDIF ! Set the bounds and increment for the angleloop. Step downdard because in cylindrical case the Nth angle ! boundary condition comes from (N+1)th angle.          NSTART    = NRA - ANGLE_INC_COUNTER + 1         NEND      = 1         NSTEP     = -ANGLE_INCREMENT         IL(:,:,:) = BBFA*RPI_SIGMA*TMPA4         ANGLE_LOOP: DO N = NSTART,NEND,NSTEP  ! Sweep through control angles ! Boundary conditions: Intensities leaving the boundaries.             WALL_LOOP1: DO IW=1,NWC               IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_LOOP1               IOR = IJKW(4,IW)               IF (DLN(IOR,N) < 0._EB) CYCLE WALL_LOOP1               II  = IJKW(1,IW)               JJ  = IJKW(2,IW)               KK  = IJKW(3,IW)               IF (.NOT.TWO_D .OR. ABS(IOR)/=2) THEN                  SELECT CASE (BOUNDARY_TYPE(IW))                     CASE (OPEN_BOUNDARY)                            IL(II,JJ,KK) = BBFA*RPI_SIGMA*TMPA4