func.f90

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

MODULE COMP_FUNCTIONS

! I/O + OS functions

USE PRECISION_PARAMETERS 
IMPLICIT NONE 
CHARACTER(255), PARAMETER :: funcid='$Id: func.f90 713 2007-09-28 21:44:10Z mcgratta $'
CHARACTER(255), PARAMETER :: funcrev='$Revision: 713 $'
CHARACTER(255), PARAMETER :: funcdate='$Date: 2007-09-28 17:44:10 -0400 (Fri, 28 Sep 2007) $'
 
CONTAINS

SUBROUTINE GET_REV_func(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') funcrev(INDEX(funcrev,':')+1:LEN_TRIM(funcrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') funcdate

END SUBROUTINE GET_REV_func

REAL(EB) FUNCTION SECOND()  ! Returns the CPU time in seconds.
REAL(FB) CPUTIME
CALL CPU_TIME(CPUTIME)
SECOND = CPUTIME
END FUNCTION SECOND


REAL(EB) FUNCTION WALL_CLOCK_TIME()  ! Returns the number of seconds since January 1, 2000, including leap years

! Thirty days hath September,
! April, June, and November.
! February has twenty-eight alone;
! All the rest have thirty-one,
! Excepting Leap-Year, that's the time
! When February's days are twenty-nine.

INTEGER :: DATE_TIME(8),WALL_CLOCK_SECONDS
CHARACTER(10) :: BIG_BEN(3)
! X_1 = common year, X_2 = leap year
INTEGER, PARAMETER :: S_PER_YEAR_1=31536000,S_PER_YEAR_2=31622400,S_PER_DAY=86400,S_PER_HOUR=3600,S_PER_MIN=60
INTEGER, PARAMETER, DIMENSION(12) :: ACCUMULATED_DAYS_1=(/0,31,59,90,120,151,181,212,243,273,304,334/), & 
                                     ACCUMULATED_DAYS_2=(/0,31,60,91,121,152,182,213,244,274,305,335/)
INTEGER :: YEAR_COUNT

CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
WALL_CLOCK_SECONDS = 0._EB
DO YEAR_COUNT=2001,DATE_TIME(1)
   !Leap year if divisible by 4 but not 100 unless by 400 (1900 no, 1904  yes, 2000 yes)
   IF (MOD(YEAR_COUNT,4)==0 .AND. (MOD(YEAR_COUNT,100)/=0 .OR. MOD(YEAR_COUNT,400)==0)) THEN
      WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_YEAR_2
   ELSE
      WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_YEAR_1
   ENDIF
ENDDO 
IF (MOD(DATE_TIME(1),4)==0 .AND. (MOD(DATE_TIME(1),100)/=0 .OR. MOD(DATE_TIME(1),400)==0 )) THEN
   WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_DAY*(ACCUMULATED_DAYS_2(DATE_TIME(2))+DATE_TIME(3))
ELSE
   WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_DAY*(ACCUMULATED_DAYS_1(DATE_TIME(2))+DATE_TIME(3))
ENDIF
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS +  S_PER_HOUR*DATE_TIME(5) + S_PER_MIN*DATE_TIME(6) + DATE_TIME(7)
WALL_CLOCK_TIME    = WALL_CLOCK_SECONDS + DATE_TIME(8)*0.001_EB

END FUNCTION WALL_CLOCK_TIME


SUBROUTINE SHUTDOWN(MESSAGE)  ! Stops the code gracefully after writing a message
USE GLOBAL_CONSTANTS, ONLY: LU_ERR
CHARACTER(*) MESSAGE
WRITE(LU_ERR,'(/A)') TRIM(MESSAGE)
STOP
END SUBROUTINE SHUTDOWN
 

!!SUBROUTINE FLUSH_BUFFER(UNIT) ! FLUSH_BUFFER flushes the buffer for the named logical unit.
!!USE IFPORT  ! For Intel compilers
!!USE GLOBAL_CONSTANTS, ONLY: FLUSH_FILE_BUFFERS
!!INTEGER, INTENT(IN) :: UNIT
!!IF (FLUSH_FILE_BUFFERS) CALL FLUSH(UNIT)
!!END SUBROUTINE FLUSH_BUFFER
 

SUBROUTINE GET_INPUT_FILE ! Read the argument after the command
USE GLOBAL_CONSTANTS, ONLY: FN_INPUT
IF (FN_INPUT=='null') CALL GETARG(1,FN_INPUT)
END SUBROUTINE GET_INPUT_FILE


SUBROUTINE CHECKREAD(NAME,LU,IOS)
 
INTEGER :: II,IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(4), INTENT(IN) :: NAME
CHARACTER(80) TEXT
 
IOS = 1
READLOOP: DO
   READ(LU,'(A)',END=10) TEXT
   TLOOP: DO II=1,72
      IF (TEXT(II:II)/='&' .AND. TEXT(II:II)/=' ') EXIT TLOOP
      IF (TEXT(II:II)=='&') THEN
         IF (TEXT(II+1:II+4)==NAME) THEN
            BACKSPACE(LU)
            IOS = 0
            EXIT READLOOP
         ELSE
            CYCLE READLOOP
         ENDIF
      ENDIF
   ENDDO TLOOP
ENDDO READLOOP
 
10 RETURN
END SUBROUTINE CHECKREAD

END MODULE COMP_FUNCTIONS



MODULE MEMORY_FUNCTIONS
USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE COMP_FUNCTIONS, ONLY : SHUTDOWN
IMPLICIT NONE

CONTAINS


SUBROUTINE ChkMemErr(CodeSect,VarName,IZERO)
 
! Memory checking routine
 
CHARACTER(*), INTENT(IN) :: CodeSect, VarName
INTEGER IZERO
CHARACTER(100) MESSAGE
 
IF (IZERO==0) RETURN
 
WRITE(MESSAGE,'(4A)') 'ERROR: Memory allocation failed for ', TRIM(VarName),' in the routine ',TRIM(CodeSect)
CALL SHUTDOWN(MESSAGE)

END SUBROUTINE ChkMemErr


SUBROUTINE RE_ALLOCATE_DROPLETS(CODE,NM,NOM,NEW_DROPS)
 
TYPE (DROPLET_TYPE), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: CODE,NM,NOM,NEW_DROPS
TYPE (MESH_TYPE), POINTER :: M
TYPE(OMESH_TYPE), POINTER :: M2
 
SELECT CASE(CODE)
   CASE(1)
      M=>MESHES(NM)
      ALLOCATE(DUMMY(1:M%NLPDIM))
      DUMMY = M%DROPLET
      DEALLOCATE(M%DROPLET)
      ALLOCATE(M%DROPLET(M%NLPDIM+NEW_DROPS))
      M%DROPLET(1:M%NLPDIM) = DUMMY(1:M%NLPDIM)
      M%NLPDIM = M%NLPDIM+NEW_DROPS
   CASE(2)
      M2=>MESHES(NM)%OMESH(NOM)
      ALLOCATE(DUMMY(1:M2%N_DROP_ORPHANS_DIM))
      DUMMY = M2%DROPLET
      DEALLOCATE(M2%DROPLET)
      ALLOCATE(M2%DROPLET(M2%N_DROP_ORPHANS_DIM+NEW_DROPS))
      M2%DROPLET(1:M2%N_DROP_ORPHANS_DIM) = DUMMY(1:M2%N_DROP_ORPHANS_DIM)
      M2%N_DROP_ORPHANS_DIM = M2%N_DROP_ORPHANS_DIM + NEW_DROPS
END SELECT
DEALLOCATE(DUMMY)

END SUBROUTINE RE_ALLOCATE_DROPLETS

 
FUNCTION REALLOCATE(P,N1,N2)          
REAL(EB), POINTER, DIMENSION(:) :: P, REALLOCATE
INTEGER, INTENT(IN) :: N1,N2
INTEGER :: NOLD, IERR
CHARACTER(100) :: MESSAGE
ALLOCATE(REALLOCATE(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE'
   CALL SHUTDOWN(MESSAGE)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE(N1:NOLD+N1-1) = P(N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P) 
END FUNCTION REALLOCATE


SUBROUTINE RE_ALLOCATE_STRINGS(NM)
 
CHARACTER(50), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M
 
M=>MESHES(NM)
ALLOCATE(DUMMY(1:M%N_STRINGS))
DUMMY = M%STRING
DEALLOCATE(M%STRING)
ALLOCATE(M%STRING(M%N_STRINGS_MAX+100))
M%STRING(1:M%N_STRINGS) = DUMMY(1:M%N_STRINGS)
M%N_STRINGS_MAX = M%N_STRINGS_MAX+100
DEALLOCATE(DUMMY)
 
END SUBROUTINE RE_ALLOCATE_STRINGS

END MODULE MEMORY_FUNCTIONS 

MODULE GEOMETRY_FUNCTIONS

! Functions for manipulating geometry

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
IMPLICIT NONE
 

CONTAINS
 

SUBROUTINE BLOCK_CELL(NM,I1,I2,J1,J2,K1,K2,IVAL,OBST_INDEX)

! Indicate which cells are blocked off
 
INTEGER :: NM,I1,I2,J1,J2,K1,K2,IVAL,I,J,K,OBST_INDEX,IC
TYPE (MESH_TYPE), POINTER :: M
 
M => MESHES(NM)
DO K=K1,K2
   DO J=J1,J2
      DO I=I1,I2
         IC = M%CELL_INDEX(I,J,K)
         SELECT CASE(IVAL)
            CASE(0) 
               M%SOLID(IC)        = .FALSE.
               M%OBST_INDEX_C(IC) = 0
            CASE(1)
               M%SOLID(IC)        = .TRUE.
               M%OBST_INDEX_C(IC) = OBST_INDEX
         END SELECT
      ENDDO
   ENDDO
ENDDO
 
END SUBROUTINE BLOCK_CELL
 
 

SUBROUTINE GET_N_LAYER_CELLS(DIFFUSIVITY,THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR,N_CELLS,DXMIN)

! Get number of wall cells in the layer

INTEGER, INTENT(OUT)  :: N_CELLS
REAL(EB), INTENT(OUT) :: DXMIN
REAL(EB), INTENT(IN)  :: DIFFUSIVITY,THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR
REAL(EB) :: DSUM
INTEGER  :: N, I

IF (THICKNESS.EQ.0._EB) THEN
   N_CELLS = 0
   DXMIN = 0._EB
   RETURN
ENDIF
SHRINK_LOOP: DO N=1,999
   DSUM = 0._EB
   SUM_LOOP: DO I=1,N
      DSUM = DSUM + STRETCH_FACTOR**(MIN(I-1,N-I))
   ENDDO SUM_LOOP
   IF ((THICKNESS/DSUM < CELL_SIZE_FACTOR*SQRT(DIFFUSIVITY)) .OR. (N==999)) THEN
      N_CELLS = N
      DXMIN = THICKNESS/DSUM
      EXIT SHRINK_LOOP
   ENDIF
ENDDO SHRINK_LOOP

END SUBROUTINE GET_N_LAYER_CELLS


SUBROUTINE GET_WALL_NODE_COORDINATES(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
         SMALLEST_CELL_SIZE,STRETCH_FACTOR,X_S)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS,N_LAYERS, N_LAYER_CELLS(N_LAYERS)
REAL(EB), INTENT(IN)    :: SMALLEST_CELL_SIZE(N_LAYERS),STRETCH_FACTOR
REAL(EB), INTENT(OUT)   :: X_S(0:N_CELLS)

INTEGER I, II, NL
REAL(EB) DX_S

   II = 0
   X_S(0) = 0._EB
   DO NL=1,N_LAYERS
      DO I=1,N_LAYER_CELLS(NL)
         II = II+1
         DX_S = SMALLEST_CELL_SIZE(NL)*STRETCH_FACTOR**(MIN(I-1,N_LAYER_CELLS(NL)-I))
         X_S(II) = X_S(II-1) + DX_S
      ENDDO
   ENDDO

END SUBROUTINE GET_WALL_NODE_COORDINATES


SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
         THICKNESS,GEOMETRY,X_S,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS, N_LAYERS, N_LAYER_CELLS(N_LAYERS), GEOMETRY
REAL(EB), INTENT(IN)    :: X_S(0:N_CELLS),THICKNESS
INTEGER, INTENT(OUT)    :: LAYER_INDEX(0:N_CELLS+1)
REAL(EB), INTENT(OUT)   :: DX(1:N_CELLS),RDX(0:N_CELLS+1),RDXN(0:N_CELLS),DX_WGT(0:N_CELLS),DXF,DXB

INTEGER I, II, NL, I_GRAD
REAL(EB) R

SELECT CASE(GEOMETRY)
CASE(SURF_CARTESIAN)
   I_GRAD = 0
CASE(SURF_CYLINDRICAL)
   I_GRAD = 1
CASE(SURF_SPHERICAL)
   I_GRAD = 2
END SELECT

   II = 0
   DO NL=1,N_LAYERS
      DO I=1,N_LAYER_CELLS(NL)
         II = II + 1
         LAYER_INDEX(II) = NL
      ENDDO
   ENDDO
   LAYER_INDEX(0) = 1
   LAYER_INDEX(N_CELLS+1) = N_LAYERS
   DXF = X_S(1)       - X_S(0)
   DXB = X_S(N_CELLS) - X_S(N_CELLS-1)

! Compute dx_weight for each node (dx_weight is the ratio of cell size to the 
! combined size of the current and next cell)

   DO I=1,N_CELLS-1
      DX_WGT(I) = (X_S(I)-X_S(I-1))/(X_S(I+1)-X_S(I-1))
   ENDDO
   DX_WGT(0)       = 0.5_EB
   DX_WGT(N_CELLS) = 0.5_EB

! Compute dx and 1/dx for each node (dx is the distance from cell boundary to cell boundary)

   DO I=1,N_CELLS
      DX(I)  = X_S(I)-X_S(I-1)
      RDX(I) = 1._EB/DX(I)
   ENDDO
  ! Adjust 1/dx_n to 1/rdr for cylindrical case and 1/r^2dr for spaherical
   IF (GEOMETRY /=SURF_CARTESIAN) THEN
      DO I=1,N_CELLS
         R = THICKNESS-0.5_EB*(X_S(I)+X_S(I-1))
         RDX(I) = RDX(I)/R**I_GRAD
      ENDDO
   ENDIF
   RDX(0)         = RDX(1)
   RDX(N_CELLS+1) = RDX(N_CELLS)

! Compute 1/dx_n for each node (dx_n is the distance from cell center to cell center)

   DO I=1,N_CELLS-1
      RDXN(I) = 2._EB/(X_S(I+1)-X_S(I-1))
   ENDDO
   RDXN(0)       = 1._EB/(X_S(1)-X_S(0))
   RDXN(N_CELLS) = 1._EB/(X_S(N_CELLS)-X_S(N_CELLS-1))

! Adjust 1/dx_n to r/dr for cylindrical case and r^2/dr for spaherical

   IF (GEOMETRY /= SURF_CARTESIAN) THEN
      DO I=0,N_CELLS
         R = THICKNESS-X_S(I)
         RDXN(I) = RDXN(I)*R**I_GRAD
      ENDDO
   ENDIF

END SUBROUTINE GET_WALL_NODE_WEIGHTS


SUBROUTINE GET_INTERPOLATION_WEIGHTS(N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS,N_LAYER_CELLS_NEW, &
            X_S,X_S_NEW,INT_WGT)

INTEGER, INTENT(IN)  :: N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS(N_LAYERS),N_LAYER_CELLS_NEW(N_LAYERS)
REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW)
REAL(EB), INTENT(OUT) :: INT_WGT(NWP_NEW,NWP)

REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,DX_NEW
INTEGER I, J, II, JJ, I_BASE, J_BASE, J_OLD,N
II = 0
JJ = 0
I_BASE = 0
J_BASE = 0

INT_WGT = 0._EB
DO N = 1,N_LAYERS
   J_OLD = 1
   DO I = 1,N_LAYER_CELLS_NEW(N)
      II       = I_BASE + I
      XUP_NEW  = X_S_NEW(II)
      XLOW_NEW = X_S_NEW(II-1)
      DX_NEW   = XUP_NEW - XLOW_NEW 
      DO J = J_OLD,N_LAYER_CELLS(N)
         JJ = J_BASE + J
         XUP =  X_S(JJ)
         XLOW = X_S(JJ-1)
         INT_WGT(II,JJ) = (MIN(XUP,XUP_NEW)-MAX(XLOW,XLOW_NEW))/DX_NEW
         IF (XUP .GE. XUP_NEW) EXIT
      ENDDO
      J_OLD = J
   ENDDO
   I_BASE = I_BASE + N_LAYER_CELLS_NEW(N)
   J_BASE = J_BASE + N_LAYER_CELLS(N)
ENDDO

END SUBROUTINE GET_INTERPOLATION_WEIGHTS


SUBROUTINE INTERPOLATE_WALL_ARRAY(NWP,NWP_NEW,INT_WGT,ARR)
INTEGER, INTENT(IN)  :: NWP,NWP_NEW
REAL(EB), INTENT(IN) :: INT_WGT(NWP_NEW,NWP)
REAL(EB) ARR(NWP),TMP(NWP)

INTEGER I,J

TMP = ARR
ARR = 0._EB
DO I = 1,NWP_NEW
DO J = 1,NWP
   ARR(I) = ARR(I) + INT_WGT(I,J)*TMP(J)
ENDDO
ENDDO

END SUBROUTINE INTERPOLATE_WALL_ARRAY

END MODULE GEOMETRY_FUNCTIONS 

 
MODULE PHYSICAL_FUNCTIONS

! Functions for physical quantities

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
IMPLICIT NONE

CONTAINS

SUBROUTINE GET_F_C(Z_1,Z_2,Z_3,F,C,Z_F)
!Returns progress variables for Mixture Fraction functions for suppression and CO production
REAL(EB) :: Z_F,WGT,Z_3_MAX,ZZ
REAL(EB), INTENT(IN) :: Z_1,Z_2,Z_3
REAL(EB), INTENT(OUT) :: F,C
INTEGER :: IZ1,IZ2