init.f90

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

MODULE INIT      
 
! Allocate numerous arrays and perform miscellaneous initializations

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE MESH_POINTERS
USE GLOBAL_CONSTANTS
USE TRAN
USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR 
USE COMP_FUNCTIONS, ONLY : SHUTDOWN
USE DEVICE_VARIABLES
IMPLICIT NONE

PRIVATE
INTEGER IZERO
CHARACTER(255), PARAMETER :: initid='$Id: init.f90 708 2007-09-28 17:17:17Z drjfloyd $'
CHARACTER(255), PARAMETER :: initrev='$Revision: 708 $'
CHARACTER(255), PARAMETER :: initdate='$Date: 2007-09-28 13:17:17 -0400 (Fri, 28 Sep 2007) $'

CHARACTER(100) MESSAGE
PUBLIC INITIALIZE_MESH_VARIABLES,INITIALIZE_GLOBAL_VARIABLES, OPEN_AND_CLOSE, GET_REV_init

TYPE (MESH_TYPE), POINTER :: M
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (MATERIAL_TYPE), POINTER :: ML
 
 
CONTAINS
 
 
SUBROUTINE INITIALIZE_MESH_VARIABLES(NM)
USE RADCONS, ONLY: UIIDIM,NSB,NRA
USE CONTROL_VARIABLES
INTEGER :: N,I,J,K,II,JJ,KK,IPTS,JPTS,KPTS,N_EDGES_DIM,IW,IC,IBC,IOR,IOPZ
INTEGER, INTENT(IN) :: NM
REAL(EB) :: MU_N
INTEGER, POINTER :: IBP1, JBP1, KBP1,IBAR, JBAR, KBAR, NDWC, N_EDGES, NWC
REAL(EB),POINTER :: XS,XF,YS,YF,ZS,ZF
TYPE (INITIALIZATION_TYPE), POINTER :: IN
TYPE (P_ZONE_TYPE), POINTER :: PZ
 
M => MESHES(NM)
IBP1 =>M%IBP1
JBP1 =>M%JBP1
KBP1 =>M%KBP1
IBAR =>M%IBAR
JBAR =>M%JBAR
KBAR =>M%KBAR
NDWC =>M%NDWC
NWC =>M%NWC
N_EDGES=>M%N_EDGES
XS=>M%XS 
YS=>M%YS 
ZS=>M%ZS
XF=>M%XF 
YF=>M%YF 
ZF=>M%ZF
ALLOCATE(  M%RHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO',IZERO)
ALLOCATE( M%RHOS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHOS',IZERO)
ALLOCATE(  M%TMP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP',IZERO)
ALLOCATE( M%FRHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FRHO',IZERO)
ALLOCATE( M%U(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','U',IZERO)
ALLOCATE( M%V(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','V',IZERO)
ALLOCATE( M%W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','W',IZERO)
ALLOCATE(M%US(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','US',IZERO)
ALLOCATE(M%VS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','VS',IZERO)
ALLOCATE(M%WS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WS',IZERO)
ALLOCATE(M%FVX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVX',IZERO)
ALLOCATE(M%FVY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVY',IZERO)
ALLOCATE(M%FVZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVZ',IZERO)
ALLOCATE(   M%H(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','H',IZERO)
ALLOCATE(   M%HP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','HP',IZERO)
ALLOCATE(M%DDDT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DDDT',IZERO)
ALLOCATE(   M%D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','D',IZERO)
ALLOCATE(  M%DS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DS',IZERO)
ALLOCATE( M%MU(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MU',IZERO)
ALLOCATE(   M%Q(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','Q',IZERO)

! Background pressure, temperature, density as a function of height (Z coordinate)

ALLOCATE(  M%PBAR(0:KBP1,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','PBAR',IZERO)
ALLOCATE(  M%PBAR_S(0:KBP1,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','PBAR_S',IZERO)
ALLOCATE(  M%R_PBAR(0:KBP1,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','R_PBAR',IZERO)
ALLOCATE(  M%D_PBAR_DT(N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','D_PBAR_DT',IZERO)
ALLOCATE(  M%D_PBAR_S_DT(N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','D_PBAR_S_DT',IZERO)
ALLOCATE(M%P_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','P_0',IZERO)
ALLOCATE(M%TMP_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_0',IZERO)
ALLOCATE(M%RHO_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO_0',IZERO)

! Leaks

ALLOCATE(  M%LEAK_AREA(0:N_ZONE,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','LEAK_AREA',IZERO)
M%LEAK_AREA = 0._EB
ALLOCATE(  M%U_LEAK(0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','U_LEAK',IZERO)
M%U_LEAK = 0._EB
ALLOCATE(  M%FDS_LEAK_AREA(0:N_ZONE,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','FDS_LEAK_AREA',IZERO)
M%FDS_LEAK_AREA = 0._EB
! Allocate species arrays
 
IF (N_SPECIES>0) THEN
   ALLOCATE( M%YY(0:IBP1,0:JBP1,0:KBP1,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','YY',IZERO)
   M%YY = 0._EB
   ALLOCATE(M%YYS(0:IBP1,0:JBP1,0:KBP1,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','YYS',IZERO)
   M%YYS = 0._EB   
   ALLOCATE(M%DEL_RHO_D_DEL_Y(0:IBP1,0:JBP1,0:KBP1,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','DEL_RHO_D_DEL_Y',IZERO)
   M%DEL_RHO_D_DEL_Y = 0._EB   
   ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','RSUM',IZERO)
   M%RSUM = RSUM0
   ALLOCATE(M%Z_SUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','Z_SUM',IZERO)
   M%Z_SUM = 0._EB
   ALLOCATE(M%Y_SUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','Y_SUM',IZERO)
   M%Y_SUM = 0._EB   
ENDIF

! Allocate water mass arrays if sprinklers are present
 
IF (DROPLET_FILE) PARTICLE_TAG = NM

IF (N_EVAP_INDICIES>0) THEN
   ALLOCATE(M%AVG_DROP_DEN(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICIES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_DEN',IZERO) 
   M%AVG_DROP_DEN=0._EB
   ALLOCATE(M%AVG_DROP_TMP(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICIES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_TMP',IZERO) 
   M%AVG_DROP_TMP=TMPM
   ALLOCATE(M%AVG_DROP_RAD(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICIES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_RAD',IZERO) 
   M%AVG_DROP_RAD=0._EB
   ALLOCATE(M%QR_W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR_W',IZERO) 
   M%QR_W = 0._EB
   ALLOCATE(M%D_VAP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','D_VAP',IZERO) 
   M%D_VAP = 0._EB
ENDIF
 
! If radiation absorption desired allocate arrays
 
ALLOCATE(M%QR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','QR',IZERO)
M%QR = 0._EB
ALLOCATE(M%KAPPA(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','KAPPA',IZERO) 
M%KAPPA = 0._EB
ALLOCATE(M%UII(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','UII',IZERO)
M%UII = 0._EB

! Work arrays
ALLOCATE(M%WORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK1',IZERO)
ALLOCATE(M%WORK2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK2',IZERO)
ALLOCATE(M%WORK3(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK3',IZERO)
ALLOCATE(M%WORK4(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK4',IZERO)
ALLOCATE(M%WORK5(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK5',IZERO)
ALLOCATE(M%WORK6(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK6',IZERO)
ALLOCATE(M%WORK7(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK7',IZERO)

ALLOCATE(M%LOGICAL_WORK(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','LOGICAL_WORK',IZERO)
 
! Boundary file patch counter
 
ALLOCATE(M%INC(-3:3,0:M%N_OBST),STAT=IZERO)
CALL ChkMemErr('INIT','INC',IZERO)

! Initialize background pressure, temperature and density

M%D_PBAR_DT   = 0._EB
M%D_PBAR_S_DT = 0._EB

IF (STRATIFICATION) THEN
   DO K=0,M%KBP1
      M%TMP_0(K) = TMPA + LAPSE_RATE*M%ZC(K)
      IF (LAPSE_RATE /= 0._EB) THEN
         M%P_0(K) = P_INF*(M%TMP_0(K)/TMPA)**(GVEC(3)/RSUM0/LAPSE_RATE)
      ELSE
         M%P_0(K) = P_INF*EXP(GVEC(3)*M%ZC(K)/(RSUM0*TMPA))
      ENDIF
   ENDDO
ELSE
   M%TMP_0(:) = TMPA
   M%P_0(:)   = P_INF
ENDIF
DO K=0,M%KBP1
   M%PBAR(K,:)   = M%P_0(K)
   M%PBAR_S(K,:) = M%P_0(K)
   M%RHO_0(K)    = M%P_0(K)/(M%TMP_0(K)*RSUM0)
ENDDO

! Initialize various time step variables
 
M%DTOLD   = M%DT
M%DTNEXT  = M%DT
M%DTINT   = M%DT

! Initialize major arrays

DO K=0,M%KBP1
   M%RHO(:,:,K) = M%RHO_0(K)
   M%TMP(:,:,K) = M%TMP_0(K)
ENDDO
M%RHOS    = M%RHO
M%RHO_AVG = RHOA
M%FRHO    = 0._EB
M%U       = U0
M%V       = V0
M%W       = W0
M%US      = U0
M%VS      = V0
M%WS      = W0
IF (NOISE) CALL INITIAL_NOISE
M%FVX   = 0._EB
M%FVY   = 0._EB
M%FVZ   = 0._EB
M%H     = H0
M%DDDT  = 0._EB
M%D     = 0._EB
M%DS    = 0._EB
M%Q     = 0._EB
IF (N_SPECIES > 0) M%DEL_RHO_D_DEL_Y = 0._EB
 
! Upper bounds on local HRR per unit volume
 
IF (TWO_D) THEN
   M%Q_UPPER = HRRPUA_SHEET*1000._EB/(M%DXMIN*M%DZMIN)**0.5_EB
ELSE
   M%Q_UPPER = HRRPUA_SHEET*1000._EB/(M%DXMIN*M%DYMIN*M%DZMIN)**ONTH
ENDIF

! Viscosity
 
MU_N = SPECIES(0)%MU(NINT(0.1_EB*TMPA))
M%MU = MU_N
IF (DNS .AND. ISOTHERMAL) THEN
   ALLOCATE(M%RREDX(M%IBAR),STAT=IZERO)
   CALL ChkMemErr('READ','RREDX',IZERO)
   ALLOCATE(M%RREDY(M%JBAR),STAT=IZERO)
   CALL ChkMemErr('READ','RREDY',IZERO)
   ALLOCATE(M%RREDZ(M%KBAR),STAT=IZERO)
   CALL ChkMemErr('READ','RREDZ',IZERO)
   DO I=1,M%IBAR
      M%RREDX(I) = (MU_N/RHOA)*M%RDX(I)
   ENDDO
   DO J=1,M%JBAR
      M%RREDY(J) = (MU_N/RHOA)*M%RDY(J)
   ENDDO
   DO K=1,M%KBAR
      M%RREDZ(K) = (MU_N/RHOA)*M%RDZ(K)
   ENDDO
ENDIF
 
! Initialize mass fraction arrays

IF (N_SPECIES > 0) THEN
   M%RSUM = RSUM0
   DO N=1,N_SPECIES
      M%YY(:,:,:,N)  = SPECIES(N)%YY0
      M%YYS(:,:,:,N) = SPECIES(N)%YY0
      IF (SPECIES(N)%MODE==MIXTURE_FRACTION_SPECIES) M%Z_SUM = M%Z_SUM + SPECIES(N)%YY0   
      IF (SPECIES(N)%MODE==GAS_SPECIES) THEN 
         M%Y_SUM = M%Y_SUM + SPECIES(N)%YY0
         M%RSUM  = M%RSUM + (SPECIES(N)%RCON - SPECIES(0)%RCON)*SPECIES(N)%YY0
      ENDIF
   ENDDO
ENDIF

! Initialize pressure ZONEs
 
ALLOCATE(M%PRESSURE_ZONE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ZONE',IZERO)
M%PRESSURE_ZONE = 0
DO N=1,N_ZONE
   PZ => P_ZONE(N)
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            IF (M%XC(I) > PZ%X1 .AND. M%XC(I) < PZ%X2 .AND. &
                M%YC(J) > PZ%Y1 .AND. M%YC(J) < PZ%Y2 .AND. &
                M%ZC(K) > PZ%Z1 .AND. M%ZC(K) < PZ%Z2) THEN 
                M%PRESSURE_ZONE(I,J,K) = N
                DO IOPZ=0,N_ZONE
                   IF (PZ%LEAK_AREA(IOPZ) > 0._EB) M%LEAK_AREA(N,IOPZ) = PZ%LEAK_AREA(IOPZ)
                   IF (PZ%LEAK_AREA(IOPZ) > 0._EB) M%LEAK_AREA(IOPZ,N) = PZ%LEAK_AREA(IOPZ)
                ENDDO
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO

! Over-ride default ambient conditions with user-prescribed INITializations
 
DO N=1,N_INIT
   IN => INITIALIZATION(N)
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            IF (M%XC(I) > IN%X1 .AND. M%XC(I) < IN%X2 .AND. &
                M%YC(J) > IN%Y1 .AND. M%YC(J) < IN%Y2 .AND. &
                M%ZC(K) > IN%Z1 .AND. M%ZC(K) < IN%Z2) THEN
               M%TMP(I,J,K)            = IN%TEMPERATURE
               M%RHO(I,J,K)            = IN%DENSITY
               IF (N_SPECIES>0) M%YY(I,J,K,1:N_SPECIES) = IN%MASS_FRACTION(1:N_SPECIES)
               IF (IN%ADJUST_DENSITY)     M%RHO(I,J,K) = M%RHO(I,J,K)*M%P_0(K)/P_INF
               IF (IN%ADJUST_TEMPERATURE) M%TMP(I,J,K) = M%TMP(I,J,K)*M%P_0(K)/P_INF
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO

M%RHOS    = M%RHO

! Radiation

M%QR    = 0._EB
M%KAPPA = 0._EB
M%UII   = 4._EB*SIGMA*TMPA4
M%WORK1 = 0._EB
M%WORK2 = 0._EB
M%WORK3 = 0._EB
M%WORK4 = 0._EB
M%WORK5 = 0._EB
M%WORK6 = 0._EB
M%WORK7 = 0._EB

! Designate each boundary cell with a reference number for wall BC's
 
NWC  = 0
NDWC = 0
 
! Determine the number of wall cells to allocate
 
OBST_LOOP_1: DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
   IF (OB%CONSUMABLE) THEN
      NDWC = NDWC + 3*(OB%I2-OB%I1+1)*(OB%J2-OB%J1+1)*(OB%K2-OB%K1+1)
   ELSEIF (OB%DEVC_INDEX>0 .OR. OB%CTRL_INDEX>0 ) THEN
      NDWC = NDWC + 3*(OB%I2-OB%I1+1)*(OB%J2-OB%J1+1)*(OB%K2-OB%K1+1)   
!      IF (.NOT. DEVICE(OB%DEVC_INDEX)%INITIAL_STATE) NDWC = NDWC + 3*(OB%I2-OB%I1+1)*(OB%J2-OB%J1+1)*(OB%K2-OB%K1+1)
!   ELSEIF (OB%CTRL_INDEX/=0) THEN
!      IF (.NOT. CONTROL(OB%CTRL_INDEX)%INITIAL_STATE) NDWC = NDWC + 3*(OB%I2-OB%I1+1)*(OB%J2-OB%J1+1)*(OB%K2-OB%K1+1)
   ELSE
      DO K=OB%K1+1,OB%K2
         DO J=OB%J1+1,OB%J2
            IC = M%CELL_INDEX(OB%I1  ,J,K)
            IF (.NOT.M%SOLID(IC)) NDWC = NDWC + 1
            IC = M%CELL_INDEX(OB%I2+1,J,K)
            IF (.NOT.M%SOLID(IC)) NDWC = NDWC + 1
         ENDDO 
      ENDDO
      DO K=OB%K1+1,OB%K2
         DO I=OB%I1+1,OB%I2
            IC = M%CELL_INDEX(I,OB%J1  ,K)
            IF (.NOT.M%SOLID(IC)) NDWC = NDWC + 1
            IC = M%CELL_INDEX(I,OB%J2+1,K)
            IF (.NOT.M%SOLID(IC)) NDWC = NDWC + 1
         ENDDO 
      ENDDO
      DO J=OB%J1+1,OB%J2
         DO I=OB%I1+1,OB%I2
            IC = M%CELL_INDEX(I,J,OB%K1  )
            IF (.NOT.M%SOLID(IC)) NDWC = NDWC + 1
            IC = M%CELL_INDEX(I,J,OB%K2+1)
            IF (.NOT.M%SOLID(IC)) NDWC = NDWC + 1
         ENDDO 
      ENDDO
   ENDIF
ENDDO OBST_LOOP_1

NDWC = NDWC + M%NEWC
ALLOCATE(M%WALL(1:NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL',IZERO)
 
ALLOCATE(M%TMP_F(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_F',IZERO) 
M%TMP_F = TMPA
ALLOCATE(M%TMP_B(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_B',IZERO) 
M%TMP_B = TMPA
ALLOCATE(M%TMP_W(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_W',IZERO) 
M%TMP_W = TMPA
IF (N_EVAP_INDICIES > 0) THEN
   ALLOCATE(M%RCP_W(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','RCP_W',IZERO) 
   M%RCP_W = 0._EB
ENDIF