init.f90

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

ALLOCATE(M%RHO_W(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RHO_W',IZERO) 
M%RHO_W = RHOA
ALLOCATE(M%RSUM_W(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RSUM_W',IZERO) 
M%RSUM_W = RSUM0
ALLOCATE(M%YY_W(NDWC,N_SPECIES),STAT=IZERO)
CALL ChkMemErr('INIT','YY_W',IZERO)  
DO IW=1,NDWC
   M%YY_W(IW,1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
ENDDO
ALLOCATE(M%E_WALL(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','E_WALL',IZERO)
M%E_WALL = 1._EB
ALLOCATE(M%QRADIN(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','QRADIN',IZERO) 
M%QRADIN = SIGMA*TMPA4
ALLOCATE(M%QRADOUT(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','QRADOUT',IZERO) 
M%QRADOUT = SIGMA*TMPA4
ALLOCATE(M%QCONF(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','QCONF',IZERO) 
M%QCONF = 0._EB
ALLOCATE(M%HEAT_TRANS_COEF(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','HEAT_TRANS_COEF',IZERO) 
M%HEAT_TRANS_COEF = 0._EB
ALLOCATE(M%XW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','XW',IZERO)
ALLOCATE(M%YW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','YW',IZERO)
ALLOCATE(M%ZW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','ZW',IZERO)
ALLOCATE(M%TW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','TW',IZERO)  
M%TW = T_BEGIN
ALLOCATE(M%EW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','EW',IZERO)  
M%EW = 0._EB
ALLOCATE(M%KW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','KW',IZERO)  
M%KW = 1._EB
ALLOCATE(M%RHODW(NDWC,N_SPECIES),STAT=IZERO)
CALL ChkMemErr('INIT','RHODW',IZERO) 
M%RHODW = 1._EB
ALLOCATE(M%AREA_ADJUST(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','AREA_ADJUST',IZERO)  
M%AREA_ADJUST = 1._EB
ALLOCATE(M%MASSFLUX(NDWC,0:N_SPECIES),STAT=IZERO)
CALL ChkMemErr('INIT','MASSFLUX',IZERO)  
M%MASSFLUX = 0._EB
ALLOCATE(M%ACTUAL_BURN_RATE(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','ACTUAL_BURN_RATE',IZERO)  
M%ACTUAL_BURN_RATE = 0._EB
ALLOCATE(M%RDN(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RDN',IZERO)
M%RDN = 1._EB
ALLOCATE(M%AW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','AW',IZERO)  
M%AW = 0._EB
ALLOCATE(M%RAW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RAW',IZERO) 
M%RAW = 0._EB
ALLOCATE(M%NPPCW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','NPPCW',IZERO) 
M%NPPCW = 1
ALLOCATE(M%UW0(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','UW0',IZERO)    
M%UW0 = 0._EB
ALLOCATE(M%UW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','UW',IZERO)    
M%UW = 0._EB
ALLOCATE(M%UWS(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','UWS',IZERO)   
M%UWS = 0._EB
ALLOCATE(M%OBST_INDEX_W(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','OBST_INDEX_W',IZERO) 
M%OBST_INDEX_W = 0
ALLOCATE(M%VENT_INDEX(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','VENT_INDEX',IZERO) 
M%VENT_INDEX = 0
ALLOCATE(M%DUWDT(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','DUWDT',IZERO) 
M%DUWDT = 0._EB
ALLOCATE(M%PRESSURE_ZONE_WALL(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ZONE_WALL',IZERO) 
M%PRESSURE_ZONE_WALL = 0
ALLOCATE(M%IJKW(12,NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','IJKW',IZERO)  
M%IJKW = 0
ALLOCATE(M%INTERPOLATION_FACTOR(M%NEWC,3),STAT=IZERO)
CALL ChkMemErr('INIT','INTERPOLATION_FACTOR',IZERO)  
ALLOCATE(M%CELL_VOLUME_RATIO(M%NEWC),STAT=IZERO)
CALL ChkMemErr('INIT','CELL_VOLUME_RATIO',IZERO)  
ALLOCATE(M%BOUNDARY_TYPE(0:NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','BOUNDARY_TYPE',IZERO) 
M%BOUNDARY_TYPE = NULL_BOUNDARY
ALLOCATE(M%MASS_LOSS(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','MASS_LOSS',IZERO)
M%MASS_LOSS = 0._EB
ALLOCATE(M%WALL_INDEX(0:M%CELL_COUNT,-3:3),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_INDEX',IZERO) 
M%WALL_INDEX = 0
ALLOCATE(M%WALL_INDEX_BACK(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_INDEX_BACK',IZERO) 
M%WALL_INDEX_BACK = 0
ALLOCATE(M%EDGE_INDEX(0:M%CELL_COUNT,1:12),STAT=IZERO)
CALL ChkMemErr('INIT','EDGE_INDEX',IZERO) 
M%EDGE_INDEX = 0
ALLOCATE(M%UVW_GHOST(0:M%CELL_COUNT,3),STAT=IZERO)
CALL ChkMemErr('INIT','UVW_GHOST',IZERO) 
M%UVW_GHOST = 0
 
! Surface water arrays
 
IF (ACCUMULATE_WATER) THEN
ALLOCATE(M%AWMPUA(NDWC,N_EVAP_INDICIES),STAT=IZERO)
CALL ChkMemErr('INIT','AWMPUA',IZERO) 
M%AWMPUA = 0._EB
ENDIF
ALLOCATE(M%WMPUA(NDWC,MAX(1,N_EVAP_INDICIES)),STAT=IZERO)
CALL ChkMemErr('INIT','WMPUA',IZERO) 
M%WMPUA = 0._EB
ALLOCATE(M%WCPUA(NDWC,MAX(1,N_EVAP_INDICIES)),STAT=IZERO)
CALL ChkMemErr('INIT','WCPUA',IZERO) 
M%WCPUA = 0._EB
 
! Surface work arrays
 
ALLOCATE(M%WALL_WORK1(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK1',IZERO) 
ALLOCATE(M%WALL_WORK2(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK2',IZERO) 
 
! Set up boundary arrays for external boundaries of the current mesh

DO K=1,KBAR
   DO J=1,JBAR
      I   = 0
      IBC = DEFAULT_SURF_INDEX
      IOR = 1
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC)
   ENDDO
ENDDO
DO K=1,KBAR
   DO J=1,JBAR
      I   = IBP1
      IBC = DEFAULT_SURF_INDEX
      IOR = -1
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC)
   ENDDO
ENDDO
 
DO K=1,KBAR
   DO I=1,IBAR
      J   = 0
      IBC = DEFAULT_SURF_INDEX
      IOR = 2
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC)
   ENDDO
ENDDO
DO K=1,KBAR
   DO I=1,IBAR
      J   = JBP1
      IBC = DEFAULT_SURF_INDEX
      IOR = -2
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC)
   ENDDO
ENDDO
 
DO J=1,JBAR
   DO I=1,IBAR
      K   = 0
      IBC = DEFAULT_SURF_INDEX
      IOR = 3
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC)
   ENDDO
ENDDO
DO J=1,JBAR
   DO I=1,IBAR
      K   = KBP1
      IBC = DEFAULT_SURF_INDEX
      IOR = -3
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC)
   ENDDO
ENDDO

! Go through all obstructions and decide which cell faces ought to be given a wall cell index and initialized
 
OBST_LOOP_2: DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
 
   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I1+1 
         IF (I==1) CYCLE   ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I-1,J,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = -1 
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC)
      ENDDO
   ENDDO
 
   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I2
         IF (I==M%IBAR) CYCLE  ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I+1,J,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = 1
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC)
      ENDDO 
   ENDDO
 
   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J1+1
         IF (J==1) CYCLE   ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J-1,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = -2
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC)
      ENDDO
   ENDDO   
 
   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J2
         IF (J==M%JBAR) CYCLE  ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J+1,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = 2
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC)
      ENDDO
   ENDDO   
 
   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K1+1
         IF (K==1) CYCLE   ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J,K-1)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = -3
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC)
      ENDDO
   ENDDO   
 
   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K2
         IF (K==M%KBAR) CYCLE  ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J,K+1)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = 3 
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC)
      ENDDO
   ENDDO   
 
ENDDO OBST_LOOP_2

! Determine back wall index for exposed surfaces

DO IW=M%NEWC+1,M%NWC  
   ! Only assign WALL_INDEX_BACK to wall cells that are not attached to the exterior boundary of the computational domain
   SF=>SURFACE(M%IJKW(5,IW))
   IF (SF%BACKING==EXPOSED) THEN
      II = M%IJKW(1,IW)
      JJ = M%IJKW(2,IW)
      KK = M%IJKW(3,IW)
      IC = M%CELL_INDEX(II,JJ,KK)
      IOR = M%IJKW(4,IW)
      IF (.NOT.M%SOLID(IC)) M%WALL_INDEX_BACK(IW) = M%WALL_INDEX(IC,IOR)
      IF (     M%SOLID(IC)) THEN
         SELECT CASE(IOR)
            CASE(-1)
               II=II+1
            CASE( 1)
               II=II-1
            CASE(-2)
               JJ=JJ+1
            CASE( 2)
               JJ=JJ-1
            CASE(-3)
               KK=KK+1
            CASE( 3)
               KK=KK-1
         END SELECT
         IC = M%CELL_INDEX(II,JJ,KK)
         M%WALL_INDEX_BACK(IW) = M%WALL_INDEX(IC,IOR)
      ENDIF
   ENDIF
ENDDO

! Set clocks and counters related to frequency of solid phase conduction updates

M%BC_CLOCK     = T_BEGIN
M%WALL_COUNTER = 0
 
! Allocate arrays for storing velocity boundary condition info
 
N_EDGES_DIM = 4*(IBP1*JBP1+IBP1*KBP1+JBP1*KBP1)
DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
   IPTS = OB%I2-OB%I1
   JPTS = OB%J2-OB%J1
   KPTS = OB%K2-OB%K1
   N_EDGES_DIM = N_EDGES_DIM + 4*(IPTS*JPTS+IPTS*KPTS+JPTS*KPTS)
ENDDO

ALLOCATE(M%IJKE(14,N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','IJKE',IZERO)   
M%IJKE  = 0
ALLOCATE(M%OME_E(N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','OME_E',IZERO)  
M%OME_E = 0._EB
ALLOCATE(M%TAU_E(N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','TAU_E',IZERO)  
M%TAU_E = 0._EB
ALLOCATE(M%ACTIVE_EDGE(N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','ACTIVE_EDGE',IZERO)  
M%ACTIVE_EDGE = .TRUE.
 
! Initialize and allocate lagrangian particle/droplet arrays

M%NLP = 0
M%NLPDIM = 1000
IF (DROPLET_FILE) THEN
   ALLOCATE(M%DROPLET(M%NLPDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','DROPLET',IZERO)
ENDIF
 
! Allocate array to hold character strings for Smokeview file
 
M%N_STRINGS     =   0
M%N_STRINGS_MAX = 100
ALLOCATE(M%STRING(M%N_STRINGS_MAX),STAT=IZERO)
CALL ChkMemErr('INIT','STRING',IZERO)
 
! Set up arrays to hold velocity boundary condition info
 
CALL INITIALIZE_EDGES
 
! Initialize Pressure solver
 
CALL INITIALIZE_POISSON_SOLVER
 
! Determine which wall cells to assign for solid phase thermocouples and profiles
 
CALL INITIALIZE_DEVC
CALL INITIALIZE_PROF
 
! Allocate and Initialize Mesh-Dependent Radiation Arrays
 
M%ANGLE_INC_COUNTER = 0
M%RAD_CALL_COUNTER  = 0
IF (RADIATION) THEN
   ALLOCATE(M%UIID(0:M%IBP1,0:M%JBP1,0:M%KBP1,1:UIIDIM),STAT=IZERO)
   M%UIID = 0.
ENDIF
CALL ChkMemErr('INIT','UIID',IZERO)
DO IW=1,M%NDWC
   IF (M%BOUNDARY_TYPE(IW)/=OPEN_BOUNDARY) THEN
      IF (RADIATION) THEN
         ALLOCATE(M%WALL(IW)%ILW(NRA,NSB),STAT=IZERO)
         CALL ChkMemErr('INIT','ILW',IZERO)
         M%WALL(IW)%ILW  = SIGMA*TMPA4*RPI
      ENDIF
      M%QRADIN(IW)  = M%E_WALL(IW)*SIGMA*TMPA4
      M%QRADOUT(IW) = M%E_WALL(IW)*SIGMA*TMPA4
   ENDIF
ENDDO

! Initialize Mesh Exchange
 
CALL INITIALIZE_INTERPOLATION
 
CONTAINS
 
 
SUBROUTINE INITIALIZE_EDGES
 
! Set up edge arrays for velocity boundary conditions
 
INTEGER I,J,K,N
 
CALL POINT_TO_MESH(NM)