radi.f90

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

                     CASE (MIRROR_BOUNDARY) 
                        WALL(IW)%ILW(N,IBND) = WALL(IW)%ILW(DLM(N,ABS(IOR)),IBND)
                        IL(II,JJ,KK) = WALL(IW)%ILW(N,IBND)
                     CASE (INTERPOLATED_BOUNDARY) 
                        IL(II,JJ,KK) = WALL(IW)%ILW(N,IBND)
                     CASE DEFAULT ! solid wall
                        WALL(IW)%ILW(N,IBND) = OUTRAD_W(IW) + RPI*(1._EB-E_WALL(IW))* INRAD_W(IW)
                  END SELECT
               ELSEIF (CYLINDRICAL) THEN
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) CYCLE WALL_LOOP1
                  IL(II,JJ,KK) = WALL(IW)%ILW(N,IBND)
               ENDIF
            ENDDO WALL_LOOP1
 
! Determine sweep direction in physical space
 
            ISTART = 1
            JSTART = 1
            KSTART = 1
            IEND   = IBAR
            JEND   = JBAR
            KEND   = KBAR
            ISTEP  = 1
            JSTEP  = 1
            KSTEP  = 1
            IF (DLX(N) < 0._EB) THEN
               ISTART = IBAR
               IEND   = 1
               ISTEP  = -1
            ENDIF
            IF (DLY(N) < 0._EB) THEN
               JSTART = JBAR
               JEND   = 1
               JSTEP  = -1
            ENDIF
            IF (DLZ(N) < 0._EB) THEN
               KSTART = KBAR
               KEND   = 1
               KSTEP  = -1
            ENDIF
 
            GEOMETRY: IF (CYLINDRICAL) THEN  ! Sweep in axisymmetric geometry

               J = 1
               CKLOOP: DO K=KSTART,KEND,KSTEP
                  CILOOP: DO I=ISTART,IEND,ISTEP
                     IC = CELL_INDEX(I,J,K)
                     IF (SOLID(IC)) CYCLE CILOOP
                     ILXU = IL(I-ISTEP,J,K)
                     ILYU = IL(I,J-JSTEP,K)
                     ILZU = IL(I,J,K-KSTEP)
                     IF (DLX(N)>=0._EB) THEN
                        RU  = R(I-1)
                        RD  = R(I)
                     ELSE
                        RU  = R(I)
                        RD  = R(I-1)
                     ENDIF
                     RP  = SQRT(0.5_EB*(RU**2+RD**2))
                     VC  = DX(I)  * RP*DPHI0 * DZ(K)
                     AXU =          RU       * DZ(K) * ABS(DLX(N))
                     AXD =          RD       * DZ(K) * ABS(DLX(N))
                     AYU = DX(I)             * DZ(K) * ABS(DLB(N))
                     AYD = DX(I)             * DZ(K) * ABS(DLY(N))
                     AZ  = DX(I)  * RP*DPHI0         * ABS(DLZ(N))
                     IF (MODULO(N,NRP(1))==1) AYD = 0._EB  ! Zero out the terms involving symmetric overhang
                     IF (MODULO(N,NRP(1))==0) AYU = 0._EB
                     IF (IC/=0) THEN
                        IW = WALL_INDEX(IC,-ISTEP)
                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                           ILXU = WALL(IW)%ILW(N,IBND)
                           AYU = 0.5*AYU
                           AZ = 0.5*AZ
                        ENDIF
                        IW = WALL_INDEX(IC,-JSTEP*2)
                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                           ILYU = WALL(IW)%ILW(N,IBND)
                           AXU = 0.5*AXU
                           AZ = 0.5*AZ
                        ENDIF
                        IW = WALL_INDEX(IC,-KSTEP*3)
                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                           ILZU = WALL(IW)%ILW(N,IBND)
                           AXU = 0.5*AXU
                           AYU = 0.5*AYU
                        ENDIF
                     ENDIF
                     RAP = 1._EB/(AXD+AYD+AZ+EXTCOE(I,J,K)*VC*RSA(N))
                     IL(I,J,K) = MAX(0._EB, RAP * (AXU*ILXU + AYU*ILYU + AZ*ILZU +  &
                                 VC*RSA(N)*RFPI*( KFST4(I,J,K)+KFST4W(I,J,K) +RSA_RAT*SCAEFF(I,J,K)*UIIOLD(I,J,K) ) ) )
                  ENDDO CILOOP
               ENDDO CKLOOP

            ELSEIF (TWO_D) THEN GEOMETRY  ! Sweep in 2D cartesian geometry

               J = 1
               K2LOOP: DO K=KSTART,KEND,KSTEP
                  I2LOOP: DO I=ISTART,IEND,ISTEP
                     IC = CELL_INDEX(I,J,K)
                     IF (SOLID(IC)) CYCLE I2LOOP
                     ILXU  = IL(I-ISTEP,J,K)
                     ILZU  = IL(I,J,K-KSTEP)
                     VC  = DX(I) * DZ(K)
                     AX  =         DZ(K) * ABS(DLX(N))
                     AZ  = DX(I)         * ABS(DLZ(N))
                     IF (IC/=0) THEN
                        IW = WALL_INDEX(IC,-ISTEP)
                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                           ILXU = WALL(IW)%ILW(N,IBND)
                           AZ = 0.5*AZ
                        ENDIF
                        IW = WALL_INDEX(IC,-KSTEP*3)
                        IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                           ILZU = WALL(IW)%ILW(N,IBND)
                           AX = 0.5*AX
                        ENDIF
                     ENDIF
                     RAP = 1._EB/(AX+AZ+EXTCOE(I,J,K)*VC*RSA(N))
                     IL(I,J,K) = MAX(0._EB, RAP * (AX*ILXU + AZ*ILZU + &
                                 VC*RSA(N)*RFPI*(KFST4(I,J,K)+KFST4W(I,J,K) +  RSA_RAT*SCAEFF(I,J,K)*UIIOLD(I,J,K) ) ) ) 
                  ENDDO I2LOOP
               ENDDO K2LOOP

            ELSE GEOMETRY  ! Sweep in 3D cartesian geometry

               KLOOP: DO K=KSTART,KEND,KSTEP
                  JLOOP: DO J=JSTART,JEND,JSTEP
                     ILOOP: DO I=ISTART,IEND,ISTEP
                        IC = CELL_INDEX(I,J,K)
                        IF (SOLID(IC)) CYCLE ILOOP
                        ILXU  = IL(I-ISTEP,J,K)
                        ILYU  = IL(I,J-JSTEP,K)
                        ILZU  = IL(I,J,K-KSTEP)
                        VC  = DX(I) * DY(J) * DZ(K)
                        AX  =         DY(J) * DZ(K) * ABS(DLX(N))
                        AY  = DX(I)         * DZ(K) * ABS(DLY(N))
                        AZ  = DX(I) * DY(J)         * ABS(DLZ(N))
                        IF (IC/=0) THEN
                           IW = WALL_INDEX(IC,-ISTEP)
                           IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                              AY = 0.5*AY
                              AZ = 0.5*AZ
                              ILXU = WALL(IW)%ILW(N,IBND)
                           ENDIF
                           IW = WALL_INDEX(IC,-JSTEP*2)
                           IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                              AX = 0.5*AX
                              AZ = 0.5*AZ
                              ILYU = WALL(IW)%ILW(N,IBND)
                           ENDIF
                           IW = WALL_INDEX(IC,-KSTEP*3)
                           IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
                              AX = 0.5*AX
                              AY = 0.5*AY
                              ILZU = WALL(IW)%ILW(N,IBND)
                           ENDIF
                        ENDIF
                        RAP = 1._EB/(AX+AY+AZ+EXTCOE(I,J,K)*VC*RSA(N))
                        IL(I,J,K) = MAX(0._EB, RAP * ( AX*ILXU + AY*ILYU + AZ*ILZU + &
                                    VC*RSA(N)*RFPI*( KFST4(I,J,K)+KFST4W(I,J,K) + RSA_RAT*SCAEFF(I,J,K)*UIIOLD(I,J,K) ) ) )
                     ENDDO ILOOP
                  ENDDO JLOOP
               ENDDO KLOOP
 
            ENDIF GEOMETRY
 
! Boundary values: Incoming radiation
 
            WALL_LOOP2: DO IW=1,NWC
               IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY)   CYCLE WALL_LOOP2     
               IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY)   CYCLE WALL_LOOP2  
               IF (BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_LOOP2  
               IOR = IJKW(4,IW)
               IF (TWO_D .AND. .NOT.CYLINDRICAL  .AND. ABS(IOR)==2) CYCLE WALL_LOOP2  ! 2-D non cylindrical
               IF (DLN(IOR,N)>=0._EB) CYCLE WALL_LOOP2     ! outgoing
               WAXIDLN = - W_AXI*DLN(IOR,N)
               IIG = IJKW(6,IW)
               JJG = IJKW(7,IW)
               KKG = IJKW(8,IW)
               INRAD_W(IW) = INRAD_W(IW) - WAXIDLN * WALL(IW)%ILW(N,IBND) ! update incoming radiation,step 1
               WALL(IW)%ILW(N,IBND) = IL(IIG,JJG,KKG)
               INRAD_W(IW) = INRAD_W(IW) + WAXIDLN * WALL(IW)%ILW(N,IBND) ! update incoming radiation,step 2
            ENDDO WALL_LOOP2
 
! Copy the Y-downwind intensities to Y-upwind in cylindrical case
 
            IF (CYLINDRICAL) THEN
               J=1
               DO K=1,KBAR
                  DO I=1,IBAR
                     IWUP   = WALL_INDEX(CELL_INDEX(I,J,K),-2)
                     IWDOWN = WALL_INDEX(CELL_INDEX(I,J,K), 2)
                     WALL(IWUP)%ILW(MAX(1,N-1),IBND) = WALL(IWDOWN)%ILW(N,IBND)
                  ENDDO
               ENDDO
            ENDIF
 
! Calculate integrated intensity UIID
 
            IF (WIDE_BAND_MODEL) THEN
               UIID(:,:,:,IBND) = UIID(:,:,:,IBND) + W_AXI*RSA(N)*IL
            ELSE
               UIID(:,:,:,ANGLE_INC_COUNTER) = UIID(:,:,:,ANGLE_INC_COUNTER) + W_AXI*RSA(N)*IL
            ENDIF
 
! Interpolate boundary intensities onto other meshes
 
            INTERPOLATION_LOOP: DO NOM=1,NMESHES
               IF (NM==NOM) CYCLE INTERPOLATION_LOOP
               IF (EVACUATION_ONLY(NOM)) CYCLE INTERPOLATION_LOOP
               IF (NIC(NOM,NM)==0) CYCLE INTERPOLATION_LOOP
               M2=>OMESH(NOM)
               OTHER_WALL_LOOP: DO IW=1,MESHES(NOM)%NEWC
                  IF (M2%IJKW(9,IW)/=NM .OR. M2%BOUNDARY_TYPE(IW)/=INTERPOLATED_BOUNDARY) CYCLE OTHER_WALL_LOOP
                  IOR = M2%IJKW(4,IW)
                  IF (DLN(IOR,N)<=0._EB) CYCLE OTHER_WALL_LOOP
                  M2%WALL(IW)%ILW(N,IBND)=IL(M2%IJKW(10,IW),M2%IJKW(11,IW),M2%IJKW(12,IW))
               ENDDO OTHER_WALL_LOOP
            ENDDO INTERPOLATION_LOOP
 
         ENDDO ANGLE_LOOP
      ENDDO UPDATE_LOOP
 
      ! Compute incoming flux on walls 
      DO IW=1,NWC
         IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) CYCLE 
         IBC = IJKW(5,IW)
         QRADIN(IW)  = QRADIN(IW) + E_WALL(IW)*(INRAD_W(IW)+SURFACE(IBC)%EXTERNAL_FLUX/NUMBER_SPECTRAL_BANDS)
      ENDDO 
   ENDIF INTENSITY_UPDATE
 
   ! Save source term for the energy equation (QR = -DIV Q)
   IF (WIDE_BAND_MODEL) THEN
      QR = QR + KAPPA*UIID(:,:,:,IBND)-KFST4
      IF (NLP>0 .AND. N_EVAP_INDICIES>0) QR_W = QR_W + KAPPAW*UIID(:,:,:,IBND) - KFST4W
   ENDIF

ENDDO BAND_LOOP

IF (UPDATE_INTENSITY) THEN
   ! Sum up the parts of the intensity
   UII = SUM(UIID, DIM = 4)
ENDIF

! Save source term for the energy equation (QR = -DIV Q). Done only in one-band (gray gas) case.
IF (.NOT. WIDE_BAND_MODEL) THEN
   QR  = KAPPA*UII - KFST4
   IF (NLP>0 .AND. N_EVAP_INDICIES>0) QR_W = QR_W + KAPPAW*UII - KFST4W
ENDIF
END SUBROUTINE RADIATION_FVM


END SUBROUTINE COMPUTE_RADIATION

REAL(EB) FUNCTION ZZ2KAPPA(ZZ,ZZF,NIYY,TTD,BAND,NS)
 
! Calculate KAPPA as a function of ZZ.
 
REAL(EB) :: ZZ,ZZF,IYY
INTEGER  :: TTD,BAND,IY1,IY2,ZINT,NIYY,NS
REAL(EB), POINTER, DIMENSION(:,:,:) :: KAPPA
KAPPA => SPECIES(NS)%KAPPA
 
IF (ZZ<=ZZF .OR. ZZF==1._EB) THEN
   IYY = ZZ*REAL(NIYY/2,EB)/ZZF
ELSE
   IYY = NIYY/2 + REAL(NIYY-NIYY/2,EB)*(ZZ-ZZF)/(1._EB-ZZF) 
ENDIF
ZINT   = MAX(0,MIN(NIYY,NINT(IYY)))
IF(IYY==ZINT) THEN
   ZZ2KAPPA=KAPPA(ZINT,TTD,BAND)
ELSE
   ZINT = MAX(0,MIN(NIYY,NINT(IYY)))
   IY2  = MIN(NIYY,ZINT + 1)
   IY1  = ZINT
   ZZ2KAPPA=KAPPA(IY1,TTD,BAND)+(IYY-REAL(IY1,EB))* (KAPPA(IY2,TTD,BAND)-KAPPA(IY1,TTD,BAND))
ENDIF
END FUNCTION ZZ2KAPPA
 
 
REAL(EB) FUNCTION IZ2ZZ(IZ,ZZF,NIZ)
 
! Function gives the mixture fraction as a function of index IZ
 
REAL(EB) :: ZZF
INTEGER  :: IZ, NIZ
 
IF (IZ<=NIZ/2) THEN
   IZ2ZZ = REAL(IZ,EB)*ZZF/REAL(NIZ/2,EB)
ELSE
   IZ2ZZ = ZZF+(1._EB-ZZF)*REAL(IZ-NIZ/2,EB)/REAL(NIZ-NIZ/2,EB)
ENDIF
 
END FUNCTION IZ2ZZ
 

REAL(EB) FUNCTION BLACKBODY_FRACTION(L1,L2,TEMP)
 
! Calculates the fraction of black body radiation between wavelengths L1 and L2 (micron) in Temperature TEMP
 
REAL(EB) :: L1,L2,TEMP,LT1,LT2,BBFLOW,BBFHIGH
INTEGER  :: IYY

LT1    =   L1 * TEMP/LTSTEP
LT2    =   L2 * TEMP/LTSTEP
IYY = MIN(NLAMBDAT,MAX(0,NINT(LT1)))
BBFLOW = BBFRAC(IYY)
IYY = MIN(NLAMBDAT,MAX(0,NINT(LT2)))
BBFHIGH = BBFRAC(IYY)
BLACKBODY_FRACTION = BBFHIGH - BBFLOW

END FUNCTION BLACKBODY_FRACTION

SUBROUTINE GET_KAPPA(Z1,Z2,Z3,YY_SUM,KAPPA,TYY,IBND)
USE PHYSICAL_FUNCTIONS, ONLY : GET_F,GET_F_C
! GET_KAPPA returns the radiative absorption
REAL(EB) :: Z_1,Z_2,Z_3,YY_SUM,F,C,ZZ,Z_F,KAPPA,KAPPA_N(N_SPECIES),Z1,Z2,Z3
INTEGER :: IBND,TYY

IF (YY_SUM >=1._EB) THEN
   KAPPA = 0._EB
   RETURN
ELSE
   YY_SUM = MAX(0._EB,YY_SUM)
   Z_1 = MAX(0._EB,Z1)/(1._EB - YY_SUM)
   Z_2 = MAX(0._EB,Z2)/(1._EB - YY_SUM)         
   Z_3 = MAX(0._EB,Z3)/(1._EB - YY_SUM)           
ENDIF
selectmethod: IF (.NOT. CO_PRODUCTION) THEN
   Z_F = REACTION(1)%Z_F
   CALL GET_F(Z_1,Z_3,F,Z_F)
   ZZ = MIN(1._EB,Z_1 + Z_3)
   KAPPA_N(I_FUEL)   = ZZ2KAPPA(ZZ,REACTION(1)%Z_F,SPECIES(I_FUEL)%NKAP_MASS,TYY,IBND,I_FUEL)
   KAPPA_N(I_PROG_F) = ZZ2KAPPA(ZZ,REACTION(2)%Z_F,SPECIES(I_PROG_F)%NKAP_MASS,TYY,IBND,I_PROG_F)
   KAPPA = (1._EB - F) * KAPPA_N(I_FUEL) + F * KAPPA_N(I_PROG_F)    
ELSE selectmethod
   CALL GET_F_C(Z_1,Z_2,Z_3,F,C,Z_F)        
   ZZ = MIN(1._EB,Z_1 + Z_2 + Z_3)
   KAPPA_N(I_FUEL)    = ZZ2KAPPA(ZZ,REACTION(1)%Z_F,SPECIES(I_FUEL)%NKAP_MASS,TYY,IBND,I_FUEL)
   KAPPA_N(I_PROG_CO) = ZZ2KAPPA(ZZ,REACTION(2)%Z_F,SPECIES(I_PROG_CO)%NKAP_MASS,TYY,IBND,I_PROG_CO)   
   KAPPA_N(I_PROG_F) = ZZ2KAPPA(ZZ,REACTION(3)%Z_F,SPECIES(I_PROG_F)%NKAP_MASS,TYY,IBND,I_PROG_F)   
   KAPPA = (1._EB - F) * ( (1._EB - C) * KAPPA_N(I_FUEL) + C * KAPPA_N(I_PROG_CO)) + F *KAPPA_N(I_PROG_F) 
ENDIF selectmethod
KAPPA = KAPPA/(1._EB - YY_SUM)

END SUBROUTINE GET_KAPPA 

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

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

END SUBROUTINE GET_REV_radi
 
END MODULE RAD