func.f90

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

 
REAL(EB), INTENT(IN) :: RAMP_INPUT,TAU
REAL(EB):: RAMP_POSITION
INTEGER,INTENT(IN)   :: RAMP_INDEX

SELECT CASE(RAMP_INDEX)
   CASE(-2)
      EVALUATE_RAMP = MAX(TANH(RAMP_INPUT/TAU),0._EB)
   CASE(-1)
      EVALUATE_RAMP = MIN( (RAMP_INPUT/TAU)**2 , 1.0_EB )
   CASE( 0)
      EVALUATE_RAMP = 1._EB
   CASE(1:)
 
      IF (RAMPS(RAMP_INDEX)%DEVC_INDEX > 0) THEN
         IF (DEVICE(RAMPS(RAMP_INDEX)%DEVC_INDEX)%CURRENT_STATE) THEN
            EVALUATE_RAMP = RAMPS(RAMP_INDEX)%VALUE
         ELSE
            RAMP_POSITION = MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,RAMP_INPUT - RAMPS(RAMP_INDEX)%T_MIN))
            EVALUATE_RAMP = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(NINT(RAMP_POSITION/RAMPS(RAMP_INDEX)%DT))
         ENDIF
      ELSEIF(RAMPS(RAMP_INDEX)%CTRL_INDEX > 0) THEN
         IF (CONTROL(RAMPS(RAMP_INDEX)%CTRL_INDEX)%CURRENT_STATE) THEN
            EVALUATE_RAMP = RAMPS(RAMP_INDEX)%VALUE
         ELSE
            RAMP_POSITION = MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,RAMP_INPUT - RAMPS(RAMP_INDEX)%T_MIN))
            EVALUATE_RAMP = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(NINT(RAMP_POSITION/RAMPS(RAMP_INDEX)%DT))
         ENDIF
      ELSE
         RAMP_POSITION = MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,RAMP_INPUT - RAMPS(RAMP_INDEX)%T_MIN))
         EVALUATE_RAMP = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(NINT(RAMP_POSITION/RAMPS(RAMP_INDEX)%DT))
      ENDIF
      RAMPS(RAMP_INDEX)%VALUE = EVALUATE_RAMP
END SELECT
 
END FUNCTION EVALUATE_RAMP

 
REAL(EB) FUNCTION ERFC(X)
 
! Complimentary ERF function
 
REAL(EB), INTENT(IN) :: X
REAL(EB) ERFCS(13), ERFCCS(24), ERC2CS(23),XSML,XMAX,SQEPS,SQRTPI,Y
DATA ERFCS( 1) /   -.049046121234691808_EB /
DATA ERFCS( 2) /   -.14226120510371364_EB /
DATA ERFCS( 3) /    .010035582187599796_EB /
DATA ERFCS( 4) /   -.000576876469976748_EB /
DATA ERFCS( 5) /    .000027419931252196_EB /
DATA ERFCS( 6) /   -.000001104317550734_EB /
DATA ERFCS( 7) /    .000000038488755420_EB /
DATA ERFCS( 8) /   -.000000001180858253_EB /
DATA ERFCS( 9) /    .000000000032334215_EB /
DATA ERFCS(10) /   -.000000000000799101_EB /
DATA ERFCS(11) /    .000000000000017990_EB /
DATA ERFCS(12) /   -.000000000000000371_EB /
DATA ERFCS(13) /    .000000000000000007_EB /
DATA ERC2CS( 1) /   -.069601346602309501_EB /
DATA ERC2CS( 2) /   -.041101339362620893_EB /
DATA ERC2CS( 3) /    .003914495866689626_EB /
DATA ERC2CS( 4) /   -.000490639565054897_EB /
DATA ERC2CS( 5) /    .000071574790013770_EB /
DATA ERC2CS( 6) /   -.000011530716341312_EB /
DATA ERC2CS( 7) /    .000001994670590201_EB /
DATA ERC2CS( 8) /   -.000000364266647159_EB /
DATA ERC2CS( 9) /    .000000069443726100_EB /
DATA ERC2CS(10) /   -.000000013712209021_EB /
DATA ERC2CS(11) /    .000000002788389661_EB /
DATA ERC2CS(12) /   -.000000000581416472_EB /
DATA ERC2CS(13) /    .000000000123892049_EB /
DATA ERC2CS(14) /   -.000000000026906391_EB /
DATA ERC2CS(15) /    .000000000005942614_EB /
DATA ERC2CS(16) /   -.000000000001332386_EB /
DATA ERC2CS(17) /    .000000000000302804_EB /
DATA ERC2CS(18) /   -.000000000000069666_EB /
DATA ERC2CS(19) /    .000000000000016208_EB /
DATA ERC2CS(20) /   -.000000000000003809_EB /
DATA ERC2CS(21) /    .000000000000000904_EB /
DATA ERC2CS(22) /   -.000000000000000216_EB /
DATA ERC2CS(23) /    .000000000000000052_EB /
DATA ERFCCS( 1) /     0.0715179310202925_EB /
DATA ERFCCS( 2) /   -.026532434337606719_EB /
DATA ERFCCS( 3) /    .001711153977920853_EB /
DATA ERFCCS( 4) /   -.000163751663458512_EB /
DATA ERFCCS( 5) /    .000019871293500549_EB /
DATA ERFCCS( 6) /   -.000002843712412769_EB /
DATA ERFCCS( 7) /    .000000460616130901_EB /
DATA ERFCCS( 8) /   -.000000082277530261_EB /
DATA ERFCCS( 9) /    .000000015921418724_EB /
DATA ERFCCS(10) /   -.000000003295071356_EB /
DATA ERFCCS(11) /    .000000000722343973_EB /
DATA ERFCCS(12) /   -.000000000166485584_EB /
DATA ERFCCS(13) /    .000000000040103931_EB /
DATA ERFCCS(14) /   -.000000000010048164_EB /
DATA ERFCCS(15) /    .000000000002608272_EB /
DATA ERFCCS(16) /   -.000000000000699105_EB /
DATA ERFCCS(17) /    .000000000000192946_EB /
DATA ERFCCS(18) /   -.000000000000054704_EB /
DATA ERFCCS(19) /    .000000000000015901_EB /
DATA ERFCCS(20) /   -.000000000000004729_EB /
DATA ERFCCS(21) /    .000000000000001432_EB /
DATA ERFCCS(22) /   -.000000000000000439_EB /
DATA ERFCCS(23) /    .000000000000000138_EB /
DATA ERFCCS(24) /   -.000000000000000048_EB /
DATA SQRTPI /1.7724538509055160_EB/
 
XSML = -200._EB
XMAX = 200._EB
SQEPS = 0.001_EB
 
IF (X<=XSML) THEN
   ERFC = 2._EB
ELSE
 
   IF (X<=XMAX) THEN
      Y = ABS(X)
      IF (Y<=1.0_EB) THEN  ! ERFC(X) = 1.0 - ERF(X) FOR -1._EB <= X <= 1.
         IF (Y<SQEPS)  ERFC = 1.0_EB - 2.0_EB*X/SQRTPI
         IF (Y>=SQEPS) ERFC = 1.0_EB - X*(1.0_EB + CSEVL (2._EB*X*X-1._EB, ERFCS, 10) )
      ELSE  ! ERFC(X) = 1.0 - ERF(X) FOR 1._EB < ABS(X) <= XMAX
         Y = Y*Y
         IF (Y<=4._EB) ERFC = EXP(-Y)/ABS(X) * (0.5_EB + CSEVL ((8._EB/Y-5._EB)/3._EB,ERC2CS, 10) )
         IF (Y>4._EB) ERFC = EXP(-Y)/ABS(X) * (0.5_EB + CSEVL (8._EB/Y-1._EB,ERFCCS, 10) )
         IF (X<0._EB) ERFC = 2.0_EB - ERFC
      ENDIF
   ELSE
      ERFC = 0._EB
   ENDIF
ENDIF
RETURN
 
END FUNCTION ERFC
 
 
SUBROUTINE GAUSSJ(A,N,NP,B,M,MP,IERROR)
 
! Solve a linear system of equations with Gauss-Jordon elimination
! Source: Press et al. "Numerical Recipes"
 
INTEGER :: M,MP,N,NP,I,ICOL,IROW,J,K,L,LL,INDXC(NP),INDXR(NP),IPIV(NP)
REAL(EB) :: A(NP,NP),B(NP,MP),BIG,DUM,PIVINV
INTEGER, INTENT(OUT) :: IERROR
 
IERROR = 0
IPIV(1:N) = 0
 
DO I=1,N
   BIG = 0._EB
   DO J=1,N
      IF (IPIV(J)/=1) THEN
         DO K=1,N
            IF (IPIV(K)==0) THEN
               IF (ABS(A(J,K))>=BIG) THEN
                  BIG = ABS(A(J,K))
                  IROW = J
                  ICOL = K
               ENDIF
            ELSE IF (IPIV(K)>1) THEN
               IERROR = 103   ! Singular matrix in gaussj
               RETURN
            ENDIF
         ENDDO
      ENDIF
   ENDDO
   IPIV(ICOL) = IPIV(ICOL) + 1
   IF (IROW/=ICOL) THEN
      DO L=1,N
         DUM = A(IROW,L)
         A(IROW,L) = A(ICOL,L)
         A(ICOL,L) = DUM
      ENDDO
      DO L=1,M
         DUM = B(IROW,L)
         B(IROW,L) = B(ICOL,L)
         B(ICOL,L) = DUM
      ENDDO
   ENDIF
   INDXR(I) = IROW
   INDXC(I) = ICOL
   IF (A(ICOL,ICOL)==0._EB) THEN
      IERROR = 103  ! Singular matrix in gaussj
      RETURN
      ENDIF
   PIVINV = 1._EB/A(ICOL,ICOL)
   A(ICOL,ICOL) = 1.
   A(ICOL,1:N) = A(ICOL,1:N) * PIVINV
   B(ICOL,1:M) = B(ICOL,1:M) * PIVINV
   DO LL=1,N
      IF (LL/=ICOL) THEN
         DUM = A(LL,ICOL)
         A(LL,ICOL) = 0._EB
         A(LL,1:N) = A(LL,1:N) - A(ICOL,1:N)*DUM
         B(LL,1:M) = B(LL,1:M) - B(ICOL,1:M)*DUM
      ENDIF
   ENDDO
ENDDO
DO L=N,1,-1
   IF (INDXR(L)/=INDXC(L)) THEN
      DO K=1,N
         DUM = A(K,INDXR(L))
         A(K,INDXR(L)) = A(K,INDXC(L))
         A(K,INDXC(L)) = DUM
      ENDDO
   ENDIF
ENDDO
 
END SUBROUTINE GAUSSJ
 
 
REAL(EB) FUNCTION CSEVL(X,CS,N)
 
REAL(EB), INTENT(IN) :: X
REAL(EB) CS(:),B1,B0,TWOX,B2
INTEGER NI,N,I
 
B1=0._EB
B0=0._EB
TWOX=2._EB*X
DO I=1,N
B2=B1
B1=B0
NI=N+1-I
B0=TWOX*B1-B2+CS(NI)
ENDDO
 
CSEVL = 0.5_EB*(B0-B2)
 
END FUNCTION CSEVL

INTEGER(2) FUNCTION TWO_BYTE_REAL(REAL_IN)
REAL(FB),INTENT(IN) :: REAL_IN
INTEGER(2) EXP,TEMP,I

IF (ABS(REAL_IN) <= 1.E-17_FB) THEN
   TWO_BYTE_REAL = 0
ELSEIF (ABS(REAL_IN) >= 1.E+16_FB) THEN
   DO I=0,14
      TWO_BYTE_REAL=IBSET(TWO_BYTE_REAL,I)
   ENDDO
ELSE
   EXP = FLOOR(LOG10(ABS(REAL_IN)))+1
   TEMP = ABS(REAL_IN * 10**(-REAL(EXP,FB)+3))
   EXP = EXP + 15
   TWO_BYTE_REAL = EXP * 2**10
   TWO_BYTE_REAL = TWO_BYTE_REAL + TEMP
ENDIF
IF (REAL_IN < 0._FB) TWO_BYTE_REAL = IBSET(TWO_BYTE_REAL,15)

END FUNCTION TWO_BYTE_REAL

INTEGER FUNCTION RLE_COMPRESSION(QIN,NQIN,QMIN,QMAX,COUT)

! Compress the array QIN(1:NQIN) by first mapping to one byte integers and then
! using Run-Length Encoding to convert runs of common integers IIIIIIII TO  #In
! where #=255 is a marker character to distinguish between literal characters 
! and runs of n repeats

INTEGER, INTENT(IN) :: NQIN
REAL, INTENT(IN), DIMENSION(NQIN) :: QIN
REAL, INTENT(IN) :: QMIN, QMAX
CHARACTER(1), DIMENSION(NQIN) :: COUT
CHARACTER(1), PARAMETER :: MARK=CHAR(255)
CHARACTER(1) :: THISCHAR, LASTCHAR
INTEGER :: IIN,IOUT,NREPEATS,IQVAL

IF (QMAX<=QMIN .OR. NQIN<=0) THEN
   RLE_COMPRESSION = 0
   RETURN
ENDIF

IOUT=1
LASTCHAR=MARK

DO IIN = 1, NQIN
   IQVAL = 254*(QIN(IIN) - QMIN)/(QMAX-QMIN)
   IF (IQVAL<0)   IQVAL=0
   IF (IQVAL>254) IQVAL=254
   THISCHAR =CHAR(IQVAL)

   IF (THISCHAR == LASTCHAR) THEN
      NREPEATS = NREPEATS + 1
   ELSE
      NREPEATS = 1
   ENDIF

   IF (NREPEATS>=1 .AND. NREPEATS<=3) THEN
      COUT(IOUT) = THISCHAR
      LASTCHAR=THISCHAR
   ELSEIF (NREPEATS>=4) THEN
      IF (NREPEATS==4) THEN
         IOUT = IOUT - 3
         COUT(IOUT) = MARK
         IOUT = IOUT + 1
         COUT(IOUT) = THISCHAR
         IOUT = IOUT + 1
      ENDIF
      IF (NREPEATS/=4) IOUT = IOUT - 1
      COUT(IOUT) = CHAR(NREPEATS)
      IF (NREPEATS==254) THEN
          NREPEATS=1
          LASTCHAR=MARK
      ENDIF
   ENDIF
   IOUT = IOUT + 1
ENDDO

RLE_COMPRESSION = IOUT - 1

END FUNCTION RLE_COMPRESSION

SUBROUTINE INTERPOLATE1D(X,Y,XI,ANS)
REAL(EB), INTENT(IN), DIMENSION(:) :: X, Y
REAL(EB), INTENT(IN) :: XI
REAL(EB), INTENT(OUT) :: ANS
INTEGER I, UX,LX

UX = UBOUND(X,1)
LX = LBOUND(X,1)

IF (XI <= X(LX)) THEN
  ANS = Y(LX)
ELSEIF (XI >= X(UX)) THEN
  ANS = Y(UX)
ELSE
  L1: DO I=LX,UX-1
    IF (XI -X(I) == 0._EB) THEN
      ANS = Y(I)
      EXIT L1
    ELSEIF (X(I+1)>XI) THEN
      ANS = Y(I)+(XI-X(I))/(X(I+1)-X(I)) * (Y(I+1)-Y(I))
      EXIT L1
    ENDIF
  ENDDO L1
ENDIF

END SUBROUTINE INTERPOLATE1D

END MODULE MATH_FUNCTIONS


MODULE TRAN 
 
! Coordinate transformation functions
 
USE PRECISION_PARAMETERS
IMPLICIT NONE
TYPE TRAN_TYPE
   REAL(EB), POINTER, DIMENSION(:,:) :: C1,C2,C3,CCSTORE,PCSTORE
   INTEGER, POINTER, DIMENSION(:,:) :: IDERIVSTORE
   INTEGER NOC(3),ITRAN(3),NOCMAX
END TYPE TRAN_TYPE
TYPE (TRAN_TYPE), ALLOCATABLE, TARGET, DIMENSION(:) :: TRANS
 
 
CONTAINS
 
 
REAL(EB) FUNCTION G(X,IC,NM)
 
! Coordinate transformation function
 
REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T
 
T => TRANS(NM)
 
N = T%NOC(IC)
IF (N==0) THEN
   G = X
   RETURN
ENDIF
 
SELECT CASE(T%ITRAN(IC))
   CASE(1)
      G = 0._EB
      DO I=1,N+1
         G = G + T%C1(I,IC)*X**I
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      G = T%C2(II-1,IC) + T%C3(II,IC)*(X-T%C1(II-1,IC))
END SELECT
 
END FUNCTION G
 
 
REAL(EB) FUNCTION GP(X,IC,NM)
 
! Derivative of the coordinate transformation function
 
REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T
 
T => TRANS(NM)
N =  T%NOC(IC)
IF (N==0) THEN
   GP = 1._EB
   RETURN
ENDIF
 
SELECT CASE(T%ITRAN(IC)) 
   CASE(1)
      GP = 0._EB
      DO I=1,N+1
         GP = GP + I*T%C1(I,IC)*X**(I-1)
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GP = T%C3(II,IC)
END SELECT
 
END FUNCTION GP
 
 
REAL(EB) FUNCTION GINV(Z,IC,NM)
 
! Inverse of the coordinate transformation function
 
REAL(EB) :: GF
INTEGER :: N,IT,II,I
REAL(EB), INTENT(IN) :: Z
INTEGER, INTENT(IN)  :: IC,NM
TYPE (TRAN_TYPE), POINTER :: T
 
T => TRANS(NM)
GINV = Z
N = T%NOC(IC)
IF (N==0) RETURN
 
SELECT CASE(T%ITRAN(IC))
   CASE(1)
      LOOP1: DO IT=1,10
         GF = G(GINV,IC,NM)-Z
         IF (ABS(GF)<0.002_EB) EXIT LOOP1
         GINV = GINV - GF/GP(GINV,IC,NM)
      ENDDO LOOP1
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (Z<=T%C2(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GINV = T%C1(II-1,IC) + (Z-T%C2(II-1,IC))/T%C3(II,IC)
END SELECT
 
END FUNCTION GINV
 
 
END MODULE TRAN