irad.f90

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

!***  EXPRESS S/D AT STP, AS IS K IN NASA SP-3080
            SDWEAK=SDWEAK*TEMP/273._EB
         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
         L120A: 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
         L102A: 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._EB)/6._EB
            IF(J/2*2==J)VBAR1=-1._EB+(V+5._EB)/6._EB
            L101A: 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 L102A
               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 L102A
               UNFLO2=EE*2._EB*GAM*F1
               IF(UNFLO2>=78._EB) CYCLE L102A
               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) CYCLE L102A
               SDSTRG=(.5_EB*G)**.5_EB*(SMINUS**.5+SPLUS**.5)/D+SDSTRG
            ENDDO L101A
         ENDDO L102A
         IF(L==2) EXIT L120A
!CONTRIBUTION TO 2.7 MICRON BAND FROM (000)-(101) AND (010)-(111) TRANS.
         ALPHA=42.3_EB
         OMPRIM=OM1+OM3
         L=2
      ENDDO L120A
!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._EB
         ENDIF
      ENDIF
   ELSEIF((OMEGA<=1975._EB).AND.(OMEGA>1100._EB)) THEN
      SDWEAK=0._EB
      GDINV=1._EB
      GDDINV=1._EB
   ELSEIF((OMEGA<=1100._EB).AND.(OMEGA>880._EB)) THEN
!CONTRIBUTION TO 10.0 MICRON BAND FROM (100)-(001) AND (020)-(001) TRANS.
      OM1=1354.91_EB
      OM2=673._EB
      OM3=2396.49_EB
      Q2=1.4388_EB
      BCNT(1)=960.8_EB
      BCNT(2)=1063.6_EB
      OMA=OM3
      OMB=(OM1+2._EB*OM2)/2._EB
      T0=300._EB
      ATOT(1)=0.0219_EB
      ATOT(2)=0.0532_EB
      BE=0.391635_EB
      DO K=1,2
         ATOT(K)=T0/TEMP*ATOT(K)*EXP(Q2*OMB*(1._EB/T0-1._EB/TEMP)) *(1._EB-EXP(-Q2*(OMA-OMB)/TEMP))/(1._EB-EXP(-Q2*OMA/TEMP)) &
                   /(1._EB-EXP(-OMB*Q2/TEMP))
      ENDDO
      SDWEAK=0._EB
      DO I=1,2
         SDWEAK=SDWEAK+ATOT(I)*Q2/(4._EB*BE*TEMP)*ABS(OMEGA-BCNT(I))  *EXP(-Q2/(4._EB*BE*TEMP)*(OMEGA-BCNT(I))**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.
   ELSEIF((OMEGA<=880._EB).AND.(OMEGA>500._EB))  THEN
!CONTRIBUTION TO 15.0 MICRON BAND FROM (000)-(010) TRANS.
      TTEMP=TEMP
      J=(OMEGA-495._EB)/5._EB
      W1=495._EB+5._EB*FLOAT(J)
      WW=(OMEGA-W1)/5
      IF(TEMP>=2400._EB)TEMP=2399.99_EB
      IF(TEMP<300._EB)TEMP=300._EB
      I=TEMP/300._EB
      IF((I>2).AND.(TEMP<1200._EB)) THEN
         I=2
         TT=(TEMP-600._EB)/600._EB
      ELSEIF((I>5).AND.(TEMP<2400._EB)) THEN
         I=5
         TT=(TEMP-1800._EB)/600._EB
      ELSE
         T1=FLOAT(I)*300._EB
         TT=(TEMP-T1)/300._EB
         IF(I>4)I=I-1
         ENDIF
      TW=TT*WW
      SDWEAK=SD15(I,J)*(1._EB-TT-WW+TW)+SD15(I+1,J)*(TT-TW)  +SD15(I,J+1)*(WW-TW)+SD15(I+1,J+1)*TW
      IF(SDWEAK==0._EB) THEN
         SDWEAK=0._EB
         GDINV=1._EB
         GDDINV=1._EB
         ELSE
!CALCULATE LINE SPACING PARAMETER FOR 15.0 MICRON BAND
         DINV1=1.2_EB
         DINV2=8.0_EB
         DINV3=30.0_EB
         TEMP1=300.0_EB
         TEMP2=550.0_EB
         TEMP3=830.0_EB
         DINV=DINV1*(TEMP-TEMP2)*(TEMP-TEMP3)/(TEMP1-TEMP2)  /(TEMP1-TEMP3)+DINV2*(TEMP-TEMP1)*(TEMP-TEMP3)&
                  /(TEMP2-TEMP1)/(TEMP2-TEMP3)+DINV3*(TEMP-TEMP1)  *(TEMP-TEMP2)/(TEMP3-TEMP1)/(TEMP3-TEMP2)
         GDINV=GC1*DINV
         GDDINV=GD*DINV
         ENDIF
         TEMP = TTEMP  ! Line added by Jason Floyd, Aug 30, 2002
   ELSE
      SDWEAK=0._EB
      GDINV=1._EB
      GDDINV=1._EB
      ENDIF
   ENDIF
END SUBROUTINE CO2

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

SUBROUTINE H2O(OMEGA,TEMP,GC2,SDWEAK,GDINV,GDDINV)

INTEGER I,J
REAL(EB) OMEGA,TEMP,GC2,SDWEAK,GDINV,GDDINV,WM,W1,WW,T1,TT,TW, D,B,DINV,TTEMP,GD

IF (OMEGA>=9300..OR.OMEGA<50._EB) THEN
   SDWEAK=0._EB
   GDINV=1._EB
   GDDINV=1._EB
   ELSE
   WM=18._EB
   GD=5.94E-6_EB*OMEGA*(TEMP/(273._EB*WM))**.5_EB
   J=(OMEGA-25._EB)/25._EB
   TTEMP=TEMP
   IF(TEMP>=2500._EB) TEMP=2499.99_EB
   IF(TEMP<300._EB) TEMP=300._EB
   I=TEMP/500._EB +1
   IF(I==2.AND.TEMP<600._EB) I=1
   W1=25._EB+25._EB*FLOAT(J)
   WW=(OMEGA-W1)/25._EB
   IF(I>2) THEN
      T1=FLOAT(I-1)*500._EB
      TT=(TEMP-T1)/500._EB
      ELSE
      IF(I==1) TT=(TEMP-300._EB)/300._EB
      IF(I==2) TT=(TEMP-600._EB)/400._EB
      ENDIF
   TW=TT*WW
   SDWEAK=SD(I,J)*(1._EB-TT-WW+TW)+SD(I+1,J)*(TT-TW)+SD(I,J+1) *(WW-TW)+SD(I+1,J+1)*TW
   D=-2.294_EB+.3004E-02_EB*TEMP-.366E-06_EB*TEMP**2
   B=SIN(.0036_EB*OMEGA-8.043_EB)
   DINV=EXP(.7941_EB*B+D)
!     DINV=EXP(0.00106*TEMP-1.21)
   GDINV=GC2*DINV
   GDDINV=GD*DINV
   TEMP=TTEMP
ENDIF
END SUBROUTINE H2O

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

SUBROUTINE CO(OMEGA,TEMP,GC4,SDWEAK,GDINV,GDDINV)

INTEGER J
REAL(EB) OMEGA,TEMP,GC4,SDWEAK,GDINV,GDDINV,AA,BB,CC,QQ,EE,FF,GG, &
         SMINUS,SPLUS,SDSTRG,B,ALPHA,A,OME,WX,WY,OMPRIM,T0, &
         Q2,WM,GD,V,GAM,OMV,DELTA,D,OMVBAR,F1,F2,TEST,DINV

IF(OMEGA<1600._EB .OR. OMEGA>2400._EB) THEN
   SDWEAK=0._EB
   GDINV=1._EB
   GDDINV=1._EB
ELSE
   B=1.93139_EB
   ALPHA=260._EB
   A=.017485_EB
   OME=2170.21_EB
   WX=13.461_EB
   WY=.0308_EB
   OMPRIM=OME-2._EB*WX+3.25_EB*WY
   T0=300._EB
   Q2=1.4388_EB
   WM=28._EB
   GD=5.94E-6_EB*OMEGA*(TEMP/(273._EB*WM))**.5_EB
   SDWEAK=1.E-99_EB
   SDSTRG=1.E-99_EB
   AA=ALPHA*B*Q2/(A*(1._EB-EXP(-OMPRIM*Q2/T0))**2)
   BB=(1._EB-EXP(-OMEGA*Q2/TEMP))*(1._EB-EXP(-OMPRIM*Q2/TEMP))**2
   CC=AA*BB*OMEGA/TEMP*T0/TEMP
   L101: DO J=1,20
   V=FLOAT(J-1)
   QQ=(V+1._EB)*EXP(-V*OME*Q2/TEMP)
   GAM=B-A*(V+1._EB)
   OMV=OME-2._EB*(V+1._EB)*WX+(3._EB*(V+1._EB)*(V+1._EB)+.25_EB)*WY
   DELTA=A*(OMEGA-OMV)
   IF(GAM*GAM<=DELTA) EXIT L101
   D=2._EB*(GAM*GAM-DELTA)**.5_EB
   OMVBAR=OMV*(1._EB-EXP(-OMV*Q2/TEMP))
   F1=GAM-D/2._EB
   F2=GAM+D/2._EB
   EE=Q2*GAM/(A*A*TEMP)
   FF=EXP(EE*DELTA*(1._EB+.5_EB*A/GAM))
   SMINUS=CC*QQ/OMVBAR*ABS(F1)*FF*EXP(-EE*2._EB*GAM*F1)
   SPLUS=CC*QQ/OMVBAR*ABS(F2)*FF*EXP(-EE*2._EB*GAM*F2)
   GG=SDWEAK
   SDWEAK=(SMINUS+SPLUS)/D+SDWEAK
   TEST=(SDWEAK-GG)/SDWEAK
   IF(TEST<.0001_EB) EXIT L101
   SDSTRG=(SMINUS**.5_EB+SPLUS**.5_EB)/D+SDSTRG
   ENDDO L101
   DINV=SDSTRG*SDSTRG/SDWEAK
   GDINV=GC4*DINV
   GDDINV=GD*DINV
!***EXPRESS S/D AT STP, AS IS K IN NASA SP-3080
   SDWEAK=SDWEAK*TEMP/273._EB
ENDIF
END SUBROUTINE CO

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

      SUBROUTINE POD(OMEGA)

!***POD CALCULATES PARTICLE OPTICAL DEPTH, XPART, OF THE VOLUME
!   FRACTION OF SOOT PARTICLES IN GAS CLOUD.  RIN AND RIK ARE
!   THE REAL AND IMAGINARY PARTS OF THE INDEX OF REFRACTION.  THE
!   PARTICLES ARE ASSUMED TO BE IN THE RAYLEIGH LIMIT.

      INTEGER J
      REAL(EB) OMEGA,ABCO,FF,LAMBDA!,RIN,RIK

      LAMBDA=10000._EB/OMEGA
!     RIK=.5
!     RIN=1.6
!      FF=36.*3.1416*RIN*RIK/LAMBDA/((RIN*RIN-RIK*RIK+2.)**2 +
!     .   (2._EB*RIN*RIK)**2) ! OLD DALZELL AND SAOFIM MODEL
!     FF=7./LAMBDA
!
!     ABSORPTION COEF. IS BASED UPON MEASUREMENTS OF WIDMANN AND
!     MULHOLLAND
      FF=8.9/LAMBDA
      DO J=1,NPT
         ABCO=FF*SVF(J)*1.E6_EB
         IF(J==1) THEN
            XPART(1)=ABCO*DD(1)
         ELSE
            XPART(J)=XPART(J-1)+ABCO*DD(J)
         ENDIF
      ENDDO
      END SUBROUTINE POD

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

SUBROUTINE FUEL(OMEGA,TEMP,PCH4,PTOT,GC3,SDWEAK,GDINV,GDDINV)

INTEGER I,J
REAL(EB) OMEGA,TEMP,PCH4,PTOT,GC3,SDWEAK,GDINV,GDDINV,BE,Q2, &
      WM,GD,OM1,OM2,OM3,OM4,COM1,COM2,COM3,COM4,DINV,PE,W1,SDB, &
      SDA,SDC