MODULE RAD ! Radiation heat transfer USE PRECISION_PARAMETERSUSE GLOBAL_CONSTANTSUSE MESH_VARIABLESUSE RADCONSIMPLICIT NONEPRIVATECHARACTER(255), PARAMETER :: radiid='$Id: radi.f90 712 2007-09-28 20:20:09Z drjfloyd $'CHARACTER(255), PARAMETER :: radirev='$Revision: 712 $'CHARACTER(255), PARAMETER :: radidate='$Date: 2007-09-28 16:20:09 -0400 (Fri, 28 Sep 2007) $'PUBLIC INIT_RADIATION,COMPUTE_RADIATION,NSB,NRA,UIIDIM,NRT,RSA,NRP,GET_REV_radi CONTAINS  SUBROUTINE INIT_RADIATIONUSE MEMORY_FUNCTIONS, ONLY : CHKMEMERRUSE MIEVUSE RADCALVREAL(EB) :: THETAUP,THETALOW,PHIUP,PHILOW,F_THETA,MW_RADCAL,PLANCK_C2,KSI,LT, &            ZZ,RCRHO,YY,BBF,AP0,AMEAN,MTOT,XLENG,YLENG,ZLENGINTEGER  :: N,I,J,K,IZERO,NN,NI,II,JJ,IIM,JJM,IBND,NS,I_RADCAL,IZTYPE(SPECIES_TYPE), POINTER :: SSTYPE(REACTION_TYPE), POINTER :: RNTYPE(PARTICLE_CLASS_TYPE), POINTER :: PC ! Determine the number of polar angles (theta) NRA = NUMBER_RADIATION_ANGLESIF (CYLINDRICAL) THEN   NRT = NINT(SQRT(REAL(NRA)))ELSEIF (TWO_D) THEN   NRT = 1ELSE   NRT = 2*NINT(0.5_EB*1.17*REAL(NRA)**(1._EB/2.26))ENDIF       ALLOCATE(NRP(1:NRT),STAT=IZERO)CALL ChkMemErr('INIT','NRP',IZERO) ! Determine number of azimuthal angles (phi) N = 0DO I=1,NRT   IF (CYLINDRICAL) THEN      NRP(I) = NINT(REAL(NRA)/(REAL(NRT)))      ELSEIF (TWO_D) THEN      NRP(I) = 4*NINT(0.25_EB*REAL(NRA))      ELSE      THETALOW = PI*REAL(I-1)/REAL(NRT)      THETAUP  = PI*REAL(I)/REAL(NRT)      NRP(I) = NINT(0.5_EB*REAL(NRA)*(COS(THETALOW)-COS(THETAUP)))      NRP(I) = MAX(4,NRP(I))      NRP(I) = 4*NINT(0.25_EB*REAL(NRP(I)))      ENDIF   N = N + NRP(I)ENDDONRA = NNUMBER_RADIATION_ANGLES = NRA ! Set the opening angle of the cylindrical geometry equal to the azimuthal angle IF (CYLINDRICAL) DPHI0 = PI/REAL(NRP(1)) NSB = NUMBER_SPECTRAL_BANDS ALLOCATE(RSA(1:NRA),STAT=IZERO)CALL ChkMemErr('INIT','RSA',IZERO)ALLOCATE(DLX(1:NRA),STAT=IZERO)CALL ChkMemErr('INIT','DLX',IZERO)ALLOCATE(DLY(1:NRA),STAT=IZERO)CALL ChkMemErr('INIT','DLY',IZERO)ALLOCATE(DLZ(1:NRA),STAT=IZERO)CALL ChkMemErr('INIT','DLZ',IZERO)IF (CYLINDRICAL) THEN   ALLOCATE(DLB(1:NRA),STAT=IZERO)   CALL ChkMemErr('INIT','DLB',IZERO)ENDIFALLOCATE(DLN(-3:3,1:NRA),STAT=IZERO)CALL ChkMemErr('INIT','DLN',IZERO)ALLOCATE(DLM(1:NRA,3),STAT=IZERO)CALL ChkMemErr('INIT','DLM',IZERO) ! Determine mean direction normals and sweeping orders N = 0DO I=1,NRT   DO J=1,NRP(I)      N = N + 1      THETALOW  = PI*REAL(I-1)/REAL(NRT)      THETAUP   = PI*REAL(I)/REAL(NRT)      F_THETA   = 0.5_EB*(THETAUP-THETALOW  - COS(THETAUP)*SIN(THETAUP) + COS(THETALOW)*SIN(THETALOW))      IF (CYLINDRICAL) THEN         PHILOW = PI*REAL(J-1)/REAL(NRP(I))         PHIUP  = PI*REAL(J)/REAL(NRP(I))      ELSEIF (TWO_D) THEN         PHILOW = TWOPI*REAL(J-1)/REAL(NRP(I)) + PIO2         PHIUP  = TWOPI*REAL(J)/REAL(NRP(I))   + PIO2      ELSE         PHILOW = TWOPI*REAL(J-1)/REAL(NRP(I))         PHIUP  = TWOPI*REAL(J)/REAL(NRP(I))      ENDIF      RSA(N) = (PHIUP-PHILOW)*(COS(THETALOW)-COS(THETAUP))      IF (CYLINDRICAL) THEN         DLX(N) = 2._EB*SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW)) *F_THETA         DLY(N) =  (-SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW))  +COS(DPHI0/2.)*(COS(PHILOW)-COS(PHIUP)))*F_THETA         DLB(N) =  (-SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW))  -COS(DPHI0/2.)*(COS(PHILOW)-COS(PHIUP)))*F_THETA         DLZ(N)    = 0.5_EB*(PHIUP-PHILOW)   * ((SIN(THETAUP))**2-(SIN(THETALOW))**2)      ELSEIF (TWO_D) THEN         DLX(N) = (SIN(PHIUP)-SIN(PHILOW))*F_THETA         DLY(N) = 0._EB         DLZ(N) = (COS(PHILOW)-COS(PHIUP))*F_THETA      ELSE         DLX(N) = (SIN(PHIUP)-SIN(PHILOW))*F_THETA         DLY(N) = (COS(PHILOW)-COS(PHIUP))*F_THETA         DLZ(N)    = 0.5_EB*(PHIUP-PHILOW)      * ((SIN(THETAUP))**2-(SIN(THETALOW))**2)      ENDIF   ENDDOENDDO ! Set (wall normal)*(angle vector) value DO N = 1,NRA   DLN(-1,N) = -DLX(N)   DLN( 1,N) =  DLX(N)   DLN(-2,N) = -DLY(N)   DLN( 2,N) =  DLY(N)   DLN(-3,N) = -DLZ(N)   DLN( 3,N) =  DLZ(N)ENDDO ! Calculate mirroring matrix N = 0DO I=1,NRT   DO J=1,NRP(I)      N = N + 1      DO K=1,3         IF (TWO_D .AND. .NOT.CYLINDRICAL) THEN            SELECT CASE(K)               CASE(1)             ! X-surfaces                  IIM = 1                  JJM = NRP(I) - J + 1               CASE(2)             ! Y-surfaces                  IIM = 1                  JJM = J               CASE(3)             ! Z-surfaces                  IIM = 1                  JJM = NRP(I)/2 - J + 1            END SELECT            JJM = MODULO(JJM,NRP(I))            IF (JJM==0) JJM = NRP(I)         ELSE            SELECT CASE(K)               CASE(1)             ! X-surfaces                  IIM = I                  JJM = NRP(I)/2 - J + 1               CASE(2)             ! Y-surfaces                  IIM = I                  JJM = NRP(I) - J + 1               CASE(3)             ! Z-surfaces                  IIM = NRT - I + 1                  JJM = J            END SELECT            IIM = MODULO(IIM,NRT)            JJM = MODULO(JJM,NRP(I))            IF (IIM==0) IIM = NRT            IF (JJM==0) JJM = NRP(I)         ENDIF          NN = 0         DO II = 1,IIM            DO JJ = 1,NRP(II)               NN = NN + 1               IF ((II==IIM).AND.(JJ==JJM)) NI = NN            ENDDO         ENDDO         DLM(N,K) = NI      ENDDO   ENDDOENDDO !-----------------------------------------------------!!            Spectral information!!-----------------------------------------------------INIT_WIDE_BAND: IF (WIDE_BAND_MODEL) THEN ! Fraction of blackbody emission in a wavelength interval    PLANCK_C2 = 14387.69_EB            ! micron.K   NLAMBDAT  = 4000_EB   LTSTEP    = 25.0_EB           ! maximum LAMBDA*T = NLANBDAT*LTSTEP   ALLOCATE(BBFRAC(0:NLAMBDAT),STAT=IZERO)   CALL ChkMemErr('INIT','BBFRAC',IZERO)   BBFRAC = 0._EB   LT = 0._EB   DO I = 1,NLAMBDAT      LT =  LT + LTSTEP      KSI = PLANCK_C2/LT      DO J = 1,50         BBFRAC(I) = BBFRAC(I) + EXP(-KSI*REAL(J))/REAL(J) * (KSI**3 + 3.*KSI**2/REAL(J) + 6.*KSI/REAL(J)**2 + 6./REAL(J)**3)      ENDDO   ENDDO   BBFRAC =  BBFRAC * 15._EB/PI**4 ! Define band limit wave lengths in micrometers    ALLOCATE(WL_LOW(1:NSB),STAT=IZERO)   CALL ChkMemErr('INIT','WL_LOW',IZERO)   ALLOCATE(WL_HIGH(1:NSB),STAT=IZERO)   CALL ChkMemErr('INIT','WL_HIGH',IZERO)   IF (CH4_BANDS) THEN      WL_LOW(1:NSB) =(/1.00_EB, 2.63_EB, 2.94_EB, 3.57_EB, 4.17_EB, 4.6_EB, 7.00_EB, 8.62_EB, 10.0_EB /)      WL_HIGH(1:NSB)=(/2.63_EB, 2.94_EB, 3.57_EB, 4.17_EB, 4.60_EB, 7.0_EB, 8.62_EB, 10.0_EB, 200._EB /)    ELSE      WL_LOW(1:NSB) =(/1.00_EB, 2.63_EB, 2.94_EB, 4.17_EB, 4.6_EB, 10.0_EB /)      WL_HIGH(1:NSB)=(/2.63_EB, 2.94_EB, 4.17_EB, 4.6_EB, 10.0_EB, 200.0_EB /)   ENDIF ENDIF INIT_WIDE_BAND !----------------------------------------------------------------------------!!     Tables for gas phase absorption coefficient!!     CONTROLLING PROGRAM FOR SUBROUTINE "RADCAL", A NARROW-BAND!     MODEL FOR CALCULATING SPECTRAL INTENSITY (W/M-2/SR/MICRON) AND!     SPECTRAL TRANSMITTANCE VERSUS WAVELENGTH (MICRONS) IN A NONISO-!     THERMAL, VARIABLE COMPOSITION  MIXTURE OF CO2, H2O, CO, N2, O2,!     CH4, AND SOOT. FOR A HOMOGENEOUS PATH, THE PROGRAM ALSO COMPUTES!     THE PLANCK-MEAN ABSORPTION COEF, AP0, THE INCIDENT-MEAN ABSORPTION!     COEFFICIENT, AIWALL, AND THE EFFECTIVE-MEAN ABSORPTION COEFFICIENT,!     AMEAN, ALL IN UNITS OF INVERSE METERS.!!     INPUT PARAMETERS:!          NPT=NUMBER OF HOMOGENEOUS ELEMENTS!          DD(J)=THICKNESS OF J TH ELEMENT, M!          RCT(J)=TEMPERATURE OF J TH ELEMENT, K.!          P(I,J)=PARTIAL PRESSURE OF GASEOUS COMPONENTS, kPa:!                  I   GASEOUS SPECIES!                  1        CO2!                  2        H2O!                  3        CH4!                  4        CO!                  5        O2!                  6        N2!          SVF(J)=SOOT VOLUME FRACTION OF J TH ELEMENT!          OMMIN=MINIMUM WAVE NUMBER IN SPECTRUM, CM-1.!          OMMAX=MAXIMUM WAVE NUMBER IN SPECTRUM, CM-1.!!------------------------------------------------------------------------- CALL RCALLOC  ! Allocate arrays for RadCal ! 20% of mean beam length, Eq 8-51, Holman, 7th Ed. Heat Transfer. Length = 3.6*Volume/Area XLENG = MESHES(1)%XF-MESHES(1)%XSYLENG = MESHES(1)%YF-MESHES(1)%YSZLENG = MESHES(1)%ZF-MESHES(1)%ZSIF (PATH_LENGTH < 0.0_EB) THEN  ! default was -1.0   IF (TWO_D) THEN ! calculate based on the geometry      PATH_LENGTH = MIN( 10._EB , 0.1_EB*3.6_EB*XLENG*ZLENG/(XLENG+ZLENG) )   ELSE      PATH_LENGTH = MIN( 10._EB , 0.1_EB*3.6_EB*XLENG*YLENG*ZLENG/(XLENG*YLENG+XLENG*ZLENG+YLENG*ZLENG) )   ENDIFENDIFDD(1) = MAX(PATH_LENGTH,1.0E-4_EB) ! Using RadCal, create look-up tables for the absorption coefficients for all gas species, mixture fraction or aerosolsSPECIES_LOOP: DO NS=1,N_SPECIES   SS => SPECIES(NS)   IF (.NOT. SS%ABSORBING) CYCLE SPECIES_LOOP    GAS_TYPE: SELECT CASE (SS%MODE)        CASE (MIXTURE_FRACTION_SPECIES) GAS_TYPE          RN => REACTION(SS%REAC_INDEX)         SS%NKAP_TEMP = 21         SS%NKAP_MASS = 40         SS%MAXMASS   = RN%Z_F         ALLOCATE(SS%KAPPA(0:SS%NKAP_MASS,0:SS%NKAP_TEMP,1:NSB),STAT=IZERO)         CALL ChkMemErr('INIT','SS%KAPPA',IZERO)         BAND_LOOP_MF: DO IBND = 1,NSB            IF (NSB>1) THEN               OMMIN = REAL(NINT(1.E4/WL_HIGH(IBND)),EB)               OMMAX = REAL(NINT(1.E4/WL_LOW(IBND)),EB)            ELSE               OMMIN = 50.               OMMAX = 10000.            ENDIF            CALL INIT_RADCAL                  T_LOOP_MF: DO K = 0,SS%NKAP_TEMP               RCT(1) = 300._EB + K*(2400._EB-300._EB)/SS%NKAP_TEMP               Z_LOOP_MF: DO J=0,SS%NKAP_MASS                  ZZ = IZ2ZZ(J,RN%Z_F,SS%NKAP_MASS)                  IZ = MIN(10000,MAX(0,NINT(ZZ*10000._EB)))                  RCRHO = SS%MW_MF(IZ)*P_INF/(R0*RCT(1))                  MTOT = SS%Y_MF(IZ,FUEL_INDEX)/RN%MW_FUEL + SS%Y_MF(IZ,O2_INDEX)/MW_O2   + SS%Y_MF(IZ,N2_INDEX)/MW_N2 +  &                         SS%Y_MF(IZ,H2O_INDEX)/MW_H2O   + SS%Y_MF(IZ,CO2_INDEX)/MW_CO2 + SS%Y_MF(IZ,CO_INDEX)/MW_CO                  SPECIE(1) = SS%Y_MF(IZ,CO2_INDEX)                  SPECIE(2) = SS%Y_MF(IZ,H2O_INDEX)                  SPECIE(3) = SS%Y_MF(IZ,FUEL_INDEX)                  SPECIE(4) = SS%Y_MF(IZ,CO_INDEX)                  SPECIE(5) = SS%Y_MF(IZ,SOOT_INDEX)*RCRHO/RHO_SOOT                  P(1,1) = SS%Y_MF(IZ,CO2_INDEX)/MW_CO2/MTOT                  P(2,1) = SS%Y_MF(IZ,H2O_INDEX)/MW_H2O/MTOT                  P(3,1) = SS%Y_MF(IZ,FUEL_INDEX)/RN%MW_FUEL/MTOT                  P(4,1) = SS%Y_MF(IZ,CO_INDEX)/MW_CO/MTOT                  P(5,1) = SS%Y_MF(IZ,O2_INDEX)/MW_O2/MTOT                  P(6,1) = SS%Y_MF(IZ,N2_INDEX)/MW_N2/MTOT                  SVF(1) = SS%Y_MF(IZ,SOOT_INDEX)*RCRHO/RHO_SOOT                  CALL RADCAL(AMEAN,AP0)                  IF (NSB==1 .AND. PATH_LENGTH > 0.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(IBND),WL_HIGH(IBND),RCT(1))                     ENDIF                     SS%KAPPA(J,K,IBND) = AP0/BBF                  ENDIF               ENDDO Z_LOOP_MF            ENDDO T_LOOP_MF         ENDDO BAND_LOOP_MF      CASE (GAS_SPECIES) GAS_TYPE         SS%NKAP_TEMP = 21         SS%NKAP_MASS = 200         SS%MAXMASS=1._EB         ALLOCATE(SS%KAPPA(0:SS%NKAP_MASS,0:SS%NKAP_TEMP,1:NSB),STAT=IZERO)         CALL ChkMemErr('RADI','KAPPA',IZERO)         IF (NS==I_CO2) THEN            I_RADCAL = 1            MW_RADCAL = MW_CO2         ELSEIF (NS==I_CO) THEN            I_RADCAL = 4            MW_RADCAL = MW_CO         ELSEIF (NS==I_WATER) THEN            I_RADCAL = 2            MW_RADCAL = MW_H2O         ELSE            IF (SS%ABSORBING) THEN               I_RADCAL = 3               MW_RADCAL = 16            ELSE               I_RADCAL = 5