init.f90

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

      IF (.NOT.SF%THERMALLY_THICK) THEN
         WRITE(MESSAGE,'(A,I3,A)') 'PROFile ',N, ' must be associated with a heat-conducting surface'
         CALL SHUTDOWN(MESSAGE)
      ENDIF
      IF (PF%QUANTITY /= 'TEMPERATURE' .AND. PF%QUANTITY /= 'DENSITY') THEN
         SUCCESS = .FALSE.
         DO NN=1,SF%N_MATL
            IF (PF%QUANTITY==SF%MATL_NAME(NN)) SUCCESS = .TRUE.
         ENDDO
         IF (.NOT.SUCCESS) THEN
            WRITE(MESSAGE,'(A,A,A)') 'QUANTITY ',TRIM(PF%QUANTITY), ' is not appropriate for the designated location'
            CALL SHUTDOWN(MESSAGE)
         ENDIF
      ENDIF
   ELSE
      WRITE(MESSAGE,'(A,I4,A)') 'Reposition PROF No.',PF%ORDINAL, '. FDS cannot determine which boundary cell to assign'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
ENDDO PROF_LOOP
 
END SUBROUTINE INITIALIZE_PROF


SUBROUTINE GET_IW(II,JJ,KK,IOR,IW)

INTEGER :: II,JJ,KK,IOR,IW

IC  = M%CELL_INDEX(II,JJ,KK)
 
IF (M%SOLID(IC)) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      II = II-1
      CASE( 1)
         IF (II<M%IBP1) II = II+1
      CASE(-2)
         IF (JJ>0)      JJ = JJ-1
      CASE( 2)
         IF (JJ<M%JBP1) JJ = JJ+1
      CASE(-3)
         IF (KK>0)      KK = KK-1
      CASE( 3)
         IF (KK<M%KBP1) KK = KK+1
   END SELECT
ENDIF
 
IC  = M%CELL_INDEX(II,JJ,KK)
IW  = M%WALL_INDEX(IC,-IOR)
 
IF (IW<=0) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      IC = M%CELL_INDEX(II-1,JJ,KK)
      CASE( 1)
         IF (II<M%IBP1) IC = M%CELL_INDEX(II+1,JJ,KK)
      CASE(-2)
         IF (JJ>0)      IC = M%CELL_INDEX(II,JJ-1,KK)
      CASE( 2)
         IF (JJ<M%JBP1) IC = M%CELL_INDEX(II,JJ+1,KK)
      CASE(-3)
         IF (KK>0)      IC = M%CELL_INDEX(II,JJ,KK-1)
      CASE( 3)
         IF (KK<M%KBP1) IC = M%CELL_INDEX(II,JJ,KK+1)
   END SELECT
   IW = M%WALL_INDEX(IC,-IOR)
ENDIF

END SUBROUTINE GET_IW



SUBROUTINE INITIALIZE_INTERPOLATION
 
! Create arrays by which info is to exchanged across meshes
 
INTEGER :: NOM,I,J,K
TYPE (MESH_TYPE), POINTER :: M2
 
IF (NM==1) RETURN
 
ALLOCATE(M%INTERPOLATED_MESH(1:M%IBAR,1:M%JBAR,1:M%KBAR),  STAT=IZERO)
CALL ChkMemErr('INIT','INTERPOLATED_MESH',IZERO)  
M%INTERPOLATED_MESH = 0
 
DO K=1,M%KBAR
   DO J=1,M%JBAR
      DO I=1,M%IBAR
         OTHER_MESH_LOOP: DO NOM=1,NM-1
            M2=>MESHES(NOM)
            IF (M%X(I-1)>=M2%XS .AND. M%X(I)<=M2%XF .AND.  M%Y(J-1)>=M2%YS .AND. M%Y(J)<=M2%YF .AND. &
                M%Z(K-1)>=M2%ZS .AND. M%Z(K)<=M2%ZF) THEN 
               M%INTERPOLATED_MESH(I,J,K) = NOM
               EXIT OTHER_MESH_LOOP
            ENDIF
         ENDDO OTHER_MESH_LOOP
      ENDDO
   ENDDO
ENDDO
 
END SUBROUTINE INITIALIZE_INTERPOLATION
 
END SUBROUTINE INITIALIZE_MESH_VARIABLES
 
 
 
SUBROUTINE INITIALIZE_GLOBAL_VARIABLES
USE CONTROL_VARIABLES, ONLY: N_CTRL
! Initialize time, printout and plot clocks
 
ALLOCATE(PART_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PART_CLOCK',IZERO) 
ALLOCATE(ISOF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','ISOF_CLOCK',IZERO) 
ALLOCATE(BNDF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','BNDF_CLOCK',IZERO) 
ALLOCATE(SLCF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','SLCF_CLOCK',IZERO) 
ALLOCATE(CORE_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','CORE_CLOCK',IZERO) 
ALLOCATE(PL3D_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PL3D_CLOCK',IZERO) 
ALLOCATE(PROF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PROF_CLOCK',IZERO) 
 
ICYC              = 0
PART_CLOCK        = T_BEGIN
DEVC_CLOCK        = T_BEGIN      
CTRL_CLOCK        = T_BEGIN      
PROF_CLOCK        = T_BEGIN      
PL3D_CLOCK        = T_BEGIN + DT_PL3D
ISOF_CLOCK        = T_BEGIN      
SLCF_CLOCK        = T_BEGIN   
BNDF_CLOCK        = T_BEGIN     
CORE_CLOCK        = T_BEGIN + DT_RESTART
HRR_CLOCK         = T_BEGIN
MINT_CLOCK        = T_BEGIN      

IF (N_DEVC==0) DEVC_CLOCK = 1.E10_EB
IF (N_CTRL==0) CTRL_CLOCK = 1.E10_EB
IF (N_PROF==0) PROF_CLOCK = 1.E10_EB
IF (N_ISOF==0) ISOF_CLOCK = 1.E10_EB
IF (N_BNDF==0) BNDF_CLOCK = 1.E10_EB

ALLOCATE(HRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','HRR',IZERO) 
HRR = 0._EB
ALLOCATE(RHRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','RHRR',IZERO) 
RHRR = 0._EB
ALLOCATE(CHRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','CHRR',IZERO) 
CHRR = 0._EB
ALLOCATE(FHRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FHRR',IZERO) 
FHRR = 0._EB
ALLOCATE(MLR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MLR',IZERO) 
MLR = 0._EB
ALLOCATE(HRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','HRR_SUM',IZERO)  
HRR_SUM=0._EB
ALLOCATE(RHRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','RHRR_SUM',IZERO) 
RHRR_SUM=0._EB
ALLOCATE(CHRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','CHRR_SUM',IZERO) 
CHRR_SUM=0._EB
ALLOCATE(FHRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FHRR_SUM',IZERO) 
FHRR_SUM=0._EB
ALLOCATE(MLR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MLR_SUM',IZERO)  
MLR_SUM=0._EB
ALLOCATE(HRR_COUNT(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MLR_COUNT',IZERO)  
HRR_COUNT=0._EB
ALLOCATE(MINT(0:MAX_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT',IZERO) 
ALLOCATE(MINT_SUM(0:MAX_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT_SUM',IZERO) 
MINT_SUM=0._EB
ALLOCATE(MINT_COUNT(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT_COUNT',IZERO) 
MINT_COUNT=0._EB
ALLOCATE(I_MIN(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','I_MIN',IZERO) 
I_MIN = -10
ALLOCATE(I_MAX(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','I_MAX',IZERO) 
I_MAX = -10
ALLOCATE(J_MIN(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','J_MIN',IZERO) 
J_MIN = -10
ALLOCATE(J_MAX(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','J_MAX',IZERO) 
J_MAX = -10
ALLOCATE(K_MIN(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','K_MIN',IZERO) 
K_MIN = -10
ALLOCATE(K_MAX(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','K_MAX',IZERO) 
K_MAX = -10
ALLOCATE(NIC(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','NIC',IZERO) 
NIC = 0
ALLOCATE(T_PER_STEP(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','T_PER_STEP',IZERO) 
T_PER_STEP = 0._EB
ALLOCATE(T_ACCUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','T_ACCUM',IZERO) 
T_ACCUM = 0._EB
ALLOCATE(NTCYC(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','NTCYC',IZERO) 
NTCYC = 0
ALLOCATE(NCYC(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','NCYC',IZERO) 
NCYC = 0

ALLOCATE(DSUM(N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','DSUM',IZERO) 
DSUM = 0._EB
ALLOCATE(PSUM(N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PSUM',IZERO) 
PSUM = 0._EB
ALLOCATE(USUM(N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','USUM',IZERO) 
USUM = 0._EB
 
END SUBROUTINE INITIALIZE_GLOBAL_VARIABLES
 
 
SUBROUTINE INIT_WALL_CELL(NM,I,J,K,I_OBST,IW,IOR,IBC)
 
INTEGER  :: NM,NOM,ICO,IBC,IOR
REAL(EB) :: PX,PY,PZ,X1,X2,Y1,Y2,Z1,Z2,T_ACTIVATE,XI,YJ,ZK, XIN,YIN,ZIN,DIST
INTEGER N,NN,I_OBST,I,J,K,IBCX,IIG,JJG,KKG,IW,II,JJ,KK,ICG,IL
LOGICAL :: VENT_FOUND
TYPE (MESH_TYPE), POINTER :: MM
TYPE (WALL_TYPE), POINTER :: WC
TYPE (OBSTRUCTION_TYPE), POINTER :: OBX
TYPE (VENTS_TYPE), POINTER :: VT
 
M=>MESHES(NM)
 
! Compute boundary cell physical coords (XW,YW,ZW) and area (AW)
 
IBCX         = IBC
M%IJKW(1,IW) = I
M%IJKW(2,IW) = J
M%IJKW(3,IW) = K
M%IJKW(4,IW) = IOR
M%IJKW(5,IW) = IBCX
M%OBST_INDEX_W(IW) = I_OBST
 
IF (ABS(IOR)==1) THEN
   IF (IOR== 1) THEN
      M%XW(IW)     = M%X(I)
      M%IJKW(6,IW) = I+1
      M%RDN(IW)    = M%RDXN(I) 
      M%AW(IW)     = M%R(I)*M%DY(J)*M%DZ(K)
      M%UW(IW)     = -U0
   ENDIF
   IF (IOR==-1) THEN
      M%XW(IW)     = M%X(I-1)
      M%IJKW(6,IW) = I-1
      M%RDN(IW)    = M%RDXN(I-1) 
      M%AW(IW)     = M%R(I-1)*M%DY(J)*M%DZ(K)
      M%UW(IW)     = U0
   ENDIF
   M%IJKW(7,IW) = J
   M%IJKW(8,IW) = K
   M%YW(IW) = 0.5_EB*(M%Y(J)+M%Y(J-1))  
   M%ZW(IW) = 0.5_EB*(M%Z(K)+M%Z(K-1))
ENDIF
IF (ABS(IOR)==2) THEN
   IF (IOR== 2) THEN
      M%YW(IW)     = M%Y(J)
      M%IJKW(7,IW) = J+1
      M%RDN(IW)    = M%RDYN(J) 
      M%UW(IW)     = -V0
   ENDIF
   IF (IOR==-2) THEN
      M%YW(IW)     = M%Y(J-1)
      M%IJKW(7,IW) = J-1
      M%RDN(IW)    = M%RDYN(J-1) 
      M%UW(IW)     = V0
   ENDIF
   M%IJKW(6,IW) = I
   M%IJKW(8,IW) = K
   M%XW(IW) = 0.5_EB*(M%X(I)+M%X(I-1))
   M%ZW(IW) = 0.5_EB*(M%Z(K)+M%Z(K-1))
   M%AW(IW) = M%DX(I)*M%DZ(K)
ENDIF
IF (ABS(IOR)==3) THEN
   IF (IOR== 3) THEN
      M%ZW(IW)     = M%Z(K)
      M%IJKW(8,IW) = K+1
      M%RDN(IW)    = M%RDZN(K) 
      M%UW(IW)     = -W0
   ENDIF
   IF (IOR==-3) THEN
      M%ZW(IW)     = M%Z(K-1)
      M%IJKW(8,IW) = K-1
      M%RDN(IW)    = M%RDZN(K-1) 
      M%UW(IW)     = W0
   ENDIF
   M%IJKW(6,IW) = I
   M%IJKW(7,IW) = J
   M%XW(IW) = 0.5_EB*(M%X(I)+M%X(I-1))
   M%YW(IW) = 0.5_EB*(M%Y(J)+M%Y(J-1))
   M%AW(IW) = M%DX(I)*M%RC(I)*M%DY(J)
ENDIF
 
IF (M%AW(IW)>0._EB) M%RAW(IW) = 1._EB/M%AW(IW)

! Gas phase cell abutting boundary cell

IIG = M%IJKW(6,IW)
JJG = M%IJKW(7,IW)
KKG = M%IJKW(8,IW)

! Use IWA to indicate the boundary cell number, IW, at the various faces of a grid cell
 
ICG = M%CELL_INDEX(IIG,JJG,KKG)
M%WALL_INDEX(ICG,-IOR) = IW

! Use BOUNDARY_TYPE to indicate whether the boundary cell is blocked or on an obstruction that is HIDDEN

M%BOUNDARY_TYPE(IW) = SOLID_BOUNDARY
IF (M%SOLID(ICG)) M%BOUNDARY_TYPE(IW) = NULL_BOUNDARY
IF (I_OBST > 0 .AND. M%OBSTRUCTION(I_OBST)%HIDDEN) M%BOUNDARY_TYPE(IW) = NULL_BOUNDARY

IF (SURFACE(IBC)%POROUS .AND. M%BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) M%BOUNDARY_TYPE(IW) = POROUS_BOUNDARY

! Assign the ZONE number to all boundary cells

M%PRESSURE_ZONE_WALL(IW) = M%PRESSURE_ZONE(IIG,JJG,KKG)
 
! Check for neighboring meshes in a multiple mesh calculation
 
CHECK_MESHES: IF (IW<=M%NEWC .AND.  .NOT.EVACUATION_ONLY(NM)) THEN
 
   XIN = M%XW(IW) 
   YIN = M%YW(IW) 
   ZIN = M%ZW(IW) 
 
   FIND_MESH: DO NOM=1,NMESHES
      IF (NOM==NM) CYCLE FIND_MESH
      IF (EVACUATION_ONLY(NOM)) CYCLE FIND_MESH
      MM=>MESHES(NOM)
      IF (XIN<MM%XS) CYCLE FIND_MESH
      IF (XIN>MM%XF) CYCLE FIND_MESH
      IF (YIN<MM%YS) CYCLE FIND_MESH
      IF (YIN>MM%YF) CYCLE FIND_MESH
      IF (ZIN<MM%ZS) CYCLE FIND_MESH
      IF (ZIN>MM%ZF) CYCLE FIND_MESH
      SELECT CASE(IOR)
         CASE( 1)  
            XIN = XIN - 0.05_EB*M%DX(0)
         CASE(-1) 
            XIN = XIN + 0.05_EB*M%DX(M%IBP1)
         CASE( 2) 
            YIN = YIN - 0.05_EB*M%DY(0)
         CASE(-2)
            YIN = YIN + 0.05_EB*M%DY(M%JBP1)
         CASE( 3)
            ZIN = ZIN - 0.05_EB*M%DZ(0)
         CASE(-3) 
            ZIN = ZIN + 0.05_EB*M%DZ(M%KBP1)
      END SELECT
      IF (XIN<MM%XS) CYCLE FIND_MESH
      IF (XIN>MM%XF) CYCLE FIND_MESH
      IF (YIN<MM%YS) CYCLE FIND_MESH
      IF (YIN>MM%YF) CYCLE FIND_MESH
      IF (ZIN<MM%ZS) CYCLE FIND_MESH
      IF (ZIN>MM%ZF) CYCLE FIND_MESH
      XI = MM%CELLSI(NINT((XIN-MM%XS)*MM%RDXINT)) + 1._EB
      YJ = MM%CELLSJ(NINT((YIN-MM%YS)*MM%RDYINT)) + 1._EB
      ZK = MM%CELLSK(NINT((ZIN-MM%ZS)*MM%RDZINT)) + 1._EB
      II = FLOOR(XI)
      JJ = FLOOR(YJ)
      KK = FLOOR(ZK)
      M%INTERPOLATION_FACTOR(IW,1) = XI-II
      M%INTERPOLATION_FACTOR(IW,2) = YJ-JJ
      M%INTERPOLATION_FACTOR(IW,3) = ZK-KK
      M%CELL_VOLUME_RATIO(IW) = M%DX(IIG)*M%RC(IIG)*M%DY(JJG)*M%DZ(KKG)/ &
                                MM%DX(II)*MM%RC(II)*MM%DY(JJ)*MM%DZ(KK)
      M%IJKW(9,IW)  = NOM
      M%IJKW(10,IW) = II 
      M%IJKW(11,IW) = JJ 
      M%IJKW(12,IW) = KK 
      IF (I_OBST==0 .AND. .NOT.M%SOLID(ICG)) THEN
         M%BOUNDARY_TYPE(IW) = INTERPOLATED_BOUNDARY
         IBCX     = INTERPOLATED_SURF_INDEX
         M%IJKW(5,IW) = IBCX
      ENDIF
      ICO = MM%CELL_INDEX(II,JJ,KK)
      IF (.NOT.MM%SOLID(ICO) .AND. .NOT.M%SOLID(ICG))   M%SOLID(M%CELL_INDEX(I,J,K)) = .FALSE.
      EXIT FIND_MESH
   ENDDO FIND_MESH
 
ENDIF CHECK_MESHES
 
! Assign internal values of temp, density, and mass fraction
 
IF (N_SPECIES>0) THEN
   M%RSUM_W(IW)= M%RSUM(IIG,JJG,KKG)
   M%RSUM(I,J,K) = M%RSUM(IIG,JJG,KKG)
   M%YY_W(IW,1:N_SPECIES)  = M%YY(IIG,JJG,KKG,1:N_SPECIES)
   M%YY(I,J,K,1:N_SPECIES) = M%YY(IIG,JJG,KKG,1:N_SPECIES)
ENDIF
 
M%RHO_W(IW) = M%RHO(IIG,JJG,KKG)
M%TMP_W(IW) = M%TMP(IIG,JJG,KKG)
 
! Assign various other quantities to the cell
 
IF (I_OBST>0) THEN
   OBX=>M%OBSTRUCTION(I_OBST)
   M%AREA_ADJUST(IW) = OBX%INPUT_AREA(ABS(IOR))/OBX%FDS_AREA(ABS(IOR))
   IF (M%AREA_ADJUST(IW)==0._EB) M%AREA_ADJUST(IW) = 1._EB
   OBX%MASS = SURFACE(IBC)%SURFACE_DENSITY*M%AW(IW)*M%AREA_ADJUST(IW)
ENDIF
 
! Prescribe exit velocity for surface cell
 
IF (SURFACE(IBCX)%VEL/=-999._EB) THEN
   M%UW0(IW) = SURFACE(IBCX)%VEL
ELSE
   M%UW0(IW) = 0._EB
ENDIF
IF (I_OBST>0 .AND. SURFACE(IBCX)%VOLUME_FLUX/=-999._EB) THEN
   OBX=>M%OBSTRUCTION(I_OBST)
   M%UW0(IW) = SURFACE(IBCX)%VOLUME_FLUX/OBX%FDS_AREA(ABS(IOR))