read.f

一个开源的火灾动力模拟的系统

      NU_CO   = CO_YIELD*MW_FUEL/MW_CO
      NU_CO2  = NU_CO2 - NU_SOOT - NU_CO
      NU_H2O  = NU_H2O - NU_H2
      NU_O2   = NU_O2 - NU_SOOT - 0.5*(NU_H2+NU_CO)
      NU_N2_AIR   = NU_O2*MW_O2*Y_N2_INFTY/(Y_O2_INFTY*MW_N2)
      NU_N2_DILUENT = NU_FUEL*MW_FUEL*Y_N2_INLET/(Y_F_INLET*MW_N2)
C
C Stoichiometric value of mixture fraction
C
      Z_F  = Y_O2_INFTY/
     .(Y_O2_INFTY+Y_F_INLET*NU_O2*MW_O2/MW_FUEL)
C
      IF (Y_F_INLET.LT.(1.-MW_N2*FUEL_N2/MW_FUEL)) 
     .         AUTOMATIC_Z = .FALSE.
      IF (DNS) AUTOMATIC_Z = .FALSE.
C
C Allocate state relation arrays
C
      ALLOCATE (MW_AVG(0:10000))
      CALL ChkMemErr('READ','MW_AVG',IZERO) ; MW_AVG=0.
      ALLOCATE (Y_STATE(0:10000,1:8))
      CALL ChkMemErr('READ','Y_STATE',IZERO) ; Y_STATE=0.
      ALLOCATE (RSUM_MF(0:10000))
      CALL ChkMemErr('READ','RSUM_MF',IZERO) ; RSUM_MF=0.
C
C Compute Y_STATE and MW_AVG
C
      Y_STATE(0,1) = 0.
      Y_STATE(0,2) = Y_O2_INFTY
      Y_STATE(0,3) = Y_N2_INFTY
      MW_AVG(0) = 1/(Y_N2_INFTY/MW_N2+Y_O2_INFTY/MW_O2)
C
      Y_STATE(10000,1) = Y_F_INLET
      Y_STATE(10000,2) = 0.
      Y_STATE(10000,3) = Y_N2_INLET
      MW_AVG(10000) = 1/(Y_F_INLET/MW_FUEL+Y_N2_INLET/MW_N2)
C
      DZZ=1./10000.
      CHI=1.0
C
      DO I=1,9999
      ZZ=I*DZZ
C
      S_TERM = NU_O2*MW_O2/MW_FUEL
      Z_TERM = ZZ*(S_TERM*Y_F_INLET+Y_O2_INFTY) - Y_O2_INFTY
      ETA = ((S_TERM-Z_TERM)*MW_FUEL-Z_TERM*NU_N2_DILUENT*MW_N2)/
     .      (NU_O2*MW_O2*(Z_TERM+1.)+NU_N2_AIR*MW_N2*Z_TERM)
      TOTAL_MASS = MW_FUEL + ETA*(NU_O2*MW_O2 + NU_N2_AIR*MW_N2)+
     .             NU_N2_DILUENT*MW_N2
C
      Y_STATE(I,1) = MW_FUEL*MAX(0._EB,1.-ETA)/TOTAL_MASS 
     .             + MW_FUEL*MIN(1._EB,ETA)*(1.-CHI)/TOTAL_MASS
      Y_STATE(I,2) = MW_O2*NU_O2*MAX(0._EB,ETA-1.)/TOTAL_MASS
     .             + MW_O2*NU_O2*MIN(1._EB,ETA)*(1.-CHI)/TOTAL_MASS
      Y_STATE(I,3) = MW_N2  *(NU_N2_AIR*ETA+NU_N2_DILUENT)/TOTAL_MASS
      Y_STATE(I,4) = MW_H2O *NU_H2O *MIN(1._EB,ETA)*CHI/TOTAL_MASS
      Y_STATE(I,5) = MW_CO2 *NU_CO2 *MIN(1._EB,ETA)*CHI/TOTAL_MASS
      Y_STATE(I,6) = MW_CO  *NU_CO  *MIN(1._EB,ETA)*CHI/TOTAL_MASS
      Y_STATE(I,7) = MW_H2  *NU_H2  *MIN(1._EB,ETA)*CHI/TOTAL_MASS
      Y_STATE(I,8) = MW_C   *NU_SOOT*MIN(1._EB,ETA)*CHI/TOTAL_MASS

      MW_AVG(I) = 1./(Y_STATE(I,1)/MW_FUEL +
     .                Y_STATE(I,2)/MW_O2   +
     .                Y_STATE(I,3)/MW_N2   +
     .                Y_STATE(I,4)/MW_H2O  +
     .                Y_STATE(I,5)/MW_CO2  +
     .                Y_STATE(I,6)/MW_CO   +
     .                Y_STATE(I,7)/MW_H2   +
     .                Y_STATE(I,8)/MW_C)
C
      ENDDO
C
C Compute R0/MW_AVG
C
      RSUM_MF = R0/MW_AVG
C
      END SUBROUTINE STATE_RELATIONSHIPS
C
C
      SUBROUTINE SET_REAC_DEFAULTS
C
      ID   = 'null'
      BOF  = 0.          ! Pre-exponential Factor (cm**3/mol-s)
      E    = 0.          ! Activation Energy (kJ/kmol)
      XNO  = 1.
      XNF  = 1.
      FUEL = 'nullspecies'
      DELTAH = -1.
      DELTAH_FUEL = 40000.  ! Energy per unit mass fuel consumed (kJ/kg)
      HRRPUA_SHEET = 200.   ! Max HRR per unit area of flame sheet
      TMP_LOWER  = 325.  ! Lower Temp for One Step Reaction
      EPUMO2     = 13100.   ! Energy per unit mass oxygen (kJ/kg)
      IF (LES) THEN
         RADIATIVE_FRACTION = 0.35
         SOOT_YIELD = 0.01
         ELSE
         RADIATIVE_FRACTION = 0.0
         SOOT_YIELD = 0.
         ENDIF
      DTSAM      = TWFIN/REAL(NFRAMES)
      Y_O2_INFTY = 0.23
      X_O2_LL    = 0.15
      CRITICAL_FLAME_TEMPERATURE = 1427. 
      Y_F_INLET  = 1.
      NU_O2      = 5.       ! Default fuel is propane
      NU_CO2     = 3.
      NU_H2O     = 4.
      NU_FUEL    = 1.
      MW_FUEL    = 44.
      MW_O2      = 32.
      MW_CO2     = 44.
      MW_C       = 12.
      MW_N2      = 28.
      MW_H2O     = 18.
      MW_CO      = 28.
      MW_H2      = 2.
      FUEL_N2    = 0.
      CO_YIELD   = -1.
      H2_YIELD   =  0.
      MASS_EXTINCTION_COEFFICIENT = 7600.     ! m2/kg
      VISIBILITY_FACTOR = 3.
      MAXIMUM_VISIBILITY = 30.
      Z_CONSTANT = 1.0
C
      END SUBROUTINE SET_REAC_DEFAULTS
C
      END SUBROUTINE READ_SPEC
C
C
      SUBROUTINE READ_PART
C
      INTEGER :: NUMBER_INITIAL_DROPLETS,NIP,NISP,NSPINS,NPDIM,
     .           NSPDIM,SAMPLING_FACTOR,DROPLETS_PER_SECOND,
     .           NPSAM,NSPSAM
      REAL(EB) :: SPECIFIC_HEAT,VAPORIZATION_TEMPERATURE,
     .            MELTING_TEMPERATURE,MASS_PER_VOLUME,DIAMETER,
     .            GAMMA_D,WMPUV,DTSAM,AGE,DELAY,INITIAL_TEMPERATURE
      CHARACTER(30) :: SPEC_ID,HEAT_ACTIVATE
      LOGICAL :: MASSLESS,STATIC,FUEL,WATER,TREE,FUEL_DROPLETS
      NAMELIST /PART/ NUMBER_INITIAL_DROPLETS,QUANTITY,FYI,
     .                DROPLETS_PER_SECOND,MASS_PER_VOLUME,
     .                SAMPLING_FACTOR,NPSAM,NSPSAM,
     .                ID,STATIC,MASSLESS,FUEL,WATER,TREE,
     .                DENSITY,VAPORIZATION_TEMPERATURE,
     .                SPECIFIC_HEAT,HEAT_OF_VAPORIZATION,
     .                MELTING_TEMPERATURE,DIAMETER,GAMMA_D,WMPUV,
     .                HEAT_OF_COMBUSTION,DTPAR,DTSPAR,
     .                NIP,NISP,NSPINS,DTSAM,AGE,NPDIM,NSPDIM,
     .                FUEL_DROPLETS,SPEC_ID,INITIAL_TEMPERATURE,XB
      EQUIVALENCE (SAMPLING_FACTOR,NPSAM,NSPSAM)
      EQUIVALENCE (NUMBER_INITIAL_DROPLETS,NISP,NIP)
      EQUIVALENCE (WMPUV,MASS_PER_VOLUME)
      EQUIVALENCE (FUEL,FUEL_DROPLETS)
C
      REAL(EB) T_ACTIVATE,T_DEACTIVATE
      NAMELIST /SPRK/ XYZ,MAKE,T_ACTIVATE,T_DEACTIVATE,FYI,
     .                ORIENTATION,ROTATION,PART_ID,LABEL,DELAY,
     .                HEAT_ACTIVATE
C
      WATER_EVAPORATION = .FALSE.
      FUEL_EVAPORATION  = .FALSE.
      EVAPORATION       = .FALSE.
C
C Look for sprinklers, if found WATER VAPOR will be a species
C
      NSPR = 0
      SPRLOOP: DO
      PART_ID = 'null'
      CALL CHECKREAD('SPRK',LU5,IOS) ; IF (IOS.EQ.1) EXIT SPRLOOP
      READ(LU5,SPRK,END=7,ERR=8,IOSTAT=IOS)
      NSPR = NSPR + 1
      IF (PART_ID.EQ.'null') WATER_EVAPORATION = .TRUE.
    8 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with SPRK line')
      ENDDO SPRLOOP
    7 REWIND(LU5)
C
C Determine total number of PART lines in the input file
C
      NPC = 0
      COUNT_PART_LOOP: DO
      CALL CHECKREAD('PART',LU5,IOS) 
      IF (IOS.EQ.1) EXIT COUNT_PART_LOOP
      READ(LU5,NML=PART,END=219,ERR=220,IOSTAT=IOS)
      NPC = NPC + 1
  220 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with PART line')
      ENDDO COUNT_PART_LOOP
  219 REWIND(LU5)
C
C Add two special PART classes, sprinkler droplets and tracers
C
      NPC = NPC + 2
C
C Allocate quantities for PART classes
C
      IPART = 0
C
      ALLOCATE(LAGRANGIAN(NPC),STAT=IZERO)
      CALL ChkMemErr('READ','LAGRANGIAN',IZERO) 

      LAGRANGIAN(1:NPC)%ADJUST_EVAPORATION = 1. 
      LAGRANGIAN(1:NPC)%SPECIES_INDEX = 0
      LAGRANGIAN(1:NPC)%SPECIES = ' '
      LAGRANGIAN(1:NPC)%COLOR_INDEX = 0
C
      READ_PART_LOOP: DO N=1,NPC
C
      LP=>LAGRANGIAN(N)
C
      DENSITY                  = 1000.   ! kg/m3
      MASS_PER_VOLUME          = 1.      ! kg/m3
      VAPORIZATION_TEMPERATURE = 100.0   ! C
      INITIAL_TEMPERATURE      = TMPA - TMPM ! C
      MELTING_TEMPERATURE      = TMPM - TMPM ! C
      SPECIFIC_HEAT            = 4.184   ! kJ/kg-K
      HEAT_OF_VAPORIZATION     = 2259.   ! kJ/kg
      HEAT_OF_COMBUSTION       = -1.     ! kJ/kg
      DROPLETS_PER_SECOND      = 1000
      DIAMETER                 = 100.    ! mu-m
      GAMMA_D                  = 2.4
      AGE                      = 1.E6    ! s
C
      ID                       = 'null'
      QUANTITY                 = 'null'
      SPEC_ID                  = 'null'
      SAMPLING_FACTOR          = -1      
      NUMBER_INITIAL_DROPLETS  = 0
      XB                       = 0.
      FUEL                     = .FALSE.
      WATER                    = .FALSE.
      STATIC                   = .FALSE.
      MASSLESS                 = .FALSE.
      TREE                     = .FALSE.
C
C Read the PART line
C
      IF (N.LT.NPC-1) THEN
      CALL CHECKREAD('PART',LU5,IOS) ; IF (IOS.EQ.1) EXIT READ_PART_LOOP
      READ(LU5,PART,END=25,IOSTAT=IOS)
C
      IF (NUMBER_INITIAL_DROPLETS.GT.0 .AND. RESTART)
     .    NUMBER_INITIAL_DROPLETS =  0
      ENDIF
C
      IF (N.EQ.NPC-1) THEN
      MASSLESS = .TRUE.
      ID       = 'tracer particles'
      ENDIF
C
      IF (SAMPLING_FACTOR.LE.0) THEN
      IF (MASSLESS) THEN
         SAMPLING_FACTOR = 1
         ELSE
         SAMPLING_FACTOR = 10
         ENDIF
         ENDIF
C
      IF (N.EQ.NPC) THEN
      WATER = .TRUE.
      ID    = 'water droplets'
      SPEC_ID = 'WATER VAPOR'
      AGE   = 60.
      ENDIF
C
      IF (TREE) FUEL   = .TRUE.
      IF (TREE) STATIC = .TRUE.
C
      IF (MASSLESS) DIAMETER = 0.
C
      LP%CLASS_NAME = ID
C
      IF (QUANTITY.EQ.'TEMPERATURE')       IPART = -1
      IF (QUANTITY.EQ.'HRRPUV')            IPART = -2
      IF (QUANTITY.EQ.'DIAMETER')          IPART = -3
      IF (QUANTITY.EQ.'VELOCITY')          IPART = -4
      IF (QUANTITY.EQ.'DROPLET_PHASE')     IPART = -5
      IF (QUANTITY.EQ.'AGE')               IPART = -6
      IF (QUANTITY.EQ.'BLACK')    LP%COLOR_INDEX = 0
      IF (QUANTITY.EQ.'YELLOW')   LP%COLOR_INDEX = 1
      IF (QUANTITY.EQ.'BLUE')     LP%COLOR_INDEX = 2
      IF (QUANTITY.EQ.'RED')      LP%COLOR_INDEX = 3
      IF (QUANTITY.EQ.'GREEN')    LP%COLOR_INDEX = 4
      IF (QUANTITY.EQ.'MAGENTA')  LP%COLOR_INDEX = 5
      IF (QUANTITY.EQ.'CYAN')     LP%COLOR_INDEX = 6
      IF (QUANTITY.EQ.'WHITE')    LP%COLOR_INDEX = 7
C
      IF (QUANTITY.EQ.'null' .AND. WATER) LP%COLOR_INDEX = 2
      IF (QUANTITY.EQ.'null' .AND. TREE)  LP%COLOR_INDEX = 4
C
      LP%SAMPLING = SAMPLING_FACTOR
C
      IF (NUMBER_INITIAL_DROPLETS.GT.0) DROPLET_FILE  = .TRUE.
C
      IF (NUMBER_INITIAL_DROPLETS.GT.0) THEN
         DO I=1,5,2
         IF (XB(I).GT.XB(I+1)) THEN
            DUMMY   = XB(I)
            XB(I)   = XB(I+1)
            XB(I+1) = DUMMY
            ENDIF
         ENDDO
         ENDIF
C
      LP%X1 = XB(1)
      LP%X2 = XB(2)
      LP%Y1 = XB(3)
      LP%Y2 = XB(4)
      LP%Z1 = XB(5)
      LP%Z2 = XB(6)
C
      LP%N_INITIAL = NUMBER_INITIAL_DROPLETS
      LP%N_INSERT  = DROPLETS_PER_SECOND*DTSPAR
C
      IF (DIAMETER.GT.0.) THEN
      ALLOCATE(LP%CDF(0:NDC),STAT=IZERO)
      CALL ChkMemErr('READ','CDF',IZERO)
      ALLOCATE(LP%R_CDF(0:NDC),STAT=IZERO)
      CALL ChkMemErr('READ','R_CDF',IZERO)
      ALLOCATE(LP%IL_CDF(NSTRATA),STAT=IZERO)
      CALL ChkMemErr('READ','IL_CDF',IZERO)
      ALLOCATE(LP%IU_CDF(NSTRATA),STAT=IZERO)
      CALL ChkMemErr('READ','IU_CDF',IZERO)
      ALLOCATE(LP%W_CDF(NSTRATA),STAT=IZERO)
      CALL ChkMemErr('READ','W_CDF',IZERO)
      ENDIF
C
C     Assign property data to LAGRANGIAN class
C
      LP%DIAMETER = DIAMETER*1.E-6
      LP%GAMMA    = GAMMA_D
      LP%SIGMA    = 1.15/GAMMA_D
      LP%TMP_INITIAL = INITIAL_TEMPERATURE + TMPM
      LP%C_P      = SPECIFIC_HEAT*1000.
      LP%H_V_0    = HEAT_OF_VAPORIZATION*1000.
      LP%DELTAH   = HEAT_OF_COMBUSTION*1000.
      LP%TMP_V    = VAPORIZATION_TEMPERATURE + TMPM
      LP%DENSITY  = DENSITY
      LP%FTPR     = FOTH*PI*DENSITY
      LP%MASS_PER_VOLUME = MASS_PER_VOLUME
      LP%TMP_MELT = MELTING_TEMPERATURE + TMPM
      LP%MASSLESS = MASSLESS
      LP%LIFETIME = AGE
      LP%TREE     = TREE
      LP%FUEL     = FUEL
      LP%WATER    = WATER
      LP%STATIC   = STATIC
      IF (FUEL)  SPEC_ID = 'MIXTURE_FRACTION'
      IF (WATER) SPEC_ID = 'WATER VAPOR'
      LP%SPECIES  = SPEC_ID
C
      IF (N.LT.NPC-1 .AND. FUEL)  FUEL_EVAPORATION = .TRUE.
      IF (N.LT.NPC-1 .AND. WATER) WATER_EVAPORATION = .TRUE.
C
      ENDDO READ_PART_LOOP
   25 REWIND(LU5)
C
      IF (FUEL_EVAPORATION .OR. WATER_EVAPORATION) EVAPORATION=.TRUE.
C
      END SUBROUTINE READ_PART
C
C
      SUBROUTINE READ_TREE
C
      INTEGER :: IPC,N_TREES_0,NM,NN,N
      REAL(EB) :: CANOPY_WIDTH,CANOPY_BASE_HEIGHT,TREE_HEIGHT
      NAMELIST /TREE/ XYZ,CANOPY_WIDTH,CANOPY_BASE_HEIGHT,
     .                     TREE_HEIGHT,PART_ID
C
C Read the TREE lines to determine how many cone shaped trees
C there will be
C
      N_TREES = 0
      COUNT_TREE_LOOP: DO
      CALL CHECKREAD('TREE',LU5,IOS) 
      IF (IOS.EQ.1) EXIT COUNT_TREE_LOOP
      READ(LU5,NML=TREE,END=11,ERR=12,IOSTAT=IOS)
      N_TREES = N_TREES + 1
   12 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with TREE line')
      ENDDO COUNT_TREE_LOOP
   11 REWIND(LU5)
C
C Sequentially read the CONE_TREE namelist to get shape and size
C parameters for each tree.
C
      IF (N_TREES.EQ.0) RETURN
C
      ALLOCATE(CANOPY_W(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','CANOPY_W',IZERO)
      ALLOCATE(CANOPY_B_H(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','CANOPY_B_H',IZERO)
      ALLOCATE(TREE_H(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','TREE_H',IZERO)
      ALLOCATE(X_TREE(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','X_TREE',IZERO)
      ALLOCATE(Y_TREE(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','Y_TREE',IZERO)
      ALLOCATE(Z_TREE(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','Z_TREE',IZERO)
      ALLOCATE(TREE_PARTICLE_CLASS(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','TREE_PARTICLE_CLASS',IZERO)
      ALLOCATE(TREE_MESH(N_TREES),STAT=IZERO)
      CALL ChkMemErr('READ','TREE_MESH',IZERO)
C
      N_TREES_0 = N_TREES
      N = 0
C
      CONE_LOOP: DO NN=1,N_TREES_0
      N = N + 1
C
      PART_ID = 'null'
C
      CALL CHECKREAD('TREE',LU5,IOS) ; IF (IOS.EQ.1) EXIT CONE_LOOP
      READ(LU5,TREE,END=25,IOSTAT=IOS)
C
      MESH_LOOP: DO NM=1,NMESHES
      IF (XYZ(1).GE.MESH(NM)%XS .AND. XYZ(1).LE.MESH(NM)%XF .AND.
     .    XYZ(2).GE.MESH(NM)%YS .AND. XYZ(2).LE.MESH(NM)%YF .AND.
     .    XYZ(3).GE.MESH(NM)%ZS .AND. XYZ(3).LE.MESH(NM)%ZF) THEN
         TREE_MESH(N) = NM
         EXIT MESH_LOOP
         ENDIF
      IF (NM.EQ.NMESHES) THEN
         N    = N-1
         N_TREES = N_TREES - 1
         CYCLE CONE_LOOP
         ENDIF
      EN