irad.f90

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

   ENDIF
   LMAX=NOM/2
   IF(LMAX*2<NOM) LMAX=LMAX+1
   DO L=1,LMAX
      K=NOM-LMAX+1-L
      J=K+LMAX
      IF(K<1) K=1
   ENDDO
ENDIF

!     THE FOLLOWING SECTION COMPUTES THE MEAN ABSORPTION COEFFICIENTS
!     IF THE SYSTEM IS HOMOGENEOUS (IE., NPT=1).

IF(NPT==1) THEN
   NM=NOM-1
   AIWALL=AB(1)*(AMBDA(1)-AMBDA(2))/2._EB*PLANCK(TWALL,AMBDA(1))
   AP0=AB(1)*(AMBDA(1)-AMBDA(2))/2._EB*PLANCK(RCT(1),AMBDA(1))
   DO KK=2,NM
      AIWALL=AIWALL+AB(KK)*(AMBDA(KK-1)-AMBDA(KK+1))/2._EB *PLANCK(TWALL,AMBDA(KK))
      AP0=AP0+AB(KK)*(AMBDA(KK-1)-AMBDA(KK+1))/2._EB*PLANCK(RCT(1),AMBDA(KK))
   ENDDO
   AP0=(AP0+AB(NOM)*(AMBDA(NM)-AMBDA(NOM))/2._EB *PLANCK(RCT(1),AMBDA(NOM)))*5.5411E7_EB/RCT(1)**4
   IF(TWALL==RCT(1).OR.TWALL==0._EB) THEN
      AIWALL=AP0
      LTERM = MAX(1E-10_EB,(5.5411E7_EB*Q-RCT(1)**4)/(-RCT(1)**4))
!      LTERM = (1._EB-5.5411E7*Q/RCT(1)**4)
      AMEAN=-1._EB/DD(1)*DLOG(LTERM)
   ELSE
      AIWALL=(AIWALL+AB(NOM)*(AMBDA(NM)-AMBDA(NOM))/2._EB*PLANCK(TWALL,AMBDA(NOM)))*5.5411E7_EB/TWALL**4
      AMEAN=-1._EB/DD(1)*DLOG((5.5411E7_EB*Q-RCT(1)**4)/(TWALL**4-RCT(1)**4))
   ENDIF
ENDIF

END SUBROUTINE RADCAL

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

SUBROUTINE CO2(OMEGA,TEMP,GC1,SDWEAK,GDINV,GDDINV)

INTEGER I,J,K,L
REAL(EB) OMEGA,TEMP,GC1,SDWEAK,GDINV,GDDINV,AA,BB,CC,QQ,EE,FF,GG, &
         SMINUS,SPLUS,SDSTRG,GD,OM1,OM2,OM3,T0,Q2,BE,COM1, &
         COM2,COM3,X13,X23,X33,XBAR,OM12,ALPHA,OMPRIM,V3,GAM, &
         OMVV3,DELTA,V,OMVBAR,F1,F2,UNFLO1,UNFLO2,UNFLO3,TEST, &
         VBAR1,OMA,OMB,TTEMP,TT,T1,TW,WW,TEMP1,TEMP2,TEMP3,WM, &
         DINV,A,B,D,G,W1,DINV1,DINV2,DINV3

IF(OMEGA>5725._EB) THEN
   SDWEAK=0._EB
   GDINV=1._EB
   GDDINV=1._EB
ELSE
   WM=44._EB
   GD=5.94E-6_EB*OMEGA*(TEMP/(273._EB*WM))**.5_EB
   IF(OMEGA>4550._EB) THEN
!CONTRIBUTION TO 2.0 MICRON BAND FROM (000)-(041),(000)-(121),AND (000)
!     -(201) TRANS.
      OM1=1354.91_EB
      OM2=673.0_EB
      OM3=2396.49_EB
      BCNT(1)=4860.5_EB
      BCNT(2)=4983.5_EB
      BCNT(3)=5109.0_EB
      T0=300._EB
      Q2=1.4388_EB
      BE=0.391635_EB
      COM1=4._EB*OM2+OM3
      COM2=OM1+2._EB*OM2+OM3
      COM3=2._EB*OM1+OM3
      ATOT(3)=0.426_EB*T0/TEMP*(1._EB-EXP(-Q2*COM3/TEMP/  (1._EB-EXP(-Q2*OM1/TEMP))))**2/(1._EB-EXP(-Q2*OM3/TEMP))
      ATOT(2)=1.01_EB *T0/TEMP*(1._EB-EXP(-Q2*COM2/TEMP))/(1._EB-EXP(-Q2*OM1/TEMP))/(1._EB-EXP(-Q2*OM2/TEMP))**2/ &
              (1._EB-EXP(-Q2*OM3/TEMP))
      ATOT(1)=0.272_EB*T0/TEMP*(1._EB-EXP(-Q2*COM1/TEMP))/(1._EB-EXP(-Q2*OM2/TEMP))**4/(1._EB-EXP(-Q2*OM3/TEMP))
      SDWEAK=0.0_EB
      DO K=1,3
         SDWEAK=SDWEAK+ATOT(K)*Q2/(4._EB*BE*TEMP)*ABS(OMEGA-BCNT(K))*EXP(-Q2/(4._EB*BE*TEMP)*(OMEGA-BCNT(K))**2)
      ENDDO
      DINV=1._EB/(4._EB*BE)
      GDINV=GC1*DINV
      GDDINV=GD*DINV
!***EXPRESS S/D AT STP, AS IS IN NASA SP-3080
      SDWEAK=SDWEAK*TEMP/273._EB
   ELSEIF((OMEGA<=4550_EB).AND.(OMEGA>3800._EB)) THEN
      SDWEAK=0._EB
      GDINV=1._EB
      GDDINV=1._EB
   ELSEIF((OMEGA<=3800._EB).AND.(OMEGA>3050._EB)) THEN
      B=.391635_EB
      A=.0030875_EB
      X13=-19.37_EB
      X23=-12.53_EB
      X33=-12.63_EB
      OM1=1354.91_EB
      OM2=673._EB
      OM3=2396.49_EB
      T0=300._EB
      Q2=1.4388_EB
      XBAR=.5_EB*(.5_EB*X13+X23)
      OM12=.5_EB*(.5_EB*OM1+OM2)
      SDWEAK=0._EB
      SDSTRG=0._EB
      IF(OMEGA<=2395._EB) THEN
         ALPHA=2700._EB
         OMPRIM=OM3
         AA=ALPHA*B*Q2/(A*(1._EB-EXP(-OM3*Q2/T0))*(1._EB-EXP(-OM12*Q2 /T0))**3*(1._EB+EXP(-OM12*Q2/T0))*(1._EB-EXP(-OMPRIM*Q2/T0)))
         BB=(1._EB-EXP(-Q2*OMEGA/TEMP))*(1._EB-EXP(-Q2*OM3/TEMP))*(1._EB-EXP(-OM12*Q2/TEMP))**3*(1._EB+EXP(-OM12*Q2/TEMP)) &
           *(1._EB-EXP(-Q2*OMPRIM/TEMP))
         CC=AA*BB*OMEGA/TEMP*T0/TEMP
         L202: DO J=1,20
            V=FLOAT(J-1)
            IF(J/2*2==J)G=(V+1._EB)*(V+3._EB)/4._EB
            IF(J/2*2/=J)G=(V+2._EB)*(V+2._EB)/4._EB
            L201: DO K=1,10
               V3=FLOAT(K-1)
               QQ=(V3+1._EB)*G*EXP(-(V3*OM3+V*OM12)*Q2/TEMP)
               GAM=B-A*(V3+1._EB)
               OMVV3=OM3+.5_EB*X13+X23+2._EB*X33+XBAR*V+2._EB*X33*V3
               DELTA=A*(OMEGA-OMVV3)
               IF(GAM*GAM<=DELTA) CYCLE L202
               D=2._EB*(GAM*GAM-DELTA)**.5_EB
               OMVBAR=OMVV3*(1._EB-EXP(-OMVV3*Q2/TEMP))
               F1=GAM-D/2
               F2=GAM+D/2._EB
               EE=Q2*GAM/(A*A*TEMP)
               UNFLO1=EE*DELTA*(1._EB+.5_EB*A/GAM)
               IF(UNFLO1<=-78.) CYCLE L202
               UNFLO2=EE*2._EB*GAM*F1
               IF(UNFLO2>=78.) CYCLE L202
               FF=EXP(EE*DELTA*(1._EB+.5_EB*A/GAM))
               SMINUS=CC*QQ/OMVBAR*ABS(F1)*FF*EXP(-EE*2._EB*GAM*F1)
               UNFLO3=EE*2._EB*GAM*F2
               IF(UNFLO3>=78.) THEN
                  SPLUS=0.
               ELSE
                  SPLUS=CC*QQ/OMVBAR*ABS(F2)*FF*EXP(-EE*2._EB*GAM*F2)
               ENDIF
               GG=SDWEAK
               SDWEAK=(SMINUS+SPLUS)/D+SDWEAK
               TEST=(SDWEAK-GG)/SDWEAK
               IF(TEST<.0001) CYCLE L202
               SDSTRG=(.5_EB*G)**.5_EB*(SMINUS**.5+SPLUS**.5)/D+SDSTRG
            ENDDO L201
         ENDDO L202
         IF(SDWEAK==0._EB) THEN
            SDWEAK=0._EB
            GDINV=1._EB
            GDDINV=1._EB
         ELSE
            DINV=SDSTRG*SDSTRG/SDWEAK
            GDINV=GC1*DINV
            GDDINV=GD*DINV
!***  EXPRESS S/D AT STP, AS IS K IN NASA SP-3080
            SDWEAK=SDWEAK*TEMP/273.
         ENDIF
      ELSE
!CALCULATE ABSORPTION COEF. AND LINE SPACING PARAMETER FOR 2.7 MICRON BAND
         L=1
!CONTRIBUTION TO 2.7 MICRON BAND FROM (000)-(021) AND (010)-(031) TRANS.
         ALPHA=28.5_EB
         OMPRIM=2._EB*OM2+OM3
         L120: DO
         AA=ALPHA*B*Q2/(A*(1._EB-EXP(-OM3*Q2/T0))*(1._EB-EXP(-OM12*Q2 /T0))**3*(1._EB+EXP(-OM12*Q2/T0))*(1._EB-EXP(-OMPRIM*Q2/T0)))
         BB=(1._EB-EXP(-Q2*OMEGA/TEMP))*(1._EB-EXP(-Q2*OM3/TEMP))* (1._EB-EXP(-OM12*Q2/TEMP))**3*(1._EB+EXP(-OM12*Q2/TEMP)) &
           *(1._EB-EXP(-Q2*OMPRIM/TEMP))
         CC=AA*BB*OMEGA/TEMP*T0/TEMP
         L102: DO J=1,20
         V=FLOAT(J-1)
         IF(J/2*2==J)G=(V+1._EB)*(V+3._EB)/4._EB
         IF(J/2*2/=J)G=(V+2._EB)*(V+2._EB)/4._EB
         VBAR1=-1._EB+(V+3._EB)*(V+4._EB)/(V+2.)/6._EB
         IF(J/2*2==J)VBAR1=-1._EB+(V+5._EB)/6._EB
         L101: DO K=1,10
         V3=FLOAT(K-1)
         QQ=(V3+1)*G*EXP(-(V3*OM3+V*OM12)*Q2/TEMP)*(VBAR1+1._EB)
         GAM=B-A*(V3+1._EB)
         IF(L==2) THEN
            OMVV3=3728._EB-5._EB*V-47._EB*V3
            IF(V==0._EB)OMVV3=3715._EB-47._EB*V3
         ELSE
            OMVV3=3598._EB-18._EB*V-47._EB*V3
            IF(V==0._EB)OMVV3=3613._EB-47._EB*V3
         ENDIF
         DELTA=A*(OMEGA-OMVV3)
         IF(GAM*GAM<=DELTA) CYCLE L102
         D=2._EB*(GAM*GAM-DELTA)**.5_EB
         OMVBAR=OMVV3*(1._EB-EXP(-OMVV3*Q2/TEMP))
         F1=GAM-D/2._EB
         F2=GAM+D/2._EB
         EE=Q2*GAM/(A*A*TEMP)
         UNFLO1=EE*DELTA*(1._EB+.5_EB*A/GAM)
         IF(UNFLO1<=-78._EB) CYCLE L102
         UNFLO2=EE*2._EB*GAM*F1
         IF(UNFLO2>=78._EB) CYCLE L102
         FF=EXP(EE*DELTA*(1._EB+.5_EB*A/GAM))
         SMINUS=CC*QQ/OMVBAR*ABS(F1)*FF*EXP(-EE*2._EB*GAM*F1)
         UNFLO3=EE*2._EB*GAM*F2
         IF(UNFLO3>=78._EB) THEN
            SPLUS=0._EB
         ELSE
            SPLUS=CC*QQ/OMVBAR*ABS(F2)*FF*EXP(-EE*2._EB*GAM*F2)
         ENDIF
         GG=SDWEAK
         SDWEAK=(SMINUS+SPLUS)/D+SDWEAK
         TEST=(SDWEAK-GG)/SDWEAK
         IF(TEST<.0001_EB) CYCLE L102
         SDSTRG=(.5_EB*G)**.5_EB*(SMINUS**.5_EB+SPLUS**.5_EB)/D+SDSTRG
         ENDDO L101
         ENDDO L102
         IF(L==2) EXIT L120
!CONTRIBUTION TO 2.7 MICRON BAND FROM (000)-(101) AND (010)-(111) TRANS.
         ALPHA=42.3_EB
         OMPRIM=OM1+OM3
         L=2
      ENDDO L120
!CALCULATE ABSORPTION COEF AND LINE SPACING PARAMETER FOR 4.3 MICRON BAND
      IF(SDWEAK==0._EB) THEN
         SDWEAK=0._EB
         GDINV=1._EB
         GDDINV=1._EB
      ELSE
         DINV=SDSTRG*SDSTRG/SDWEAK
         GDINV=GC1*DINV
         GDDINV=GD*DINV
!***EXPRESS S/D AT STP, AS IS K IN NASA SP-3080
         SDWEAK=SDWEAK*TEMP/273.
      ENDIF
   ENDIF
ELSEIF((OMEGA<=3050._EB).AND.(OMEGA>2474.)) THEN
   SDWEAK=0._EB
   GDINV=1._EB
   GDDINV=1._EB
ELSEIF((OMEGA<=2474.).AND.(OMEGA>1975.)) THEN
   B=.391635_EB
   A=.0030875_EB
   X13=-19.37_EB
   X23=-12.53_EB
   X33=-12.63_EB
   OM1=1354.91_EB
   OM2=673._EB
   OM3=2396.49_EB
   T0=300._EB
   Q2=1.4388_EB
   XBAR=.5_EB*(.5_EB*X13+X23)
   OM12=.5_EB*(.5_EB*OM1+OM2)
   SDWEAK=0._EB
   SDSTRG=0._EB
   IF(OMEGA<=2395._EB) THEN
         ALPHA=2700._EB
         OMPRIM=OM3
         AA=ALPHA*B*Q2/(A*(1._EB-EXP(-OM3*Q2/T0))*(1._EB-EXP(-OM12*Q2 /T0))**3*(1._EB+EXP(-OM12*Q2/T0))*(1._EB-EXP(-OMPRIM*Q2/T0)))
         BB=(1._EB-EXP(-Q2*OMEGA/TEMP))*(1._EB-EXP(-Q2*OM3/TEMP))*(1._EB-EXP(-OM12*Q2/TEMP))**3*(1._EB+EXP(-OM12*Q2/TEMP)) &
           *(1._EB-EXP(-Q2*OMPRIM/TEMP))
         CC=AA*BB*OMEGA/TEMP*T0/TEMP
         L202A: DO J=1,20
            V=FLOAT(J-1)
            IF(J/2*2==J)G=(V+1._EB)*(V+3._EB)/4._EB
            IF(J/2*2/=J)G=(V+2._EB)*(V+2._EB)/4._EB
            L201A: DO K=1,10
               V3=FLOAT(K-1)
               QQ=(V3+1._EB)*G*EXP(-(V3*OM3+V*OM12)*Q2/TEMP)
               GAM=B-A*(V3+1._EB)
               OMVV3=OM3+.5_EB*X13+X23+2._EB*X33+XBAR*V+2._EB*X33*V3
               DELTA=A*(OMEGA-OMVV3)
               IF(GAM*GAM<=DELTA) CYCLE L202A
               D=2._EB*(GAM*GAM-DELTA)**.5_EB
               OMVBAR=OMVV3*(1._EB-EXP(-OMVV3*Q2/TEMP))
               F1=GAM-D/2._EB
               F2=GAM+D/2._EB
               EE=Q2*GAM/(A*A*TEMP)
               UNFLO1=EE*DELTA*(1._EB+.5_EB*A/GAM)
               IF(UNFLO1<=-78.) CYCLE L202A
               UNFLO2=EE*2._EB*GAM*F1
               IF(UNFLO2>=78.) CYCLE L202A
               FF=EXP(EE*DELTA*(1._EB+.5_EB*A/GAM))
               SMINUS=CC*QQ/OMVBAR*ABS(F1)*FF*EXP(-EE*2._EB*GAM*F1)
               UNFLO3=EE*2._EB*GAM*F2
               IF(UNFLO3>=78._EB) THEN
                  SPLUS=0._EB
                  ELSE
                  SPLUS=CC*QQ/OMVBAR*ABS(F2)*FF*EXP(-EE*2._EB*GAM*F2)
                  ENDIF
               GG=SDWEAK
               SDWEAK=(SMINUS+SPLUS)/D+SDWEAK
               TEST=(SDWEAK-GG)/SDWEAK
               IF(TEST<.0001_EB) CYCLE L202A
               SDSTRG=(.5_EB*G)**.5_EB*(SMINUS**.5+SPLUS**.5)/D+SDSTRG
            ENDDO L201A
         ENDDO L202A
         IF(SDWEAK==0._EB) THEN
            SDWEAK=0._EB
            GDINV=1._EB
            GDDINV=1._EB
         ELSE
            DINV=SDSTRG*SDSTRG/SDWEAK
            GDINV=GC1*DINV
            GDDINV=GD*DINV