read.f90

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

MODULE READ_INPUT
 
USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
USE TRAN
USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY: SECOND,CHECKREAD, SHUTDOWN
USE MEMORY_FUNCTIONS, ONLY: ChkMemErr
USE COMP_FUNCTIONS, ONLY: GET_INPUT_FILE
 
IMPLICIT NONE
PRIVATE
CHARACTER(255), PARAMETER :: readid='$Id: read.f90 719 2007-10-01 17:09:23Z mcgratta $'
CHARACTER(255), PARAMETER :: readrev='$Revision: 719 $'
CHARACTER(255), PARAMETER :: readdate='$Date: 2007-10-01 13:09:23 -0400 (Mon, 01 Oct 2007) $'

PUBLIC READ_DATA, GET_REV_read

CHARACTER(30) :: LABEL,MB
CHARACTER(100) :: MESSAGE,FYI
CHARACTER(30) :: ID,SURF_DEFAULT,BACKGROUND_SPECIES,EVAC_SURF_DEFAULT
LOGICAL :: SUCCESS,EX,THICKEN_OBSTRUCTIONS,BNDF_DEFAULT,BAD
LOGICAL :: WATER_EVAPORATION=.FALSE.,FUEL_EVAPORATION=.FALSE.
REAL(EB) :: XB(6),TEXTURE_ORIGIN(3)
REAL(EB) :: PBX,PBY,PBZ,MW_BACKGROUND,HUMIDITY
REAL(EB) :: MU_USER(0:MAX_SPECIES),K_USER(0:MAX_SPECIES),D_USER(MAX_SPECIES),EPSK(0:MAX_SPECIES),SIG(0:MAX_SPECIES),MW_MIN,MW_MAX
INTEGER  :: I,J,K,IZERO,IOS
TYPE (MESH_TYPE), POINTER :: M
TYPE(OBSTRUCTION_TYPE), POINTER :: OB
TYPE (VENTS_TYPE), POINTER :: VT
TYPE(SPECIES_TYPE), POINTER :: SS,SS0
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(MATERIAL_TYPE), POINTER :: ML
TYPE(REACTION_TYPE), POINTER :: RN
 
 
CONTAINS
 
 
SUBROUTINE READ_DATA(MYID)

INTEGER, INTENT(IN) :: MYID

! Create an array of output QUANTITY names that are included in the various NAMELIST groups
 
CALL FIXED_OUTPUT_QUANTITIES

! Get the name of the input file by reading the command line argument

CALL GET_INPUT_FILE

! If no input file is given, just print out the version number and stop

IF (FN_INPUT(1:1)==' ') THEN
   IF (MYID==0) THEN
      WRITE(LU_ERR,'(/A)') "Fire Dynamics Simulator"
      IF (SERIAL)   WRITE(LU_ERR,'(/A,A,A)') "Version: ",TRIM(VERSION_STRING)," Serial"
      IF (PARALLEL) WRITE(LU_ERR,'(/A,A,A)') "Version: ",TRIM(VERSION_STRING)," Parallel"
      WRITE(LU_ERR,'(A,I4)') "SVN Revision Number: ",SVN_REVISION_NUMBER
      WRITE(LU_ERR,'(A,A)') "Compile Date: ",TRIM(COMPILE_DATE)
      WRITE(LU_ERR,'(/A)')  "Consult FDS Users Guide Chapter, Running FDS, for further instructions."
      WRITE(LU_ERR,'(/A)')  "Hit Enter to Escape..."
   ENDIF
   READ(5,*)
   STOP
ENDIF

! Stop FDS if the input file cannot be found in the current directory

INQUIRE(FILE=FN_INPUT,EXIST=EX)
IF (.NOT.EX) THEN
   IF (MYID==0) WRITE(LU_ERR,'(A,A,A)') "ERROR: The file, ", TRIM(FN_INPUT),", does not exist in the current directory"
   STOP
ENDIF

! Open the input file

OPEN(LU_INPUT,FILE=FN_INPUT,ACTION='READ')

! Read the input file, NAMELIST group by NAMELIST group

CALL READ_DEAD    ! Scan input file looking for old NAMELIST groups, and stop the run if they exist
CALL READ_HEAD
CALL READ_MESH
CALL READ_TRAN
CALL READ_TIME
CALL READ_MISC
CALL READ_RADI
CALL READ_PROP
CALL READ_PART
CALL READ_DEVC
CALL READ_CTRL
CALL READ_TREE
CALL READ_MATL
CALL READ_SURF
CALL READ_OBST
CALL READ_VENT
CALL READ_REAC
CALL READ_SPEC
CALL PROC_SPEC    ! Set up various SPECies constructs
CALL PROC_SURF_1  ! Set up SURFace constructs for species
CALL READ_RAMP    ! Read in all RAMPs, assuming they have all been identified previously
CALL READ_TABL   ! Read in all TABLs, assuming they have all been identified previously
CALL PROC_MATL    ! Set up various MATeriaL constructs
CALL PROC_SURF_2  ! Set up remaining SURFace constructs
CALL READ_DUMP
CALL READ_CLIP
CALL READ_INIT
CALL READ_ZONE
CALL PROC_WALL    ! Set up grid for 1-D heat transfer in solids
CALL PROC_CTRL    ! Set up various ConTRoL constructs
CALL PROC_PROP    ! Set up various PROPerty constructs
CALL PROC_DEVC    ! Set up various DEViCe constructs
CALL READ_PROF
CALL READ_SLCF
CALL READ_ISOF
CALL READ_BNDF

! Close the input file, and never open it again
 
CLOSE (LU_INPUT)

! Set QUANTITY ambient values
CALL SET_QUANTITIES_AMBIENT
 
END SUBROUTINE READ_DATA



SUBROUTINE READ_DEAD
 
! Look for outdated NAMELIST groups and stop the run if any are found
 
REWIND(LU_INPUT)
CALL CHECKREAD('GRID',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: GRID is no longer a valid NAMELIST group. Read User Guide discussion on MESH.')
REWIND(LU_INPUT)
CALL CHECKREAD('HEAT',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: HEAT is no longer a valid NAMELIST group. Read User Guide discussion on PROP and DEVC.')
REWIND(LU_INPUT)
CALL CHECKREAD('PDIM',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: PDIM is no longer a valid NAMELIST group. Read User Guide discussion on MESH.')
REWIND(LU_INPUT)
CALL CHECKREAD('PIPE',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: PIPE is no longer a valid NAMELIST group. Read User Guide discussion on PROP and DEVC.')
REWIND(LU_INPUT)
CALL CHECKREAD('PL3D',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: PL3D is no longer a valid NAMELIST group. Read User Guide discussion on DUMP.')
REWIND(LU_INPUT)
CALL CHECKREAD('SMOD',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: SMOD is no longer a valid NAMELIST group. Read User Guide discussion on DEVC.')
REWIND(LU_INPUT)
CALL CHECKREAD('SPRK',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: SPRK is no longer a valid NAMELIST group. Read User Guide discussion on PROP and DEVC.')
REWIND(LU_INPUT)
CALL CHECKREAD('THCP',LU_INPUT,IOS)
IF (IOS==0) CALL SHUTDOWN('ERROR: THCP is no longer a valid NAMELIST group. Read User Guide discussion on DEVC.')

REWIND(LU_INPUT)
 
END SUBROUTINE READ_DEAD

 
SUBROUTINE READ_HEAD
 
NAMELIST /HEAD/ TITLE,CHID,FYI
 
CHID    = 'output'
TITLE   = '      '
 
REWIND(LU_INPUT)
HEAD_LOOP: DO
   CALL CHECKREAD('HEAD',LU_INPUT,IOS)
   IF (IOS==1) EXIT HEAD_LOOP
   READ(LU_INPUT,HEAD,END=13,ERR=14,IOSTAT=IOS)
   14 IF (IOS>0) CALL SHUTDOWN('ERROR: Problem with HEAD line')
ENDDO HEAD_LOOP
13 REWIND(LU_INPUT)
 
CLOOP: DO I=1,39
   IF (CHID(I:I)=='.') CALL SHUTDOWN('ERROR: No periods allowed in CHID')
   IF (CHID(I:I)==' ') EXIT CLOOP
ENDDO CLOOP
 
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
 
END SUBROUTINE READ_HEAD
 
 
SUBROUTINE READ_MESH

INTEGER :: IJK(3),NM
INTEGER :: IBAR2,JBAR2,KBAR2,POISSON_BC(6),IC,JC,KC,RGB(3)
LOGICAL :: EVACUATION, EVAC_HUMANS
CHARACTER(25) :: COLOR
REAL(EB) :: XB(6)
NAMELIST /MESH/ IJK,FYI,ID,SYNCHRONIZE,EVACUATION,EVAC_HUMANS,POISSON_BC, &
                IBAR2,JBAR2,KBAR2,CYLINDRICAL,XB,RGB,COLOR
TYPE (MESH_TYPE), POINTER :: M
 
NMESHES = 0
 
REWIND(LU_INPUT)
COUNT_MESH_LOOP: DO
   CALL CHECKREAD('MESH',LU_INPUT,IOS)
   IF (IOS==1) EXIT COUNT_MESH_LOOP
   READ(LU_INPUT,MESH,END=15,ERR=16,IOSTAT=IOS)
   NMESHES = NMESHES + 1
   16 IF (IOS>0) CALL SHUTDOWN('ERROR: Problem with MESH line.')
ENDDO COUNT_MESH_LOOP
15 CONTINUE


! Allocate parameters associated with the mesh.
 
ALLOCATE(MESHES(NMESHES),STAT=IZERO)
CALL ChkMemErr('READ','MESHES',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(INTERPOLATED(NMESHES),STAT=IZERO)
CALL ChkMemErr('READ','INTERPOLATED',IZERO)
ALLOCATE(CHANGE_TIME_STEP(NMESHES),STAT=IZERO)
CALL ChkMemErr('READ','CHANGE_TIME_STEP',IZERO)
CHANGE_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)

! Read in the Mesh lines from Input file

REWIND(LU_INPUT)
 
MESH_LOOP: DO NM=1,NMESHES
   IJK(1)=10
   IJK(2)=10
   IJK(3)=10
   IBAR2=1
   JBAR2=1
   KBAR2=1
   TWO_D = .FALSE.
   XB(1) = 0._EB
   XB(2) = 1._EB
   XB(3) = 0._EB
   XB(4) = 1._EB
   XB(5) = 0._EB
   XB(6) = 1._EB
   RGB   = -1
   COLOR = 'null'
   CYLINDRICAL = .FALSE.
   ID = 'null'
   SYNCHRONIZE = .TRUE.
   EVACUATION  = .FALSE.
   EVAC_HUMANS = .FALSE.
   POISSON_BC  = -1
   WRITE(MESH_NAME(NM),'(A,I3)') 'MESH',NM
   CALL CHECKREAD('MESH',LU_INPUT,IOS)
   IF (IOS==1) EXIT MESH_LOOP
   READ(LU_INPUT,MESH)
   M => MESHES(NM)
   M%IBAR = IJK(1)
   M%JBAR = IJK(2)
   M%KBAR = IJK(3)
   M%IBAR2 = IBAR2
   M%JBAR2 = JBAR2
   M%KBAR2 = KBAR2
   M%NEWC = 2*M%IBAR*M%JBAR+2*M%IBAR*M%KBAR+2*M%JBAR*M%KBAR
   IF (SYNCHRONIZE) SYNC_TIME_STEP(NM)  = .TRUE.
   IF (EVACUATION)  EVACUATION_ONLY(NM) = .TRUE.
   IF (EVAC_HUMANS) EVACUATION_GRID(NM) = .TRUE.
   IF (M%JBAR==1) TWO_D = .TRUE.
   IF (TWO_D .AND. M%JBAR/=1) THEN
      WRITE(MESSAGE,'(A)') 'ERROR: IJK(2) must be 1 for all grids in 2D Calculation'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   
   ! Mesh boundary colors
   
   IF (ANY(RGB<0) .AND. COLOR=='null') COLOR = 'BLACK'
   IF (COLOR /= 'null') CALL COLOR2RGB(RGB,COLOR)
   ALLOCATE(M%RGB(3))
   M%RGB = RGB
   
   ! Mesh Geometry and Name
   
   IF (NMESHES > 1 .AND. CYLINDRICAL) THEN
      WRITE(MESSAGE,'(A)') 'ERROR: Cannot use more than 1 MESH when CYLINDRICAL=.TRUE.'
      CALL SHUTDOWN(MESSAGE)
   ENDIF 
   
   IF (ID/='null') MESH_NAME(NM) = ID
   
   ! Kevin's experimental Pressure code
   
   PBC(:,NM) = POISSON_BC(:)
   IF (MOD(M%IBAR,IBAR2)/=0 .OR. MOD(M%JBAR,JBAR2)/=0 .OR. MOD(M%KBAR,KBAR2)/=0) THEN
      WRITE(MESSAGE,'(A)') 'ERROR: IBAR2, JBAR2 or KBAR2 not right'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   ALLOCATE(M%I_LO(IBAR2))
   ALLOCATE(M%I_HI(IBAR2))
   ALLOCATE(M%J_LO(JBAR2))
   ALLOCATE(M%J_HI(JBAR2))
   ALLOCATE(M%K_LO(KBAR2))
   ALLOCATE(M%K_HI(KBAR2))
   DO I=1,IBAR2
      M%I_LO(I) = NINT((I-1)*REAL(M%IBAR,EB)/REAL(IBAR2,EB)) + 1
      M%I_HI(I) = NINT( I   *REAL(M%IBAR,EB)/REAL(IBAR2,EB))
   ENDDO
   DO J=1,JBAR2
      M%J_LO(J) = NINT((J-1)*REAL(M%JBAR,EB)/REAL(JBAR2,EB)) + 1
      M%J_HI(J) = NINT( J   *REAL(M%JBAR,EB)/REAL(JBAR2,EB))
   ENDDO
   DO K=1,KBAR2
      M%K_LO(K) = NINT((K-1)*REAL(M%KBAR,EB)/REAL(KBAR2,EB)) + 1
      M%K_HI(K) = NINT( K   *REAL(M%KBAR,EB)/REAL(KBAR2,EB))
   ENDDO

   ! Process Physical Coordinates
   IF (XB(1) > XB(2)) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: XMIN > XMAX on MESH ', NM
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (XB(3) > XB(4)) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: YMIN > YMAX on MESH ', NM
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (XB(5) > XB(6)) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: ZMIN > ZMAX on MESH ', NM
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   M%XS    = XB(1)
   M%XF    = XB(2)
   M%YS    = XB(3)
   M%YF    = XB(4)
   M%ZS    = XB(5)
   M%ZF    = XB(6)
   M%DXI   = (M%XF-M%XS)/REAL(M%IBAR,EB)
   M%DETA  = (M%YF-M%YS)/REAL(M%JBAR,EB)
   M%DZETA = (M%ZF-M%ZS)/REAL(M%KBAR,EB)
   M%RDXI  = 1._EB/M%DXI
   M%RDETA = 1._EB/M%DETA
   M%RDZETA= 1._EB/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
ENDDO MESH_LOOP
REWIND(LU_INPUT)
 
! Set up coarse grid arrays for Kevin's Pressure Code
 
NCGC = 0
MESH_LOOP_2: DO NM=1,NMESHES
   IF(EVACUATION_ONLY(NM)) CYCLE MESH_LOOP_2
   M=>MESHES(NM)
   ALLOCATE(M%CGI(M%IBAR,M%JBAR,M%KBAR))
   ALLOCATE(M%CGI2(M%IBAR2,M%JBAR2,M%KBAR2))
   DO KC=1,M%KBAR2
      DO JC=1,M%JBAR2
         DO IC=1,M%IBAR2
            NCGC = NCGC+1
            M%CGI2(IC,JC,KC) = NCGC
            M%CGI(M%I_LO(IC):M%I_HI(IC), M%J_LO(JC):M%J_HI(JC), M%K_LO(KC):M%K_HI(KC)) = NCGC
         ENDDO
      ENDDO
   ENDDO
ENDDO MESH_LOOP_2
 
! Start the timing arrays
 
TUSED      = 0._EB
TUSED(1,:) = SECOND()
 
END SUBROUTINE READ_MESH



SUBROUTINE READ_TRAN
USE MATH_FUNCTIONS, ONLY : GAUSSJ
!
! Compute the polynomial transform function for the vertical coordinate
!
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: A,XX
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ND
REAL(EB) :: PC,CC,COEF,XI,ETA,ZETA
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
!