read.f

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

         T%C2(T%NOC(1)+1,1) = M%XF-M%XS
      CASE(2)
         T%C1(T%NOC(2)+1,2) = M%YF-M%YS
         T%C2(T%NOC(2)+1,2) = M%YF-M%YS
      CASE(3)
         T%C1(T%NOC(3)+1,3) = M%ZF-M%ZS
         T%C2(T%NOC(3)+1,3) = M%ZF-M%ZS
      END SELECT
C
      DO N=1,T%NOC(IC)+1
      T%C3(N,IC) = (T%C2(N,IC)-T%C2(N-1,IC))/(T%C1(N,IC)-T%C1(N-1,IC))
      ENDDO
C
      ENDIF 
C
      ENDDO ICLOOP
C
      DEALLOCATE(A)
      DEALLOCATE(XX)
      DEALLOCATE(ND)
C
C Set up grid stretching arrays
C
      ALLOCATE(M%R(0:M%IBAR),STAT=IZERO)
      CALL ChkMemErr('READ','R',IZERO)
      ALLOCATE(M%RC(0:M%IBAR+1),STAT=IZERO)
      CALL ChkMemErr('READ','RC',IZERO) ; M%RC = 1.
      ALLOCATE(M%RRN(0:M%IBP1),STAT=IZERO)
      CALL ChkMemErr('READ','RRN',IZERO) ; M%RRN = 1.
      ALLOCATE(M%X(0:M%IBAR),STAT=IZERO)
      CALL ChkMemErr('READ','X',IZERO)
      ALLOCATE(M%XC(0:M%IBP1),STAT=IZERO)
      CALL ChkMemErr('READ','XC',IZERO)
      ALLOCATE(M%HX(0:M%IBP1),STAT=IZERO)
      CALL ChkMemErr('READ','HX',IZERO)
      ALLOCATE(M%DX(0:M%IBP1),STAT=IZERO)
      CALL ChkMemErr('READ','DX',IZERO)
      ALLOCATE(M%RDX(0:M%IBP1),STAT=IZERO)
      CALL ChkMemErr('READ','RDX',IZERO)
      ALLOCATE(M%DXN(0:M%IBAR),STAT=IZERO)
      CALL ChkMemErr('READ','DXN',IZERO)
      ALLOCATE(M%RDXN(0:M%IBAR),STAT=IZERO)
      CALL ChkMemErr('READ','RDXN',IZERO)
      ALLOCATE(M%Y(0:M%JBAR),STAT=IZERO)
      CALL ChkMemErr('READ','Y',IZERO)
      ALLOCATE(M%YC(0:M%JBP1),STAT=IZERO)
      CALL ChkMemErr('READ','YC',IZERO)
      ALLOCATE(M%HY(0:M%JBP1),STAT=IZERO)
      CALL ChkMemErr('READ','HY',IZERO)
      ALLOCATE(M%DY(0:M%JBP1),STAT=IZERO)
      CALL ChkMemErr('READ','DY',IZERO)
      ALLOCATE(M%RDY(0:M%JBP1),STAT=IZERO)
      CALL ChkMemErr('READ','RDY',IZERO)
      ALLOCATE(M%DYN(0:M%JBAR),STAT=IZERO)
      CALL ChkMemErr('READ','DYN',IZERO)
      ALLOCATE(M%RDYN(0:M%JBAR),STAT=IZERO)
      CALL ChkMemErr('READ','RDYN',IZERO)
      ALLOCATE(M%Z(0:M%KBAR),STAT=IZERO)
      CALL ChkMemErr('READ','Z',IZERO)
      ALLOCATE(M%ZC(0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('READ','ZC',IZERO)
      ALLOCATE(M%HZ(0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('READ','HZ',IZERO)
      ALLOCATE(M%DZ(0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('READ','DZ',IZERO)
      ALLOCATE(M%RDZ(0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('READ','RDZ',IZERO)
      ALLOCATE(M%DZN(0:M%KBAR),STAT=IZERO)
      CALL ChkMemErr('READ','DZN',IZERO)
      ALLOCATE(M%RDZN(0:M%KBAR),STAT=IZERO)
      CALL ChkMemErr('READ','RDZN',IZERO)
C
C Define X grid stretching terms
C
      M%DXMIN = 1000.
      DO I=1,M%IBAR
      XI    = (REAL(I,EB)-.5)*M%DXI
      M%HX(I) = GP(XI,1,NM)
      M%DX(I) = M%HX(I)*M%DXI
      M%DXMIN = MIN(M%DXMIN,M%DX(I))
      IF (M%HX(I).LE.0.) THEN
         WRITE(MESSAGE,'(A,I2)') 
     .       'ERROR: x transformation not monotonic, mesh ',NM
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      M%RDX(I) = 1./M%DX(I)
      ENDDO
C
      M%HX(0)    = M%HX(1)
      M%HX(M%IBP1) = M%HX(M%IBAR)
      M%DX(0)    = M%DX(1)
      M%DX(M%IBP1) = M%DX(M%IBAR)
      M%RDX(0)    = 1./M%DX(1)
      M%RDX(M%IBP1) = 1./M%DX(M%IBAR)
C
      DO I=0,M%IBAR
      XI     = I*M%DXI
      M%X(I) = M%XS + G(XI,1,NM)
      IF (CYLINDRICAL) THEN ; M%R(I) = M%X(I)
                       ELSE ; M%R(I) = 1. ; ENDIF
      M%DXN(I)  = 0.5*(M%DX(I)+M%DX(I+1))
      M%RDXN(I) = 1./M%DXN(I)
      ENDDO
      M%X(0)      = M%XS
      M%X(M%IBAR) = M%XF
C
      DO I=1,M%IBAR
      M%XC(I) = 0.5*(M%X(I)+M%X(I-1))
      ENDDO
      M%XC(0)      = M%XS - 0.5*M%DX(0)
      M%XC(M%IBP1) = M%XF + 0.5*M%DX(M%IBP1)
C
      IF (CYLINDRICAL) THEN  
      DO I=1,M%IBAR
      M%RRN(I) = 2./(M%R(I)+M%R(I-1))
      M%RC(I)  = 0.5*(M%R(I)+M%R(I-1))
      ENDDO
      M%RRN(0)    = M%RRN(1)
      M%RRN(M%IBP1) = M%RRN(M%IBAR)
      ENDIF
C
C Define Y grid stretching terms
C
      M%DYMIN = 1000.
      DO J=1,M%JBAR
      ETA   = (REAL(J,EB)-.5)*M%DETA
      M%HY(J) = GP(ETA,2,NM)
      M%DY(J) = M%HY(J)*M%DETA
      M%DYMIN = MIN(M%DYMIN,M%DY(J))
      IF (M%HY(J).LE.0.) THEN
         WRITE(MESSAGE,'(A,I2)') 
     .       'ERROR: y transformation not monotonic, mesh ',NM
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      M%RDY(J) = 1./M%DY(J)
      ENDDO
C
      M%HY(0)    = M%HY(1)
      M%HY(M%JBP1) = M%HY(M%JBAR)
      M%DY(0)    = M%DY(1)
      M%DY(M%JBP1) = M%DY(M%JBAR)
      M%RDY(0)    = 1./M%DY(1)
      M%RDY(M%JBP1) = 1./M%DY(M%JBAR)
C
      DO J=0,M%JBAR
      ETA     = J*M%DETA
      M%Y(J)    = M%YS + G(ETA,2,NM)
      M%DYN(J)  = 0.5*(M%DY(J)+M%DY(J+1))
      M%RDYN(J) = 1./M%DYN(J)
      ENDDO
C
      M%Y(0)      = M%YS
      M%Y(M%JBAR) = M%YF
C
      DO J=1,M%JBAR
      M%YC(J) = 0.5*(M%Y(J)+M%Y(J-1))
      ENDDO
      M%YC(0)      = M%YS - 0.5*M%DY(0)
      M%YC(M%JBP1) = M%YF + 0.5*M%DY(M%JBP1)
C
C Define Z grid stretching terms
C
      M%DZMIN = 1000.
      DO K=1,M%KBAR
      ZETA  = (REAL(K,EB)-.5)*M%DZETA
      M%HZ(K) = GP(ZETA,3,NM)
      M%DZ(K) = M%HZ(K)*M%DZETA
      M%DZMIN = MIN(M%DZMIN,M%DZ(K))
      IF (M%HZ(K).LE.0.) THEN
         WRITE(MESSAGE,'(A,I2)') 
     .       'ERROR: z transformation not monotonic, mesh ',NM
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      M%RDZ(K) = 1./M%DZ(K)
      ENDDO
C
      M%HZ(0)    = M%HZ(1)
      M%HZ(M%KBP1) = M%HZ(M%KBAR)
      M%DZ(0)    = M%DZ(1)
      M%DZ(M%KBP1) = M%DZ(M%KBAR)
      M%RDZ(0)    = 1./M%DZ(1)
      M%RDZ(M%KBP1) = 1./M%DZ(M%KBAR)
C
      DO K=0,M%KBAR
      ZETA      = K*M%DZETA
      M%Z(K)    = M%ZS + G(ZETA,3,NM)
      M%DZN(K)  = 0.5*(M%DZ(K)+M%DZ(K+1))
      M%RDZN(K) = 1./M%DZN(K)
      ENDDO
C
      M%Z(0)      = M%ZS
      M%Z(M%KBAR) = M%ZF
C
      DO K=1,M%KBAR
      M%ZC(K) = 0.5*(M%Z(K)+M%Z(K-1))
      ENDDO
      M%ZC(0)      = M%ZS - 0.5*M%DZ(0)
      M%ZC(M%KBP1) = M%ZF + 0.5*M%DZ(M%KBP1)
C
C Set up arrays that will return coordinate positions
C
      NIPX   = 100*M%IBAR
      NIPY   = 100*M%JBAR
      NIPZ   = 100*M%KBAR
      NIPXS  = NINT(NIPX*M%DX(0)/(M%XF-M%XS))
      NIPXF  = NINT(NIPX*M%DX(M%IBP1)/(M%XF-M%XS))
      NIPYS  = NINT(NIPY*M%DY(0)/(M%YF-M%YS))
      NIPYF  = NINT(NIPY*M%DY(M%JBP1)/(M%YF-M%YS))
      NIPZS  = NINT(NIPZ*M%DZ(0)/(M%ZF-M%ZS))
      NIPZF  = NINT(NIPZ*M%DZ(M%KBP1)/(M%ZF-M%ZS))
      M%RDXINT = REAL(NIPX,EB)/(M%XF-M%XS)
      M%RDYINT = REAL(NIPY,EB)/(M%YF-M%YS)
      M%RDZINT = REAL(NIPZ,EB)/(M%ZF-M%ZS)
C
      ALLOCATE(M%CELLSI(-NIPXS:NIPX+NIPXF),STAT=IZERO)
      CALL ChkMemErr('READ','CELLSI',IZERO)
      ALLOCATE(M%CELLSJ(-NIPYS:NIPY+NIPYF),STAT=IZERO)
      CALL ChkMemErr('READ','CELLSJ',IZERO)
      ALLOCATE(M%CELLSK(-NIPZS:NIPZ+NIPZF),STAT=IZERO)
      CALL ChkMemErr('READ','CELLSK',IZERO)
C
      DO I=-NIPXS,NIPX+NIPXF
      M%CELLSI(I) = GINV(REAL(I,EB)/M%RDXINT,1,NM)*M%RDXI
      M%CELLSI(I) = MAX(M%CELLSI(I),-0.9_EB)
      M%CELLSI(I) = MIN(M%CELLSI(I),REAL(M%IBAR)+0.9_EB)
      ENDDO
      DO J=-NIPYS,NIPY+NIPYF
      M%CELLSJ(J) = GINV(REAL(J,EB)/M%RDYINT,2,NM)*M%RDETA
      M%CELLSJ(J) = MAX(M%CELLSJ(J),-0.9_EB)
      M%CELLSJ(J) = MIN(M%CELLSJ(J),REAL(M%JBAR)+0.9_EB)
      ENDDO
      DO K=-NIPZS,NIPZ+NIPZF
      M%CELLSK(K) = GINV(REAL(K,EB)/M%RDZINT,3,NM)*M%RDZETA
      M%CELLSK(K) = MAX(M%CELLSK(K),-0.9_EB)
      M%CELLSK(K) = MIN(M%CELLSK(K),REAL(M%KBAR)+0.9_EB)
      ENDDO
C
      ENDDO MESH_LOOP
C
C
      CONTAINS
C
      INTEGER FUNCTION IFAC(II,N)
      INTEGER II,N
      IFAC = 1
      DO I=II-N+1,II
      IFAC = IFAC*I
      ENDDO
      END FUNCTION IFAC

      END SUBROUTINE READ_TRAN
C
C
      SUBROUTINE READ_TIME
C
      REAL(EB) :: DT,VEL_CHAR
      INTEGER :: NM
      NAMELIST /TIME/ DT,TWFIN,FYI,WALL_INCREMENT,SYNCHRONIZE,
     .                EVAC_DT_FLOWFIELD,EVAC_DT_STEADY_STATE
      TYPE (MESH_TYPE), POINTER :: M
C
      DT             =-1.
      TWFIN          = 1.0
      WALL_INCREMENT = 2
      SET_UP         = .FALSE.
      SYNCHRONIZE    = .FALSE.
      EVAC_DT_FLOWFIELD = 0.01
      EVAC_DT_STEADY_STATE = 0.05
C
      TIME_LOOP: DO
      CALL CHECKREAD('TIME',LU5,IOS) ; IF (IOS.EQ.1) EXIT TIME_LOOP
      READ(LU5,TIME,END=21,ERR=22,IOSTAT=IOS)
   22 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with TIME line')
      ENDDO TIME_LOOP
   21 REWIND(LU5)
C
      IF (TWFIN.LE.0.) SET_UP = .TRUE.
C
      IF (SYNCHRONIZE) SYNC_TIME_STEP = .TRUE.
      IF (ANY(SYNC_TIME_STEP)) SYNCHRONIZE = .TRUE.
C
      MESH_LOOP: DO NM=1,NMESHES
      M=>MESH(NM)
      IF (DT.GT.0.) THEN
         M%DT = DT
      ELSE
         VEL_CHAR = 0.2*SQRT(10.*(M%ZF-M%ZS))
         M%DT     = (M%DXMIN*M%DYMIN*M%DZMIN)**(1./3.)/VEL_CHAR
      ENDIF
      IF (EVACUATION_ONLY(NM)) THEN
         SYNC_TIME_STEP(NM) = .FALSE.
         M%DT = EVAC_DT_FLOWFIELD
         ENDIF
      ENDDO MESH_LOOP
C
      END SUBROUTINE READ_TIME
C
C
      SUBROUTINE READ_MISC
C
      REAL(EB) MAX_OVER_PRESSURE,DTSAM_PART,GAUGE_TEMPERATURE
      NAMELIST /MISC/ PR,SC,TMPA,TMPO,GVEC,RF,NFRAMES,FYI,
     .                CSMAG,RAMP_G,BAROCLINIC,
     .                DT0DZ,ISOTHERMAL,INCOMPRESSIBLE,
     .                PINF,DATABASE,SURF_DEFAULT,EVAC_SURF_DEFAULT,
     .                DENSITY,DATABASE_DIRECTORY,RENDER_FILE,
     .                C_FORCED,C_VERTICAL,C_HORIZONTAL,RESTART,
     .                DTCORE,BACKGROUND_SPECIES,MW,LES,DNS,
     .                VISCOSITY,THERMAL_CONDUCTIVITY,NOISE,
     .                RADIATION,GAMMA,BNDF_DEFAULT,REACTION,
     .                MAX_OVER_PRESSURE,AUTOMATIC_Z,U0,V0,W0,HUMIDITY,
     .                POROUS_FLOOR,SUPPRESSION,CHECK_POISSON,
     .                TEXTURE_ORIGIN,NSTRATA,SMOKE3D,
     .                THICKEN_OBSTRUCTIONS,LEAK_AREA,
     .                DROP_VERTICAL_VELOCITY,
     .                DROP_HORIZONTAL_VELOCITY,RHO_SOOT,
     .                DTSAM_PART,NPPS,DTPAR,DTSPAR,DEBUG,TIMING,
     .                GAUGE_TEMPERATURE,CHARACTERISTIC_VELOCITY,
     .                FLUSH_FILE_BUFFERS,MAXIMUM_DROPLETS,
     .                EVAC_PRESSURE_ITERATIONS,EVAC_TIME_ITERATIONS,
     .                TEXTURE_DIRECTORY
C
C Physical constants
C
      R0      = 8314.3    ! Universal Gas Constant (J/kmol/K)
      R1      = 1.9862E-3 ! Universal Gas Constant (kcal/mol/K)
      TMPA    = 20.       ! Ambient temperature (C)
      TMPO    = -10000. 
      GRAV    = 9.81      ! Acceleration of gravity (m/s**2)
      GAMMA   = 1.4       ! Heat capacity ratio for air
      PINF    = 101325.   ! Ambient pressure (Pa)
      TMPM    = 273.15    ! Melting temperature of water (K)
      SIGMA   = 5.67E-8   ! Stefan-Boltzmann constant (W/m**2/K**4)
      HUMIDITY= -1.       ! Relative Humidity
      RHO_SOOT= 1850.     ! Density of soot particle (kg/m3)
C
C Empirical constants
C
      C_VERTICAL   = 1.31 ! Vertical free convection (Holman, Table 7-2)
      C_HORIZONTAL = 1.52 ! Horizontal free convection 
      C_FORCED     = 0.037 ! Forced convection coefficient
C
C Miscellanious parameters
C
      PI      = 4.*ATAN(1.0_EB)
      RPI     = 1./PI
      TWOPI   = 2.*PI
      PIO2    = PI/2.
      ONTH    = 1./3.
      THFO    = 3./4.
      ONSI    = 1./6.
      TWTH    = 2./3.
      FOTH    = 4./3.
      RFPI    = 1./(4.*PI)
C
C Numerical Parameters
C
      U0 = 0.
      V0 = 0.
      W0 = 0.
      BACKGROUND_SPECIES = 'AIR'
      VISCOSITY = -1.
      THERMAL_CONDUCTIVITY = -1.
      MU_USER = -1.
      K_USER  = -1.
      MW      = 0.   
      DENSITY = -1.
      RESTART = .FALSE.
      RADIATION      = .TRUE.
      SUPPRESSION    = .TRUE.
      CHECK_POISSON  = .FALSE.
      BAROCLINIC     = .FALSE.
      NOISE          = .TRUE.
      ISOTHERMAL     = .FALSE.  
      INCOMPRESSIBLE = .FALSE.
      BNDF_DEFAULT   = .TRUE.
      LES            = .TRUE.
      DNS            = .FALSE.
      COMBUSTION_MODEL = 1
      AUTOMATIC_Z    = .TRUE.
      MAX_OVER_PRESSURE = 1000000.
      POROUS_FLOOR = .TRUE.
      TEXTURE_ORIGIN(1) = 0.
      TEXTURE_ORIGIN(2) = 0.
      TEXTURE_ORIGIN(3) = 0.
      LEAK_AREA = 0.
      GAUGE_TEMPERATURE = -1.
      CHARACTERISTIC_VELOCITY = 1.
      MAXIMUM_DROPLETS = 500000
C
C EVACuation parameters
C
      EVAC_PRESSURE_ITERATIONS = 50
      EVAC_TIME_ITERATIONS     = 50
C
C LES parameters
C
      CSMAG                = 0.20     ! Smagorinsky constant
      PR                   = -1.0     ! Turbulent Prandtl number
      SC                   = -1.0     ! Turbulent Schmidt number
C
C Misc
C
      DATABASE             = 'null'
      DATABASE_DIRECTORY   = 'null'
      TEXTURE_DIRECTORY    = 'null'
      RAMP_G               = 'null'
      SURF_DEFAULT         = 'INERT'
      EVAC_SURF_DEFAULT    = 'INERT'
      RENDER_FILE          = 'null'
      REACTION             = 'null'
      DTCORE               = 1000000.
      GVEC(1)              = 0.        ! x-component of gravity 
      GVEC(2)              = 0.        ! y-component of gravity 
      GVEC(3)              = -GRAV     ! z-component of gravity 
      DT0DZ                = -1000000.
      NFRAMES              = 1000      ! Number of output data sets
      DTSAM_PART           = -1.
      NPPS                 = 100000    ! Number Particles Per Set
      DTPAR                = 0.05      ! Particle Insertion Interval
      DTSPAR               = 0.05      ! Droplet Insertion Interval
      RF                   = 1.00      ! Relaxation factor for no-flux
      DROPLET_FILE         = .FALSE.
      NSTRATA              = 7         ! Number bins for drop dist.
      IF (.NOT.TWO_D) SMOKE3D = .TRUE.
      IF (     TWO_D) SMOKE3D = .FALSE.