read.f90

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

 
 
SUBROUTINE READ_TIME
 
REAL(EB) :: DT,VEL_CHAR,TWFIN
INTEGER :: NM
NAMELIST /TIME/ DT,T_BEGIN,T_END,TWFIN,FYI,WALL_INCREMENT,SYNCHRONIZE, &
                EVAC_DT_FLOWFIELD,EVAC_DT_STEADY_STATE,TIME_SHRINK_FACTOR
TYPE (MESH_TYPE), POINTER :: M
 
DT                   = -1._EB
EVAC_DT_FLOWFIELD    = 0.01_EB
EVAC_DT_STEADY_STATE = 0.05_EB
SYNCHRONIZE          = .FALSE.
TIME_SHRINK_FACTOR   = 1._EB
T_BEGIN              = 0._EB
T_END                = -9999999_EB
TWFIN                = -99999999._EB
WALL_INCREMENT       = 2
 
REWIND(LU_INPUT)
READ_TIME_LOOP: DO
   CALL CHECKREAD('TIME',LU_INPUT,IOS)
   IF (IOS==1) EXIT READ_TIME_LOOP
   READ(LU_INPUT,TIME,END=21,ERR=22,IOSTAT=IOS)
   22 IF (IOS>0) CALL SHUTDOWN('ERROR: Problem with TIME line')
   IF (TWFIN > T_END) T_END=TWFIN
ENDDO READ_TIME_LOOP
21 REWIND(LU_INPUT)
 
IF (T_END<=T_BEGIN) SET_UP = .TRUE.
T_END = T_BEGIN + (T_END-T_BEGIN)/TIME_SHRINK_FACTOR
 
IF (SYNCHRONIZE) SYNC_TIME_STEP = .TRUE.
IF (ANY(SYNC_TIME_STEP)) SYNCHRONIZE = .TRUE.
 
MESH_LOOP: DO NM=1,NMESHES
   M=>MESHES(NM)
   IF (DT>0._EB) THEN
      M%DT = DT
   ELSE
      VEL_CHAR = 0.2_EB*SQRT(10._EB*(M%ZF-M%ZS))
      IF (TWO_D) THEN
         M%DT=(M%DXMIN*M%DZMIN)**(1._EB/2._EB)/VEL_CHAR
      ELSE
         M%DT=(M%DXMIN*M%DYMIN*M%DZMIN)**(1._EB/3._EB)/VEL_CHAR
      ENDIF
   ENDIF
   IF (EVACUATION_ONLY(NM)) THEN
      SYNC_TIME_STEP(NM) = .FALSE.
      M%DT = EVAC_DT_FLOWFIELD
   ENDIF
ENDDO MESH_LOOP
 
END SUBROUTINE READ_TIME
 
 
SUBROUTINE READ_MISC
USE MATH_FUNCTIONS, ONLY: GET_RAMP_INDEX
 
REAL(EB) :: X_H2O_TMPA,X_H2O_40_C,C_HORIZONTAL,C_VERTICAL,MW,VISCOSITY,CONDUCTIVITY
CHARACTER(30) :: RAMP_GX,RAMP_GY,RAMP_GZ
NAMELIST /MISC/ PR,SC,TMPA,GVEC,RELAXATION_FACTOR,FYI, &
                CSMAG,RAMP_GX,RAMP_GY,RAMP_GZ,BAROCLINIC, &
                LAPSE_RATE,ISOTHERMAL, &
                P_INF,SURF_DEFAULT,EVAC_SURF_DEFAULT, &
                C_FORCED,C_VERTICAL,C_HORIZONTAL,H_FIXED,RESTART,ASSUMED_GAS_TEMPERATURE, &
                BACKGROUND_SPECIES,MW,LES,DNS, &
                VISCOSITY,CONDUCTIVITY,NOISE, &
                RADIATION,GAMMA,BNDF_DEFAULT, &
                U0,V0,W0,HUMIDITY, &
                POROUS_FLOOR,SUPPRESSION,CO_PRODUCTION, &
                TEXTURE_ORIGIN,NSTRATA, &
                THICKEN_OBSTRUCTIONS, &
                EVAC_PRESSURE_ITERATIONS,EVAC_TIME_ITERATIONS, &
                PRESSURE_CORRECTION,CHECK_POISSON,STRATIFICATION,RESTART_CHID, &
                CFL_MAX,CFL_MIN,VN_MAX,VN_MIN,SOLID_PHASE_ONLY
 
! Physical constants
 
R0      = 8314.472_EB     ! Universal Gas Constant (J/kmol/K) (NIST Physics Constants)
R1      = 1.986257E-03_EB ! Universal Gas Constant (kcal/mol/K)
TMPA    = 20._EB        ! Ambient temperature (C)
GRAV    = 9.80665_EB    ! Acceleration of gravity (m/s**2)
GAMMA   = 1.4_EB        ! Heat capacity ratio for air
P_INF   = 101325._EB    ! Ambient pressure (Pa)
P_STP   = 101325._EB    ! Standard pressure (Pa)
TMPM    = 273.15_EB     ! Melting temperature of water (K)
SIGMA   = 5.67E-8_EB    ! Stefan-Boltzmann constant (W/m**2/K**4)
C_P_W   = 4184._EB           ! Specific Heat of Water (J/kg/K)
H_V_W   = 2259._EB*1000._EB  ! Heat of Vap of Water (J/kg)
MW_AIR  = 28.8_EB        ! g/mol
MW_SOOT = 0.9_EB * 12._EB + 0.1_EB * 1._EB
HUMIDITY= -1._EB                           ! Relative Humidity
RHO_SOOT= 1850._EB                         ! Density of soot particle (kg/m3)
RESTART_CHID = CHID
 
! Empirical heat transfer constants
 
C_VERTICAL   = 1.31_EB  ! Vertical free convection (Holman, Table 7-2)
C_HORIZONTAL = 1.52_EB  ! Horizontal free convection 
C_FORCED     = 0.037_EB ! Forced convection coefficient
H_FIXED                 = -1.      ! Fixed heat transfer coefficient, used for diagnostics
ASSUMED_GAS_TEMPERATURE = -1000.   ! Assumed gas temperature, used for diagnostics
 
! Often used numbers
 
PI      = 4._EB*ATAN(1.0_EB)
SQRTPI  = SQRT(PI)
RPI     = 1._EB/PI
TWOPI   = 2._EB*PI
PIO2    = PI/2._EB
ONTH    = 1._EB/3._EB
THFO    = 3._EB/4._EB
ONSI    = 1._EB/6._EB
TWTH    = 2._EB/3._EB
FOTH    = 4._EB/3._EB
RFPI    = 1._EB/(4._EB*PI)
 
! Background parameters
 
U0 = 0._EB
V0 = 0._EB
W0 = 0._EB
BACKGROUND_SPECIES = 'AIR'
VISCOSITY = -1._EB
CONDUCTIVITY = -1._EB
MU_USER = -1._EB
K_USER  = -1._EB
MW      = 0._EB   

! Logical constants

RESTART        = .FALSE.
RADIATION      = .TRUE.
SUPPRESSION    = .TRUE.
CO_PRODUCTION  = .FALSE.
CHECK_POISSON  = .FALSE.
BAROCLINIC     = .FALSE.
NOISE          = .TRUE.
ISOTHERMAL     = .FALSE.  
BNDF_DEFAULT   = .TRUE.
LES            = .TRUE.
DNS            = .FALSE.
POROUS_FLOOR   = .TRUE.
PRESSURE_CORRECTION  = .FALSE.
STRATIFICATION = .TRUE.
SOLID_PHASE_ONLY = .FALSE.

TEXTURE_ORIGIN(1) = 0._EB
TEXTURE_ORIGIN(2) = 0._EB
TEXTURE_ORIGIN(3) = 0._EB
 
! EVACuation parameters
 
EVAC_PRESSURE_ITERATIONS = 50
EVAC_TIME_ITERATIONS     = 50
 
! LES parameters
 
CSMAG                = 0.20_EB  ! Smagorinsky constant
PR                   = -1.0_EB  ! Turbulent Prandtl number
SC                   = -1.0_EB  ! Turbulent Schmidt number
 
! Misc
 
RAMP_GX              = 'null'
RAMP_GY              = 'null'
RAMP_GZ              = 'null'
SURF_DEFAULT         = 'INERT'
EVAC_SURF_DEFAULT    = 'INERT'
GVEC(1)              = 0._EB        ! x-component of gravity 
GVEC(2)              = 0._EB        ! y-component of gravity 
GVEC(3)              = -GRAV        ! z-component of gravity 
LAPSE_RATE           = 0._EB       
RELAXATION_FACTOR    = 1.00_EB      ! Relaxation factor for no-flux
NSTRATA              = 7            ! Number bins for drop dist.
THICKEN_OBSTRUCTIONS = .FALSE.
CFL_MAX              = 1.0_EB       ! Stability bounds
CFL_MIN              = 0.8_EB
VN_MAX               = 1.0_EB
VN_MIN               = 0.8_EB
 
REWIND(LU_INPUT)
MISC_LOOP: DO 
   CALL CHECKREAD('MISC',LU_INPUT,IOS)
   IF (IOS==1) EXIT MISC_LOOP
   READ(LU_INPUT,MISC,END=23,ERR=24,IOSTAT=IOS)
   24 IF (IOS>0) CALL SHUTDOWN('ERROR: Problem with MISC line')
ENDDO MISC_LOOP
23 REWIND(LU_INPUT)
 
! Temperature conversions

H0    = 0.5_EB*(U0**2+V0**2+W0**2)
TMPA  = TMPA + TMPM
TMPA4 = TMPA**4
 
! Humidity (40% by default, but limited for high temps)
 
IF (HUMIDITY < 0._EB) THEN
   X_H2O_TMPA = MIN( 1._EB , EXP(-(H_V_W*MW_H2O/R0)*(1._EB/TMPA     -1._EB/373.15_EB)) )
   X_H2O_40_C =              EXP(-(H_V_W*MW_H2O/R0)*(1._EB/313.15_EB-1._EB/373.15_EB))
   HUMIDITY = 40._EB*MIN( 1._EB , X_H2O_40_C/X_H2O_TMPA )
ENDIF
 
! Miscellaneous
 
MW_BACKGROUND = MW
MU_USER(0) = VISCOSITY
K_USER(0)  = CONDUCTIVITY
HCH    = C_HORIZONTAL
HCV    = C_VERTICAL
IF (ISOTHERMAL) RADIATION = .FALSE.
C_FORCED = C_FORCED*(1012._EB)*(1.8E-5_EB)**0.2_EB / (0.7_EB)**(2._EB/3._EB)
ASSUMED_GAS_TEMPERATURE = ASSUMED_GAS_TEMPERATURE + TMPM
TEX_ORI = TEXTURE_ORIGIN
 
IF (PRESSURE_CORRECTION) THEN
   SYNCHRONIZE = .TRUE.
   SYNC_TIME_STEP = .TRUE.
ENDIF
 
! Gravity ramp
 
I_RAMP_GX   = 0
I_RAMP_GY   = 0
I_RAMP_GZ   = 0
N_RAMP      = 0
IF (RAMP_GX/='null') CALL GET_RAMP_INDEX(RAMP_GX,'TIME',I_RAMP_GX)
IF (RAMP_GY/='null') CALL GET_RAMP_INDEX(RAMP_GY,'TIME',I_RAMP_GY)
IF (RAMP_GZ/='null') CALL GET_RAMP_INDEX(RAMP_GZ,'TIME',I_RAMP_GZ)
 
! Prandtl and Schmidt numbers
 
IF (DNS) THEN
   BAROCLINIC       = .TRUE.
   LES              = .FALSE.
   IF (PR<0._EB) PR = 0.7_EB
   IF (SC<0._EB) SC = 1.0_EB
ELSE
   IF (PR<0._EB) PR = 0.5_EB
   IF (SC<0._EB) SC = 0.5_EB
ENDIF
 
RSC = 1._EB/SC
RPR = 1._EB/PR
 
! Check for a restart file
 
APPEND = .FALSE.
IF (RESTART .AND. RESTART_CHID == CHID) APPEND = .TRUE.
IF (RESTART) NOISE  = .FALSE.
 
! Min and Max values of species and temperature
 
TMPMIN = TMPM
IF (LAPSE_RATE < 0._EB) TMPMIN = MIN(TMPMIN,TMPA+LAPSE_RATE*MESHES(1)%ZF)
TMPMAX = 3000._EB
YYMIN  = 0._EB
YYMAX  = 1._EB

IF (.NOT. SUPPRESSION .AND. CO_PRODUCTION) THEN
   WRITE(MESSAGE,'(A)')  'Cannot set SUPPRESSION=.FALSE. when CO_PRODUCTION=.TRUE.'
   CALL SHUTDOWN(MESSAGE)
ENDIF 
 
END SUBROUTINE READ_MISC

SUBROUTINE READ_DUMP

! Read parameters associated with output files
 
INTEGER :: N,ND
NAMELIST /DUMP/ RENDER_FILE,SMOKE3D,SMOKE3D_COMPRESSION,FLUSH_FILE_BUFFERS,MASS_FILE, &
                DT_CTRL,DT_PART,DT_MASS,DT_HRR,DT_DEVC,DT_PROF,DT_SLCF,DT_PL3D,DT_ISOF,DT_BNDF, &
                NFRAMES,DT_RESTART,DEBUG,TIMING,COLUMN_DUMP_LIMIT,MAXIMUM_DROPLETS,WRITE_XYZ,PLOT3D_QUANTITY
 
RENDER_FILE          = 'null'
NFRAMES              = 1000 
SMOKE3D              = .TRUE.
IF (TWO_D .OR. N_REACTIONS==0 .OR. SOLID_PHASE_ONLY .OR. .NOT. MIXTURE_FRACTION) SMOKE3D = .FALSE.
SMOKE3D_COMPRESSION  = 'RLE'
DEBUG                = .FALSE.
TIMING               = .FALSE.
FLUSH_FILE_BUFFERS   = .TRUE.
PLOT3D_QUANTITY(1) = 'TEMPERATURE'
PLOT3D_QUANTITY(2) = 'U-VELOCITY'
PLOT3D_QUANTITY(3) = 'V-VELOCITY'
PLOT3D_QUANTITY(4) = 'W-VELOCITY'
PLOT3D_QUANTITY(5) = 'HRRPUV'
IF (ISOTHERMAL) PLOT3D_QUANTITY(5) = 'PRESSURE'
MASS_FILE = .FALSE.
WRITE_XYZ            = .FALSE.
MAXIMUM_DROPLETS        = 500000
COLUMN_DUMP_LIMIT = .TRUE.  ! Limit csv files to 255 columns
 
DT_BNDF = -1.
DT_RESTART = 1000000._EB
DT_DEVC = -1.
DT_HRR  = -1.
DT_ISOF = -1.
DT_MASS = -1.
DT_PART = -1.
DT_PL3D = -1.
DT_PROF = -1.
DT_SLCF = -1.
DT_CTRL = -1.
 
REWIND(LU_INPUT)
DUMP_LOOP: DO 
   CALL CHECKREAD('DUMP',LU_INPUT,IOS)
   IF (IOS==1) EXIT DUMP_LOOP
   READ(LU_INPUT,DUMP,END=23,ERR=24,IOSTAT=IOS)
   24 IF (IOS>0) CALL SHUTDOWN('ERROR: Problem with DUMP line')
ENDDO DUMP_LOOP
23 REWIND(LU_INPUT)
IF (DT_BNDF < 0._EB) DT_BNDF = 2._EB * (T_END - T_BEGIN)/REAL(NFRAMES,EB) 
IF (DT_DEVC < 0._EB) DT_DEVC = (T_END - T_BEGIN)/REAL(NFRAMES,EB) 
IF (DT_HRR  < 0._EB) DT_HRR  = (T_END - T_BEGIN)/REAL(NFRAMES,EB)
IF (DT_ISOF < 0._EB) DT_ISOF = (T_END - T_BEGIN)/REAL(NFRAMES,EB) 
IF (DT_MASS < 0._EB) DT_MASS = (T_END - T_BEGIN)/REAL(NFRAMES,EB) 
IF (DT_PART < 0._EB) DT_PART = (T_END - T_BEGIN)/REAL(NFRAMES,EB) 
IF (DT_PL3D < 0._EB) DT_PL3D = (T_END - T_BEGIN)/5._EB 
IF (DT_PROF < 0._EB) DT_PROF = (T_END - T_BEGIN)/REAL(NFRAMES,EB) 
IF (DT_SLCF < 0._EB) DT_SLCF = (T_END - T_BEGIN)/REAL(NFRAMES,EB)
IF (DT_CTRL < 0._EB) DT_CTRL = (T_END - T_BEGIN)/REAL(NFRAMES,EB)

! Check Plot3D QUANTITIES

PLOOP: DO N=1,5
   DO ND=-N_OUTPUT_QUANTITIES,N_OUTPUT_QUANTITIES
      IF (PLOT3D_QUANTITY(N)==OUTPUT_QUANTITY(ND)%NAME) THEN
         PLOT3D_QUANTITY_INDEX(N) = ND
         IF (OUTPUT_QUANTITY(ND)%MIXTURE_FRACTION_ONLY .AND. .NOT.MIXTURE_FRACTION) THEN
            WRITE(MESSAGE,'(3A)')  'ERROR: PLOT3D quantity ',TRIM(PLOT3D_QUANTITY(N)), ' not appropriate for non-mixture fraction'
            CALL SHUTDOWN(MESSAGE)
         ENDIF 
         IF (OUTPUT_QUANTITY(ND)%SOLID_PHASE) THEN
            WRITE(MESSAGE,'(3A)') 'ERROR: PLOT3D quantity ',TRIM(PLOT3D_QUANTITY(N)), ' not appropriate for gas phase'
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         IF (OUTPUT_QUANTITY(ND)%PART_APPROPRIATE) THEN
            WRITE(MESSAGE,'(3A)')  'ERROR: PLOT3D quantity ',TRIM(PLOT3D_QUANTITY(N)), ' not appropriate as Plot3D'
            CALL SHUTDOWN(MESSAGE)
         ENDIF 
         CYCLE PLOOP
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(3A)') 'ERROR: PLOT3D quantity ', TRIM(PLOT3D_QUANTITY(N)),' not found'
   CALL SHUTDOWN(MESSAGE)
ENDDO PLOOP

END SUBROUTINE READ_DUMP


 
SUBROUTINE READ_SPEC
 
REAL(EB) :: MASS_FRACTION_0,MW, SIGMALJ,EPSILONKLJ,XVAP,VISCOSITY,CONDUCTIVITY,DIFFUSIVITY
INTEGER  :: N_SPEC_READ,N_SPEC_EXTRA,N_MIX,N,NN
LOGICAL  :: ABSORBING
NAMELIST /SPEC/ MASS_FRACTION_0,MW,FYI,ID,SIGMALJ, EPSILONKLJ,CONDUCTIVITY, VISCOSITY,DIFFUSIVITY, ABSORBING
 
! Zero out indices of various species
 
I_WATER   = 0
I_FUEL    = 0
I_PROG_F  = 0
I_PROG_CO = 0
I_CO2     = 0
I_CO      = 0
I_O2      = 0
I_SOOT    = 0
 
! Count SPEC lines and check for errors
 
N_SPECIES = 0
N_SPEC_READ = 0
REWIND(LU_INPUT)
COUNT_SPEC_LOOP: DO
   CALL CHECKREAD('SPEC',LU_INPUT,IOS) 
   IF (IOS==1) EXIT COUNT_SPEC_LOOP
   ID = 'null'
   READ(LU_INPUT,NML=SPEC,END=29,ERR=30,IOSTAT=IOS)
   IF (ID=='null') THEN
      WRITE(MESSAGE,'(A,I2,A)') 'ERROR: Species',N_SPECIES+1, 'needs a name (ID=...)'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   N_SPECIES = N_SPECIES + 1
   N_SPEC_READ = N_SPEC_READ + 1   
   IF (ID=='WATER VAPOR')     I_WATER = N_SPECIES
   IF (ID=='CARBON DIOXIDE')  I_CO2   = N_SPECIES
   IF (ID=='CARBON MONOXIDE') I_CO    = N_SPECIES      
   IF (ID=='OXYGEN')          I_O2    = N_SPECIES            
   IF (ID=='SOOT')            I_SOOT  = N_SPECIES            
   IF (MIXTURE_FRACTION .AND. (ID=='MIXTURE_FRACTION_1' .OR. ID=='MIXTURE_FRACTION_2' .OR. &
      (CO_PRODUCTION .AND. ID=='MIXTURE_FRACTION_3'))) N_SPECIES = N_SPECIES - 1
   30 IF (IOS>0) THEN
         WRITE(MESSAGE,'(A,I2)') 'ERROR: Problem with SPECies number',N_SPECIES+1
         CALL SHUTDOWN(MESSAGE)
      ENDIF
ENDDO COUNT_SPEC_LOOP
29 REWIND(LU_INPUT)
N_SPEC_EXTRA = N_SPECIES 
! Add other "SPECies" to the list
 
N_MIX = 0
DO N=1,N_REACTIONS