irad.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

MODULE RADCALV

! Module wrapper for RadCal subroutine

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS, ONLY: PI,SQRTPI
IMPLICIT NONE

PRIVATE
CHARACTER(255), PARAMETER :: iradid='$Id: irad.f90 710 2007-09-28 19:52:12Z drjfloyd $'
CHARACTER(255), PARAMETER :: iradrev='$Revision: 710 $'
CHARACTER(255), PARAMETER :: iraddate='$Date: 2007-09-28 15:52:12 -0400 (Fri, 28 Sep 2007) $'

PUBLIC OMMAX,OMMIN,DD,SPECIE,SVF,PLANCK,P,RCT,RCALLOC,INIT_RADCAL,RADCAL,RCDEALLOC,GET_REV_irad
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: P,X,UUU,GC,GAMMA,SD15, SD,SD7,SD3
REAL(EB), ALLOCATABLE, DIMENSION(:) :: SVF,XPART,RCT,SPECIE,TAU,TAUS,AC,AD,QW,TTAU,XTOT,XT,XSTAR,AB,PKPA,AMBDA,ATOT,BCNT,DD
REAL(EB) OMMIN,OMMAX,TWALL
INTEGER NPT,NPRINT,NOM

CONTAINS

SUBROUTINE GET_REV_irad(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') iradrev(INDEX(iradrev,':')+1:LEN_TRIM(iradrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') iraddate

END SUBROUTINE GET_REV_irad

SUBROUTINE INIT_RADCAL

NPT=1
TWALL=0._EB
!      TWALL=293.15
IF(OMMAX<1100._EB) THEN 
   NOM=INT((OMMAX-OMMIN)/5._EB)
ELSEIF(OMMIN>5000._EB) THEN
   NOM=INT((OMMAX-OMMIN)/50._EB)
ELSEIF(OMMIN<1100._EB.AND.OMMAX>5000._EB) THEN
   NOM=INT((1100._EB-OMMIN)/5._EB)+INT((5000._EB-1100._EB)/25._EB) +INT((OMMAX-5000._EB)/50._EB)
ELSEIF(OMMIN<1100._EB) THEN
   NOM=INT((1100._EB-OMMIN)/5._EB)+INT((OMMAX-1100._EB)/25._EB)
ELSEIF(OMMAX>5000._EB) THEN
   NOM=INT((5000._EB-OMMIN)/25._EB)+INT((OMMAX-5000._EB)/50._EB)
ELSE
   NOM=INT((OMMAX-OMMIN)/25._EB)
ENDIF
NPRINT=1

END SUBROUTINE INIT_RADCAL

!**************************************************************************

SUBROUTINE RADCAL(AMEAN,AP0)
!
INTEGER I,II,J,K,L,KK,NM,N,MM,KMAX,LMAX,KMIN
REAL (EB) DOM,ABGAS,PTOT,TEMP,UK,GKD,GKDD,XD,YD,XX,ENN,ARG,ARGNEW, &
           RSL,RSS,ABLONG,ABSHRT,ABIL,ABIS,OMEGA,WL,DAMBDA,AP0, &
           SDWEAK,GDINV,GDDINV,YC,Y,AIWALL,AMEAN,XC,AOM,Q, &
           LTERM

!  [NOTE: THE TOTAL INTENSITY CALCULATED IS THAT WHICH LEAVES INTERVAL J=1.
!  P(I,J) IS PARTIAL PRESSURE, ATM, OF SPECIES I IN  INTERVAL J.
!  I=1,2,3,4,5, OR 6 IMPLIES SPECIES IS CO2, H2O, CH4, CO, O2, OR N2, RESP.]

DOM=5.0_EB
OMEGA=OMMIN-DOM
NM=NOM-1

!     LOOP 1000 COMPUTES EACH SPECTRAL CONTRIBUTION
!     *********************************************
L1000: DO KK=1,NOM
   OMEGA=OMEGA+DOM
   IF(OMEGA>1100._EB) OMEGA=OMEGA+20._EB
   IF(OMEGA>5000._EB) OMEGA=OMEGA+25._EB
   AMBDA(KK)=10000._EB/OMEGA
   ABGAS=0._EB
!
!          LOOP 200 COMPUTES THE CONTRIBUTION OF EACH SPECIES TO TAU
!	   *********************************************************
   L200: DO I=1,4
!     IF SPECIE(I) IS SET TO 0., THAT PARTICULAR RADIATING SPECIES IS
!     NOT PRESENT.  THE SPECIES CONSIDERED ARE
!          I   SPECIES
!          1     CO2
!          2     H2O
!          3     CH4
!          4     CO
!          5     PARTICULATES
      IF(SPECIE(I)==0._EB) CYCLE L200
!
!               LOOP 100 IS FOR EACH ELEMENT ALONG PATH
!               ***************************************
      L100: DO J=1,NPT
!     (CALCULATION PROCEEDS IN ACCORDANCE WITH THE SLG MODEL, TABLE 5-18
!     IN NASA SP-3080._EB)
         IF(KK<=1) THEN
            UUU(I,J)=273./RCT(J)*P(I,J)*100.*DD(J)
            GC(I,J)=0._EB
            PTOT=0._EB
            DO II=1,6
               PTOT=P(II,J)+PTOT
               GC(I,J)=GC(I,J)+GAMMA(I,II)*P(II,J)*(273./RCT(J))**.5
            ENDDO
            GC(I,J)=GC(I,J)+GAMMA(I,7)*P(I,J)*273./RCT(J)
         ENDIF
         IF(P(I,J)==0._EB) THEN
            IF(J>1) THEN
               XSTAR(J)=XSTAR(J-1)
               AC(J)=AC(J-1)
               AD(J)=AD(J-1)
            ELSE
               XSTAR(1)=1.E-34_EB
               AC(1)=1._EB
               AD(1)=1._EB
            ENDIF
            X(I,J)=XSTAR(J)
         ELSE
            TEMP=RCT(J)
            SELECT CASE (I)
            CASE (1)
            CALL CO2(OMEGA,TEMP,GC(1,J),SDWEAK,GDINV,GDDINV)
            CASE (2)
            CALL H2O(OMEGA,TEMP,GC(2,J),SDWEAK,GDINV,GDDINV)
            CASE (3)
            CALL FUEL(OMEGA,TEMP,P(3,J),PTOT,GC(3,J),SDWEAK,GDINV,GDDINV)
            CASE (4)
            CALL CO(OMEGA,TEMP,GC(4,J),SDWEAK,GDINV,GDDINV)   
            END SELECT
            UK=SDWEAK*UUU(I,J)
            IF(J==1) THEN
               XSTAR(1)=UK+1.E-34_EB
               ABGAS=UK/DD(1)+ABGAS
               AD(1)=GDDINV
               AC(1)=GDINV
            ELSE
               GKD=UK*GDINV
               GKDD=UK*GDDINV
               XSTAR(J)=XSTAR(J-1)+UK
               AD(J)=(XSTAR(J-1)*AD(J-1)+GKDD)/XSTAR(J)
               AC(J)=(XSTAR(J-1)*AC(J-1)+GKD)/XSTAR(J)
            ENDIF
            IF(XSTAR(J)>=1.E-6_EB) THEN
               XD=1.7_EB*AD(J)*(DLOG(1._EB+(XSTAR(J)/1.7_EB/AD(J))**2))**.5_EB
               YD=1._EB-(XD/XSTAR(J))**2
               XC=XSTAR(J)/(1._EB+XSTAR(J)/4._EB/AC(J))**.5_EB
!
!   THE FOLLOWING LOOP COMPUTES THE OPTICAL THICKNESS, XC, FOR METHANE USING
!   THE GODSON EQUATION AND AN APPROXIMATION TO THE LADENBERG-REICHE
!   FUNCTION AS RECOMMENDED BY BROSMER AND TIEN (JQSRT 33,P 521).  THE
!   ERROR FUNCTION IS FOUND FROM ITS SERIES EXPANSION.
!
               IF((I==3._EB).AND.(XC<=10)) THEN
                  AOM=XC
                  XX=.5_EB*SQRTPI*XC
                  IF(XX<=3._EB) THEN
                     ENN=1._EB
                     DO N=1,30
                        ENN=ENN*REAL(N,EB)
                        MM=2*N+1
                        ARG=1.128379_EB*(-1._EB)**N*((.88622693_EB*XC)**MM)/(REAL(MM,EB)*ENN)
                        ARGNEW=ARG+AOM
!     IF(ABS(ARG/ARGNEW)<.000001)N=30
                        AOM=ARGNEW
                     ENDDO
                  ELSE
                     AOM=1._EB-EXP(-XX**2)/(SQRTPI*XX)
                  ENDIF
                  IF(AOM>=1._EB)AOM=.9999999_EB
                  XC=-DLOG(1._EB-AOM)
               ENDIF
     
               YC=1._EB-(XC/XSTAR(J))**2
               Y=MAX(1._EB/YC**2+1._EB/YD**2-1._EB,1._EB)
               X(I,J)=XSTAR(J)*((1._EB-(Y**(-.5_EB)))**.5_EB)
            ELSE
               X(I,J)=XSTAR(J)
            ENDIF
         ENDIF
      ENDDO L100

   ENDDO L200

!  DETERMINE OPTICAL DEPTH OF SOOT

   IF(SPECIE(5)==0._EB) THEN
      DO J=1,NPT
         XPART(J)=0._EB
      ENDDO
   ELSE
      CALL POD(OMEGA)
   ENDIF
AB(KK)=ABGAS+XPART(1)/DD(1)
     
!     EVALUATE THE COMBINED SPECTRAL TRANSMITTANCE AND RADIANCE
!     *********************************************************
   L500: DO J=1,NPT
      XTOT(J)=0._EB
      DO I=1,4
         IF(SPECIE(I)==0._EB) X(I,J)=0._EB
         XTOT(J)=X(I,J)+XTOT(J)
      ENDDO
      XTOT(J)=XTOT(J)+XPART(J)
      IF(XTOT(J)>=99._EB) THEN
         TAU(J)=0._EB
      ELSE
         TAU(J)=EXP(-XTOT(J))
      ENDIF
      IF(J==1) THEN
         QW(KK)=-(TAU(1)-1._EB)*PLANCK(RCT(1),AMBDA(KK))
      ELSE
         QW(KK)=QW(KK)-(TAU(J)-TAU(J-1))*PLANCK(RCT(J),AMBDA(KK))
      ENDIF
   ENDDO L500
   XT(KK)=XTOT(NPT)
   TTAU(KK)=TAU(NPT)
   QW(KK)=QW(KK)+TTAU(KK)*PLANCK(TWALL,AMBDA(KK))
ENDDO L1000
     
!     INTEGRATE THE RADIANCE OVER THE SPECTRUM
     
Q=QW(1)*(AMBDA(1)-AMBDA(2))
DO KK=2,NM
   Q=Q+QW(KK)*(AMBDA(KK-1)-AMBDA(KK+1))/2._EB
ENDDO
Q=Q+QW(NOM)*(AMBDA(NOM-1)-AMBDA(NOM))
     
!     DETERMINE SOOT RADIANCE FOR SHORT AND LONG WAVELENGTHS.
     
RSL=0._EB
RSS=0._EB
ABLONG=0._EB
ABSHRT=0._EB
ABIL=0._EB
ABIS=0._EB
   IF(.NOT.(SPECIE(5)==0._EB .AND. TWALL==0._EB)) THEN
   KMAX=OMMIN/5*5
   DO KK=5,KMAX,5
      OMEGA=FLOAT(KK)
      WL=10000._EB/OMEGA
      DAMBDA=10000._EB/(OMEGA-2.5_EB)-10000._EB/(OMEGA+2.5_EB)
      CALL POD(OMEGA)
   
      DO J=1,NPT
         IF(XPART(J)>=33._EB) THEN
            TAUS(J)=0._EB
         ELSE
            TAUS(J)=EXP(-XPART(J))
         ENDIF
         IF(J==1) THEN
            RSL=RSL-(TAUS(1)-1._EB)*PLANCK(RCT(1),WL)*DAMBDA
            ABLONG=ABLONG+XPART(1)/DD(1)*PLANCK(RCT(1),WL)*DAMBDA*5.5411E7_EB/(RCT(1))**4
            ABIL=ABIL+XPART(1)/DD(1)*PLANCK(TWALL,WL)*DAMBDA*5.5411E7_EB /(TWALL+.000001_EB)**4
         ELSE
            RSL=RSL-(TAUS(J)-TAUS(J-1))*PLANCK(RCT(J),WL)*DAMBDA
         ENDIF
      ENDDO

      RSL=RSL+TAUS(NPT)*PLANCK(TWALL,WL)*DAMBDA
   ENDDO
   KMIN=OMMAX/100*100
   DO KK=KMIN,25000,100
      OMEGA=FLOAT(KK)
      WL=10000._EB/OMEGA
      DAMBDA=10000._EB/(OMEGA-50._EB)-10000._EB/(OMEGA+50._EB)
      CALL POD(OMEGA)
      DO J=1,NPT
         IF(XPART(J)>=33._EB) THEN
            TAUS(J)=0._EB
         ELSE
            TAUS(J)=EXP(-XPART(J))
         ENDIF
         IF(J/=1) THEN
            RSS=RSS-(TAUS(J)-TAUS(J-1))*PLANCK(RCT(J),WL)*DAMBDA
         ELSE
            RSS=RSS-(TAUS(1)-1._EB)*PLANCK(RCT(1),WL)*DAMBDA
            ABSHRT=ABSHRT+XPART(1)/DD(1)*PLANCK(RCT(1),WL)*DAMBDA*5.5411E7_EB/(RCT(1))**4
            ABIS=ABIS+XPART(1)/DD(1)*PLANCK(TWALL,WL)*DAMBDA*5.5411E7_EB/(TWALL+.000001_EB)**4
         ENDIF
      ENDDO
      RSS=RSS+TAUS(NPT)*PLANCK(TWALL,WL)*DAMBDA
   ENDDO
ENDIF
Q=Q+RSS+RSL


IF (NPRINT/=0) THEN
   IF(NPRINT==1) THEN
      NPRINT=2
      DO J=1,NPT
         DO I=1,6
            PKPA(I)=P(I,J)*101._EB
         ENDDO
      ENDDO