read.f

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

      MODULE READ 
C
      USE PREC
      USE VARS
      USE CONS
      USE TRAN
      USE PACKER
C
      IMPLICIT NONE
C
      PRIVATE
      PUBLIC READ_DATA
      CHARACTER(60) FYI
      CHARACTER(30) QUANTITY,MAKE,LABEL,CB
      CHARACTER(30) COLOR_QUANTITY
      CHARACTER(80) MESSAGE
      CHARACTER(6)  PHASE
      CHARACTER(26) PART_ID,SURF_ID,SURF_IDS(3),SURF_ID6(6),
     .              ID,AKA,AKA_NAME(0:200),
     .              RAMP_MF(0:20),RAMP_Q,RAMP_V,RAMP_G,SURF_DEFAULT,
     .              BACKGROUND_SPECIES,REACTION,RAMP_KS,RAMP_C_P,
     .              RAMP_KS_CHAR,RAMP_C_P_CHAR,EVAC_SURF_DEFAULT
      CHARACTER(20) PROFILE,BACKING,PARTICLE_COLOR
      LOGICAL SUCCESS,BNDF_FACE(-3:3),BNDF_BLOCK,EX,
     .        THICKEN_OBSTRUCTIONS,PARTICLES
      LOGICAL PAPER_MODEL,BNDF_DEFAULT,ADIABATIC,BAD,
     .        BURN_AWAY,LEAKING,OUTLINE
      REAL(EB) KS,TAU_MF(0:20)
      REAL(EB) XYZ(3),XB(6),VEL_T(2),TEXTURE_ORIGIN(3)
      CHARACTER(60) TEXTURE_MAP
      REAL(EB) TM,TAU_Q,TAU_V,TIGN,F,HRRPUA,EMISSIVITY,TEXTURE_WIDTH,
     .         TEXTURE_HEIGHT,
     .         E_COEFFICIENT,RADIATIVE_FRACTION,VOLUME_FLUX,
     .         C_HORIZONTAL,C_VERTICAL,
     .         ALPHA,C_P,TMPWAL,TMPWAL0,DELTA,VEL,C_DELTA_RHO,VBC,
     .         TMPIGN,TMPEVAP,VBC0,DTSAM,PBX,PBY,PBZ,
     .         T,MASS_FLUX(0:20),MASS_FRACTION(0:20),XI,ETA,ZETA,
     .         Z0,PLE,DELTAH,SURFACE_DENSITY,
     .         MW_BACKGROUND,MW,DENSITY,DUMMY,
     .         THERMAL_CONDUCTIVITY,VISCOSITY,HEAT_FLUX
      REAL(EB) DIFFUSION_COEFFICIENT,MU_USER(0:20),K_USER(0:20),
     .         D_USER(0:20),POROSITY,
     .         HEAT_OF_VAPORIZATION,HEAT_OF_COMBUSTION,
     .         HEAT_OF_ABLATION,ABLATION_TEMPERATURE,ABLATION_RATE,
     .         BURNING_RATE_MAX,X_O2_LL,IN_DEPTH_COEFFICIENT,
     .         GEN_TO_YLD,HUMIDITY,MW_MIN,MW_MAX,
     .         ORIENTATION(3),ROTATION,RGB(3),RGB4(4),RADIUS,
     .         MASS_FLUX_CRITICAL,A,
     .         E,MOISTURE_FRACTION,FUEL_FRACTION,
     .         CHAR_DENSITY, C_P_CHAR, KS_CHAR,
     .         CRITICAL_FLAME_TEMPERATURE,DX_SOLID,
     .         EXTERNAL_FLUX,TMP_BACK
      INTEGER NSPRO,NTCO,II,NSFO,NNN,NR,WALL_POINTS,
     .        NPPC,IOR,NSPC,ND,NN,
     .        I1,I2,J1,J2,K1,K2,N,I,J,K,IZERO,
     .        IOS,NVO,LUDUM,ITER,NFILES,NITER
      TYPE (MESH_TYPE), POINTER :: M
      TYPE(OBSTRUCTION_TYPE), POINTER :: OB,OB2
      TYPE (VENTS_TYPE), POINTER :: VT
      TYPE(LAGRANGIAN_TYPE), POINTER :: LP
C
C
      CONTAINS
C
C
      SUBROUTINE READ_DATA
C
      CALL READ_HEAD
      CALL READ_GRID
      CALL READ_PDIM
      CALL READ_TRAN
      CALL READ_TIME
      CALL READ_MISC
      CALL READ_OBST
      CALL READ_VENT
      CALL READ_PART
      CALL READ_TREE
      CALL READ_SPEC
      CALL READ_SURF
      CALL READ_INIT
      CALL READ_RAMP
      CALL READ_SPRK
      CALL READ_HEAT
      CALL READ_SMOD
      CALL READ_THCP
      CALL READ_PL3D
      CALL READ_SLCF
      CALL READ_ISOF
      CALL READ_BNDF
C
      END SUBROUTINE READ_DATA
C
C
C
      SUBROUTINE READ_HEAD
C
      NAMELIST /HEAD/ TITLE,CHID,FYI
C
      CHID    = 'output'
      TITLE   = '      '
C
      HEAD_LOOP: DO
      CALL CHECKREAD('HEAD',LU5,IOS) ; IF (IOS.EQ.1) EXIT HEAD_LOOP
      READ(LU5,HEAD,END=13,ERR=14,IOSTAT=IOS)
   14 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with HEAD line')
      ENDDO HEAD_LOOP
   13 REWIND(LU5)
C
      CLOOP: DO I=1,39
      IF (CHID(I:I).EQ.'.')
     .   CALL SHUTDOWN('ERROR: No periods allowed in CHID')
      IF (CHID(I:I).EQ.' ') EXIT CLOOP
      ENDDO CLOOP
C
      INQUIRE(FILE=TRIM(CHID)//'.stop',EXIST=EX)
      IF (EX) THEN
         WRITE(MESSAGE,'(A,A,A)') "ERROR: Remove the file, ",
     .   TRIM(CHID)//'.stop',", from the current directory"
         CALL SHUTDOWN(MESSAGE)
         ENDIF
C
      END SUBROUTINE READ_HEAD
C
C
      SUBROUTINE READ_GRID
C
      INTEGER :: IBAR,JBAR,KBAR,NM,POISSON_BC(6)
      LOGICAL :: EVACUATION, EVAC_HUMANS
      NAMELIST /GRID/ IBAR,JBAR,KBAR,FYI,ID,SYNCHRONIZE,
     .     EVACUATION,EVAC_HUMANS,POISSON_BC
      TYPE (MESH_TYPE), POINTER :: M
C
      NMESHES = 0
C
      GRID_LOOP: DO
      CALL CHECKREAD('GRID',LU5,IOS) ; IF (IOS.EQ.1) EXIT GRID_LOOP
      READ(LU5,GRID,END=15,ERR=16,IOSTAT=IOS)
      NMESHES = NMESHES + 1
   16 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with GRID line')
      ENDDO GRID_LOOP
   15 REWIND(LU5)
C
      ALLOCATE(MESH(NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','MESH',IZERO)
      ALLOCATE(MESH_NAME(NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','MESH_NAME',IZERO)
      ALLOCATE(TUSED(N_TIMERS,NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','TUSED',IZERO)
      ALLOCATE(SYNC_TIME_STEP(NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','SYNC_TIME_STEP',IZERO)
      SYNC_TIME_STEP = .FALSE.
      ALLOCATE(EVACUATION_ONLY(NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','EVACUATION_ONLY',IZERO)
      EVACUATION_ONLY = .FALSE.
      ALLOCATE(EVACUATION_GRID(NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','EVACUATION_GRID',IZERO)
      EVACUATION_GRID = .FALSE.
      ALLOCATE(PBC(6,NMESHES),STAT=IZERO)
      CALL ChkMemErr('READ','PBC',IZERO)
C
C Read and store GRID dimensions
C
      MESH_LOOP: DO NM=1,NMESHES
      IBAR=10 ; JBAR=10 ; KBAR=10
      TWO_D = .FALSE.
      ID = 'null'
      SYNCHRONIZE = .FALSE.
      EVACUATION  = .FALSE.
      EVAC_HUMANS = .FALSE.
      POISSON_BC  = -1
      WRITE(MESH_NAME(NM),'(A,I3)') 'MESH',NM
      CALL CHECKREAD('GRID',LU5,IOS) ; IF (IOS.EQ.1) EXIT MESH_LOOP
      READ(LU5,GRID,END=115)
      M => MESH(NM)
      M%IBAR = IBAR
      M%JBAR = JBAR
      M%KBAR = KBAR
      M%NEWC = 2*IBAR*JBAR+2*IBAR*KBAR+2*JBAR*KBAR
      IF (SYNCHRONIZE) SYNC_TIME_STEP(NM)  = .TRUE.
      IF (EVACUATION)  EVACUATION_ONLY(NM) = .TRUE.
      IF (EVAC_HUMANS) EVACUATION_GRID(NM) = .TRUE.
      IF (JBAR.EQ.1) TWO_D = .TRUE.
      IF (TWO_D .AND. JBAR.NE.1) THEN
         WRITE(MESSAGE,'(A)') 'ERROR: JBAR must be 1 for all grids'
         CALL SHUTDOWN(MESSAGE)
         ENDIF
      IF (ID.NE.'null') MESH_NAME(NM) = ID
      PBC(:,NM) = POISSON_BC(:)
      ENDDO MESH_LOOP
  115 REWIND(LU5)
C
C Start the timing arrays
C
      TUSED      = 0.
      TUSED(1,:) = SECOND()
C
      END SUBROUTINE READ_GRID
C
C
      SUBROUTINE READ_PDIM
C
      REAL(EB) :: XBAR,YBAR,ZBAR,XBAR0,YBAR0,ZBAR0
      INTEGER :: NM
C
      REAL(EB) :: RBAR0,RBAR
      NAMELIST /PDIM/ XBAR,YBAR,ZBAR,XBAR0,YBAR0,ZBAR0,FYI,RBAR0,RBAR
C
      MESH_LOOP: DO NM=1,NMESHES
C
      RBAR0 = 0. ; RBAR = -1.
      XBAR0 = 0. ; XBAR =  1.
      YBAR0 = 0. ; YBAR =  1.
      ZBAR0 = 0. ; ZBAR =  1.
      CYLINDRICAL = .FALSE.
C
      CALL CHECKREAD('PDIM',LU5,IOS) ; IF (IOS.EQ.1) EXIT MESH_LOOP
      READ(LU5,PDIM,END=19,ERR=20,IOSTAT=IOS) 
   20 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with PDIM line')
C
      M => MESH(NM)
C
      IF (RBAR.GT.0.) THEN
         CYLINDRICAL = .TRUE.
         XBAR        = RBAR
         XBAR0       = RBAR0
         ENDIF
C
      M%XS    = XBAR0
      M%XF    = XBAR
      M%YS    = YBAR0
      M%YF    = YBAR
      M%ZS    = ZBAR0
      M%ZF    = ZBAR
      M%DXI   = (XBAR-XBAR0)/REAL(M%IBAR,EB)
      M%DETA  = (YBAR-YBAR0)/REAL(M%JBAR,EB)
      M%DZETA = (ZBAR-ZBAR0)/REAL(M%KBAR,EB)
      M%RDXI  = 1./M%DXI
      M%RDETA = 1./M%DETA
      M%RDZETA= 1./M%DZETA
      M%IBM1  = M%IBAR-1
      M%JBM1  = M%JBAR-1
      M%KBM1  = M%KBAR-1
      M%IBP1  = M%IBAR+1
      M%JBP1  = M%JBAR+1
      M%KBP1  = M%KBAR+1
C
      ENDDO MESH_LOOP
   19 REWIND(LU5)
C
      END SUBROUTINE READ_PDIM
C
C
      SUBROUTINE READ_TRAN
C
C Compute the polynomial transform function for the vertical coordinate
C
      REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: A,XX
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ND
      REAL(EB)  PC,CC,COEF
      INTEGER  IEXP,IC,IDERIV,N,K,IERROR,IOS,I,MESH_NUMBER,
     .         NIPX,NIPY,NIPZ,NIPXS,NIPYS,NIPZS,NIPXF,NIPYF,NIPZF,NM
      TYPE (MESH_TYPE), POINTER :: M
      TYPE (TRAN_TYPE), POINTER :: T
      NAMELIST /TRNX/ IDERIV,CC,PC,FYI,MESH_NUMBER
      NAMELIST /TRNY/ IDERIV,CC,PC,FYI,MESH_NUMBER
      NAMELIST /TRNZ/ IDERIV,CC,PC,FYI,MESH_NUMBER
C
C Scan the input file, counting the number of NAMELIST entries
C
      ALLOCATE(TRANS(NMESHES))
C
      MESH_LOOP: DO NM=1,NMESHES
      M => MESH(NM)
      T => TRANS(NM)
C
      DO N=1,3
      T%NOC(N) = 0
      TRNLOOP: DO
      IF (N.EQ.1) THEN
         CALL CHECKREAD('TRNX',LU5,IOS) ; IF (IOS.EQ.1) EXIT TRNLOOP
         MESH_NUMBER = 1
         READ(LU5,NML=TRNX,END=17,ERR=18,IOSTAT=IOS)
         IF (MESH_NUMBER.NE.NM) CYCLE TRNLOOP
         ENDIF
      IF (N.EQ.2) THEN
         CALL CHECKREAD('TRNY',LU5,IOS) ; IF (IOS.EQ.1) EXIT TRNLOOP
         MESH_NUMBER = 1
         READ(LU5,NML=TRNY,END=17,ERR=18,IOSTAT=IOS)
         IF (MESH_NUMBER.NE.NM) CYCLE TRNLOOP
         ENDIF
      IF (N.EQ.3) THEN
         CALL CHECKREAD('TRNZ',LU5,IOS) ; IF (IOS.EQ.1) EXIT TRNLOOP
         MESH_NUMBER = 1
         READ(LU5,NML=TRNZ,END=17,ERR=18,IOSTAT=IOS)
         IF (MESH_NUMBER.NE.NM) CYCLE TRNLOOP
         ENDIF
      T%NOC(N) = T%NOC(N) + 1
   18 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with TRN* line')
      ENDDO TRNLOOP
   17 REWIND(LU5)
      ENDDO
C
      T%NOCMAX = MAX(T%NOC(1),T%NOC(2),T%NOC(3))
      ALLOCATE(A(T%NOCMAX+1,T%NOCMAX+1))
      ALLOCATE(XX(T%NOCMAX+1,3))
      ALLOCATE(ND(T%NOCMAX+1,3))
      ALLOCATE(T%C1(0:T%NOCMAX+1,3))
      T%C1               = 0.
      T%C1(1,1:3)        = 1.
      ALLOCATE(T%C2(0:T%NOCMAX+1,3))
      ALLOCATE(T%C3(0:T%NOCMAX+1,3))
      ALLOCATE(T%CCSTORE(T%NOCMAX,3))
      ALLOCATE(T%PCSTORE(T%NOCMAX,3))
      ALLOCATE(T%IDERIVSTORE(T%NOCMAX,3))
C
      T%ITRAN  = 0
C
      DO IC=1,3
      NLOOP:  DO N=1,T%NOC(IC)
      IDERIV = -1
      IF (IC.EQ.1) THEN
         LOOP1: DO
         CALL CHECKREAD('TRNX',LU5,IOS) ; IF (IOS.EQ.1) EXIT NLOOP
         MESH_NUMBER = 1
         READ(LU5,TRNX,END=1,ERR=2)
         IF (MESH_NUMBER.EQ.NM) EXIT LOOP1
         ENDDO LOOP1
         ENDIF
      IF (IC.EQ.2) THEN
         LOOP2: DO
         CALL CHECKREAD('TRNY',LU5,IOS) ; IF (IOS.EQ.1) EXIT NLOOP
         MESH_NUMBER = 1
         READ(LU5,TRNY,END=1,ERR=2)
         IF (MESH_NUMBER.EQ.NM) EXIT LOOP2
         ENDDO LOOP2
         ENDIF
      IF (IC.EQ.3) THEN
         LOOP3: DO
         CALL CHECKREAD('TRNZ',LU5,IOS) ; IF (IOS.EQ.1) EXIT NLOOP
         MESH_NUMBER = 1
         READ(LU5,TRNZ,END=1,ERR=2)
         IF (MESH_NUMBER.EQ.NM) EXIT LOOP3
         ENDDO LOOP3
         ENDIF
      T%CCSTORE(N,IC) = CC
      T%PCSTORE(N,IC) = PC
      T%IDERIVSTORE(N,IC) = IDERIV
      IF (IDERIV.GE.0) T%ITRAN(IC) = 1
      IF (IDERIV.LT.0) T%ITRAN(IC) = 2
    2 ENDDO NLOOP
    1 REWIND(LU5)
      ENDDO 
C
      ICLOOP: DO IC=1,3
C
      IF (T%ITRAN(IC).EQ.1) THEN
C
C If ITRAN=1, do a polynomial transformation 
C
      ND(1,IC)  = 0
      SELECT CASE(IC)
      CASE(1)
         XX(1,IC)    = M%XF-M%XS
         T%C1(1,IC)  = M%XF-M%XS
      CASE(2)
         XX(1,IC)    = M%YF-M%YS
         T%C1(1,IC)  = M%YF-M%YS
      CASE(3)
         XX(1,IC)    = M%ZF-M%ZS
         T%C1(1,IC)  = M%ZF-M%ZS
      END SELECT
 
      NNLOOP:  DO N=2,T%NOC(IC)+1
      IDERIV = T%IDERIVSTORE(N-1,IC)
      IF (IC.EQ.1) CC = T%CCSTORE(N-1,IC)-M%XS
      IF (IC.EQ.2) CC = T%CCSTORE(N-1,IC)-M%YS
      IF (IC.EQ.3) CC = T%CCSTORE(N-1,IC)-M%ZS
      IF (IC.EQ.1 .AND. IDERIV.EQ.0) PC = T%PCSTORE(N-1,IC)-M%XS
      IF (IC.EQ.2 .AND. IDERIV.EQ.0) PC = T%PCSTORE(N-1,IC)-M%YS
      IF (IC.EQ.3 .AND. IDERIV.EQ.0) PC = T%PCSTORE(N-1,IC)-M%ZS
      IF (IC.EQ.1 .AND. IDERIV.GT.0) PC = T%PCSTORE(N-1,IC)
      IF (IC.EQ.2 .AND. IDERIV.GT.0) PC = T%PCSTORE(N-1,IC)
      IF (IC.EQ.3 .AND. IDERIV.GT.0) PC = T%PCSTORE(N-1,IC)
        ND(N,IC) = IDERIV
        XX(N,IC) = CC
      T%C1(N,IC) = PC
      ENDDO NNLOOP
C
      DO K=1,T%NOC(IC)+1
      DO N=1,T%NOC(IC)+1
      COEF = IFAC(K,ND(N,IC))
      IEXP = K-ND(N,IC)
      IF (IEXP.LT.0) A(N,K) = 0.
      IF (IEXP.EQ.0) A(N,K) = COEF
      IF (IEXP.GT.0) A(N,K) = COEF*XX(N,IC)**IEXP
      ENDDO
      ENDDO
C
      IERROR = 0
      CALL GAUSSJ(A,T%NOC(IC)+1,T%NOCMAX+1,T%C1(1:T%NOCMAX+1,IC),
     .            1,1,IERROR)
      IF (IERROR.NE.0)
     .   CALL SHUTDOWN('ERROR: Problem with grid transformation')
C
      ENDIF
C
      IF (T%ITRAN(IC).EQ.2) THEN
C
C If ITRAN=2, do a linear transformation
C
      T%C1(0,IC) = 0.
      T%C2(0,IC) = 0.
      DO N=1,T%NOC(IC)
      IF (IC.EQ.1) CC = T%CCSTORE(N,IC)-M%XS
      IF (IC.EQ.2) CC = T%CCSTORE(N,IC)-M%YS
      IF (IC.EQ.3) CC = T%CCSTORE(N,IC)-M%ZS
      IF (IC.EQ.1) PC = T%PCSTORE(N,IC)-M%XS
      IF (IC.EQ.2) PC = T%PCSTORE(N,IC)-M%YS
      IF (IC.EQ.3) PC = T%PCSTORE(N,IC)-M%ZS
      T%C1(N,IC) = CC
      T%C2(N,IC) = PC
      ENDDO
C
      SELECT CASE(IC)
      CASE(1)
         T%C1(T%NOC(1)+1,1) = M%XF-M%XS