read.f

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

      THICKEN_OBSTRUCTIONS = .FALSE.
      DROP_VERTICAL_VELOCITY   = 0.20  
      DROP_HORIZONTAL_VELOCITY = 0.50  
      DEBUG                = .FALSE.
      TIMING               = .FALSE.
      FLUSH_FILE_BUFFERS   = .TRUE.
C
      MISC_LOOP: DO 
      CALL CHECKREAD('MISC',LU5,IOS) ; IF (IOS.EQ.1) EXIT MISC_LOOP
      READ(LU5,MISC,END=23,ERR=24,IOSTAT=IOS)
   24 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with MISC line')
      ENDDO MISC_LOOP
   23 REWIND(LU5)
C
C Look to see if DATABASE exists
C
      IF (DATABASE.EQ.'null' .AND. DATABASE_DIRECTORY.NE.'null') THEN
      DATABASE = TRIM(DATABASE_DIRECTORY)//'database4.data'
      ENDIF
C
      NFILES = 1
      IF (DATABASE.NE.'null') THEN
      INQUIRE(FILE=DATABASE,EXIST=EX)
      IF (.NOT.EX) THEN
         WRITE(MESSAGE,'(A)') 'ERROR: DATABASE '//TRIM(DATABASE)//
     .                        ' not found'
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      NFILES = 2
      OPEN(LU80,FILE=DATABASE,FORM='FORMATTED',STATUS='OLD',ERR=190,
     .     ACTION='READ',IOSTAT=IOS)
  190 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem opening DATABASE')
      ENDIF
C
      H0    = 0.5*(U0**2+V0**2+W0**2)
      TMPA  = TMPA + TMPM
      TMPA4 = TMPA**4
C
      IF (TMPO.GT.-1000.) THEN
         TMPO = TMPO + TMPM
      ELSE 
         TMPO = TMPA
      ENDIF
C
C Humidity (40% by default, but limited for high temps)
C
      IF (HUMIDITY.LT.0.)
     .HUMIDITY=0.4*MIN(1._EB,EXP(-4890.5*(1/TMPA-1/373.15))/0.0812)
C
C Miscellaneous
C
      IF (DENSITY.GT.0.) MW = TMPA*DENSITY*R0/PINF
      MW_BACKGROUND = MW
      MU_USER(0) = VISCOSITY
      K_USER(0)  = THERMAL_CONDUCTIVITY
      HCH    = C_HORIZONTAL
      HCV    = C_VERTICAL
      SURFNAME(0) = SURF_DEFAULT
      IF (ISOTHERMAL) RADIATION = .FALSE.
      C_FORCED = C_FORCED*(1012.)*(1.8E-5)**0.2 / (0.7)**(2./3.)
      TEX_ORI = TEXTURE_ORIGIN
      IF (GAUGE_TEMPERATURE.LT.0.) THEN
         TMP_GAUGE = TMPA
      ELSE
         TMP_GAUGE = GAUGE_TEMPERATURE + TMPM
      ENDIF
C
      IF (DTSAM_PART.LT.0.) THEN
      WPAR = TWFIN/REAL(NFRAMES,EB)
      ELSE
      WPAR = DTSAM_PART
      ENDIF
C
C Background Pressure Stuff
C
      P0_MAX = PINF + MAX_OVER_PRESSURE*PINF
C
C Gravity ramp
C
      IRAMPG = 0
      NRAMP  = 0
      IF (RAMP_G.NE.'null') THEN
      NRAMP = NRAMP + 1
      RAMPID(NRAMP) = RAMP_G
      RAMPTYPE(NRAMP) = 'TIME'
      IRAMPG = NRAMP
      ENDIF
C
C Prandtl and Schmidt numbers
C
      IF (DNS) THEN
         BAROCLINIC = .TRUE.
         LES = .FALSE.
         IF (PR.LT.0.) PR = 0.7
         IF (SC.LT.0.) SC = 1.0
         ELSE
         IF (PR.LT.0.) PR = 0.5
         IF (SC.LT.0.) SC = 0.5
         ENDIF
C
      RSC = 1./SC
      RPR = 1./PR
C
C Check for a restart file
C
      APPEND = .FALSE.
      IF (RESTART) APPEND = .TRUE.
C
C Min and Max values of species and temperature
C
      TMPMIN = MIN(TMPA,TMPO)
      TMPMAX = 3000.
      YYMIN  = 0.
      YYMAX  = 1.
C
      END SUBROUTINE READ_MISC
C
C
      SUBROUTINE READ_SPEC
C
      REAL(EB) EPSIJ,SIGMA2,AMW,MASS_FRACTION_0,MW,OMEGA,TSTAR,
     .         EPSK(0:5),SIG(0:5),SIGMALJ,EPSILONKLJ,NU,MU_N,K_N,D_N0,
     .         ZZ,DZZ,TE,WGT,EPSK_N,SIG_N,MW_N
      REAL(EB) CP_N2,CP_F,CP_CO2,CP_O2,CP_H2O,XVAP,MAXIMUM_VISIBILITY
      REAL(EB) RCON_MF
      INTEGER IPC,NN
      LOGICAL READ_MIX
      INTEGER  ITMP,IYY,NSPEC_READ
      CHARACTER(26) FUEL,STATE_SPECIES(10)
      NAMELIST /SURF/ VBC,TMPIGN,TMPWAL,TMPWAL0,ALPHA,C_P,KS,DELTA,
     .                C_DELTA_RHO,MASS_FRACTION,VEL,VEL_T,NPPC,
     .                E_COEFFICIENT,HEAT_FLUX,WALL_POINTS,
     .                TIGN,TAU_Q,TAU_V,RAMP_Q,SURFACE_DENSITY,TAU_MF,
     .                RAMP_MF,PART_ID,RAMP_V,VOLUME_FLUX,RAMP_KS,
     .                RAMP_KS_CHAR,RAMP_C_P_CHAR,
     .                PROFILE,PLE,Z0,ID,AKA,MASS_FLUX,RAMP_C_P,
     .                FYI,POROSITY,PAPER_MODEL,
     .                BACKING,TMP_BACK,HRRPUA,EMISSIVITY,PHASE,RADIUS,
     .                HEAT_OF_VAPORIZATION,DENSITY,IN_DEPTH_COEFFICIENT,
     .                HEAT_OF_ABLATION,ABLATION_TEMPERATURE,
     .                ABLATION_RATE,ADIABATIC,HEAT_OF_COMBUSTION,
     .                BURNING_RATE_MAX,PARTICLES,PARTICLE_COLOR,
     .                TEXTURE_MAP,TEXTURE_WIDTH,TEXTURE_HEIGHT,RGB,RGB4,
     .                MASS_FLUX_CRITICAL,A,
     .                E,MOISTURE_FRACTION,FUEL_FRACTION,
     .                CHAR_DENSITY, C_P_CHAR, KS_CHAR,BURN_AWAY,LEAKING,
     .                DX_SOLID,EXTERNAL_FLUX,TMPEVAP
      NAMELIST /SPEC/ MASS_FRACTION_0,MW,FYI,ID,SIGMALJ,
     .                EPSILONKLJ,NU,DENSITY,THERMAL_CONDUCTIVITY,
     .                VISCOSITY,DIFFUSION_COEFFICIENT
      NAMELIST /REAC/ E,BOF,XNO,XNF,DELTAH,FYI,FUEL,EPUMO2,ID,
     .                RADIATIVE_FRACTION,HRRPUA_SHEET,TMP_LOWER,
     .                Y_O2_INFTY,Y_F_INLET,NU_FUEL,NU_O2,
     .                NU_H2O,NU_CO2,MW_FUEL,FUEL_N2,
     .                H2_YIELD,SOOT_YIELD,CO_YIELD,
     .                MASS_EXTINCTION_COEFFICIENT,DTSAM,
     .                X_O2_LL,Z_CONSTANT,VISIBILITY_FACTOR,
     .                MAXIMUM_VISIBILITY,
     .                CRITICAL_FLAME_TEMPERATURE
C
C Scan all active SURF lines to see if mixture fraction combustion
C model is desired
C
      IWATER  = 0
      IFUEL   = 0
      IOXYGEN = 0     
      ICO2    = 0
      MIXTURE_FRACTION = .FALSE.
C
      IF (FUEL_EVAPORATION) MIXTURE_FRACTION = .TRUE.
C
      N = -1
      SURF_LOOP: DO 
      N = N+1
      IF (N.GT.NBT) EXIT SURF_LOOP
C
      AKA_NAME(N) = 'null'
      IF (SURFNAME(N).EQ.'INERT') CYCLE SURF_LOOP
C
      FILE_LOOP: DO ITER=1,NFILES
C
      IF (ITER.EQ.1) LUDUM = LU5
      IF (ITER.EQ.2) LUDUM = LU80
      REWIND(LUDUM)
C
      LOOP1: DO
      CALL CHECKREAD('SURF',LUDUM,IOS) ; IF (IOS.EQ.1) CYCLE FILE_LOOP
      CALL SET_SURF_DEFAULTS
      READ(LUDUM,SURF,ERR=34,IOSTAT=IOS)
C
      NN = -1
      FIND_BACKING: DO
      NN = NN+1
      IF (NN.GT.NBT) EXIT FIND_BACKING
      IF (BACKING.EQ.SURFNAME(NN)) EXIT FIND_BACKING
      IF (NN.EQ.NBT .AND. BACKING.NE.'VOID' .AND. BACKING.NE.'EXPOSED'
     .    .AND. BACKING.NE.'INSULATED') THEN
         NBT = NBT+1
         SURFNAME(NBT) = BACKING
         EXIT FIND_BACKING
         ENDIF
      ENDDO FIND_BACKING

      IF (ID.EQ.SURFNAME(N) .OR. ID.EQ.AKA_NAME(N)) THEN
         IF (AKA.NE.'null') THEN
             AKA_NAME(N) = AKA
             CYCLE LOOP1
             ENDIF
         IF (HRRPUA.GT.0. .OR. HEAT_OF_VAPORIZATION.GT.0.) THEN
             MIXTURE_FRACTION = .TRUE.
             ENDIF 
         EXIT FILE_LOOP
         ENDIF
   34 IF (IOS.GT.0) THEN
         WRITE(MESSAGE,'(A)') 'ERROR: Problem with SURF '//
     .       TRIM(SURFNAME(N))//' or some SURF line preceding it'
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      ENDDO LOOP1
C
      ENDDO FILE_LOOP
C
      IF (ID.NE.SURFNAME(N) .AND. ID.NE.AKA_NAME(N)) THEN
      WRITE(MESSAGE,'(A)') 'ERROR: SURF '//TRIM(SURFNAME(N))//
     .                     ' not found'
      CALL SHUTDOWN(MESSAGE)
      ENDIF
C
      REWIND(LU5)
      IF (DATABASE.NE.'null') REWIND(LU80)
C
      ENDDO SURF_LOOP
C
C Check SPEC lines for errors and determine number of extra species
C
      ALLOCATE(SPECIES_ID(0:5),STAT=IZERO)
      CALL ChkMemErr('READ','SPECIES_ID',IZERO)
      SPECIES_ID(0) = BACKGROUND_SPECIES
      SPECIES_ID(1:5) = 'nullspecies'
C
      NSPEC = 0
      READ_MIX = .FALSE.
C
      SPECLOOP: DO
      CALL CHECKREAD('SPEC',LU5,IOS) ; IF (IOS.EQ.1) EXIT SPECLOOP
      ID = 'nullspecies'
      NU = 0.
      READ(LU5,NML=SPEC,END=29,ERR=30,IOSTAT=IOS)
C
      IF (ID.EQ.'nullspecies') THEN
         WRITE(MESSAGE,'(A,I2,A)') 'ERROR: Species',N,
     .                             'needs a name (ID=...)'
         CALL SHUTDOWN(MESSAGE)
         ENDIF
C
      NSPEC = NSPEC + 1
      SPECIES_ID(NSPEC) = ID
C
      IF (NU.NE.0.) MIXTURE_FRACTION = .FALSE.
      IF (ID.EQ.'MIXTURE_FRACTION') THEN
         READ_MIX         = .TRUE.
         MIXTURE_FRACTION = .TRUE.
         IFUEL            = NSPEC
         ENDIF
      IF (ID.EQ.'WATER VAPOR') IWATER = NSPEC
      IF (ID.EQ.'OXYGEN') IOXYGEN = NSPEC
      IF (ID.EQ.'CARBON DIOXIDE') ICO2 = NSPEC
C
   30 IF (IOS.GT.0) THEN
         WRITE(MESSAGE,'(A,I2)') 
     .       'ERROR: Problem with SPECies number',NSPEC+1
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      ENDDO SPECLOOP
   29 REWIND(LU5)
C
      NSPEC_READ = NSPEC
C
C Check if MIXTURE_FRACTION needs to be added to SPEC list
C
      IF (MIXTURE_FRACTION .AND. IFUEL.EQ.0) THEN
         NSPEC = NSPEC + 1
         SPECIES_ID(NSPEC) = 'MIXTURE_FRACTION'
         IFUEL = NSPEC
         ENDIF
C
C Check if water vapor should be included in the calculation
C
      IF (WATER_EVAPORATION .AND. IWATER.EQ.0) THEN
         NSPEC = NSPEC + 1
         SPECIES_ID(NSPEC) = 'WATER VAPOR'
         IWATER = NSPEC
         ENDIF
C
C Allocate species-related arrays
C
      ALLOCATE(NUN(0:NSPEC),STAT=IZERO)
      CALL ChkMemErr('READ','NUN',IZERO)
C
      ALLOCATE(YY0(0:NSPEC),STAT=IZERO)
      CALL ChkMemErr('READ','YY0',IZERO)
      ALLOCATE(RCON(0:NSPEC),STAT=IZERO)
      CALL ChkMemErr('READ','RCON',IZERO)
C
      IF (MIXTURE_FRACTION) THEN
      ALLOCATE(CP(0:100,500),STAT=IZERO)
      CALL ChkMemErr('READ','CP',IZERO)
      ALLOCATE(RCP(0:100,500),STAT=IZERO)
      CALL ChkMemErr('READ','RCP',IZERO)
      ALLOCATE(HH(0:100,0:500),STAT=IZERO)
      CALL ChkMemErr('READ','HH',IZERO)
      ALLOCATE(RCON_STATE(0:10),STAT=IZERO)
      CALL ChkMemErr('READ','RCON_STATE',IZERO)
      ALLOCATE(MWN(0:10),STAT=IZERO)
      CALL ChkMemErr('READ','MWN',IZERO)
      ELSE
      ALLOCATE(CP(0:NSPEC,500),STAT=IZERO)
      CALL ChkMemErr('READ','CP',IZERO)
      ALLOCATE(RCP(0:NSPEC,500),STAT=IZERO)
      CALL ChkMemErr('READ','RCP',IZERO)
      ALLOCATE(HH(0:NSPEC,0:500),STAT=IZERO)
      CALL ChkMemErr('READ','HH',IZERO)
      ALLOCATE(MWN(0:NSPEC),STAT=IZERO)
      CALL ChkMemErr('READ','MWN',IZERO)
      ENDIF
C
      YY0     = 0.
      MWN     = 0.
      MWN(0)  = MW_BACKGROUND
      EPSK    = 0.
      SIG     = 0.
      NUN     = 0.
C
C Initialize variables for MIX_FRAC and WATER_VAPOR, if necessary
C
      IF (MIXTURE_FRACTION .AND. .NOT.READ_MIX) THEN
         YY0(IFUEL) = 0.
         MWN(IFUEL) = 0.
         EPSK(IFUEL) = 0.
         SIG(IFUEL)  = 0.
         NUN(IFUEL)  = 0.
         ENDIF
C
      IF (IWATER.GT.0) THEN
         LP=>LAGRANGIAN(NPC)
         XVAP  = MIN(1._EB,EXP(LP%H_V_0*18./R0*(1./LP%TMP_V-1./
     .                     MIN(TMPA,LP%TMP_V))))
         YY0(IWATER) = HUMIDITY*XVAP/(29./18. + (1.-29./18.)*XVAP)
         MWN(IWATER) = 0.
         EPSK(IWATER) = 0.
         SIG(IWATER)  = 0.
         NUN(IWATER)  = 0.
         ENDIF
C
C Read SPEC info for species other than MIXTURE_FRACTION and WATER VAPOR
C
      SPEC_LOOP: DO N=1,NSPEC_READ
C
      MASS_FRACTION_0 = -1.
      MW  = 0.
      DENSITY = -1.
      EPSILONKLJ = 0.
      SIGMALJ    = 0.
      NU         = 0.
      THERMAL_CONDUCTIVITY  = -1.
      VISCOSITY     = -1.
      DIFFUSION_COEFFICIENT = -1.
C
      CALL CHECKREAD('SPEC',LU5,IOS) ; IF (IOS.EQ.1) CYCLE SPEC_LOOP
      READ(LU5,SPEC,END=31,ERR=32)
C
      IF (MASS_FRACTION_0.LT.0.) THEN
         YY0(N) = 0. 
         IF (ID.EQ.'OXYGEN') YY0(N) = 0.23
         ELSE
         YY0(N) = MASS_FRACTION_0
         ENDIF
C
      IF (DENSITY.GT.0.) MW = TMPA*DENSITY*R0/PINF
      MU_USER(N) = VISCOSITY
      K_USER(N) = THERMAL_CONDUCTIVITY
      D_USER(N) = DIFFUSION_COEFFICIENT
      MWN(N) = MW
      EPSK(N) = EPSILONKLJ
      SIG(N)  = SIGMALJ
      NUN(N)  = NU
   32 ENDDO SPEC_LOOP
   31 REWIND(LU5)
C
C Set default REACtion parameters
C
      CALL SET_REAC_DEFAULTS
C
C Search for user specified REACtion parameters
C
      REAC_LOOP: DO ITER=1,NFILES
C
      IF (ITER.EQ.2 .AND. REACTION.EQ.'null') EXIT REAC_LOOP
C
      IF (ITER.EQ.1) LUDUM = LU5
      IF (ITER.EQ.2) LUDUM = LU80
      REWIND(LUDUM)
C
      SCANLOOP: DO
      CALL CHECKREAD('REAC',LUDUM,IOS) ; IF (IOS.EQ.1) CYCLE REAC_LOOP
      CALL SET_REAC_DEFAULTS
      READ(LUDUM,REAC,ERR=434,IOSTAT=IOS)
      IF (ID.EQ.REACTION) EXIT REAC_LOOP
  434 IF (IOS.GT.0) THEN
         WRITE(MESSAGE,'(A)') 'ERROR: Problem with REAC '//
     .              TRIM(REACTION)//' or some REAC line preceding it'
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      ENDDO SCANLOOP
C
      ENDDO REAC_LOOP
C
      IF (ID.NE.REACTION .AND. REACTION.NE.'null') THEN
      WRITE(MESSAGE,'(A,A,A)') 
     .    'ERROR: REAC ',TRIM(REACTION),' not found'
      CALL SHUTDOWN(MESSAGE)
      ENDIF
C
      REWIND(LU5)
      IF (DATABASE.NE.'null') REWIND(LU80)
C