pois.f90

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

DATARY=.TRUE.
SCALE=1._EB 
IFWRD = 1
100 CONTINUE

IF (N/=1) THEN
  TPOSE=.FALSE.
  IF (IFWRD==2) TPOSE=.TRUE.
  
!                               TRANSFORM IN Z
  IF (DATARY) THEN
    CALL FSH26S(IGRID,IFWRD,NP,L,N,M,LDIMFT,F,FT,CFZ,WSAVEZ)
    DATARY=OUTARY

  ELSE
    CALL FSH26S(IGRID,IFWRD,NP,L,N,M,LDIMFT,FT,F,CFZ,WSAVEZ)
    DATARY=.NOT.OUTARY
  END IF
  
END IF
SELECT CASE(IFWRD)
CASE(1)   ; GO TO 490
CASE(2)   ; GO TO 510
END SELECT

490 CONTINUE
IF (M/=1) THEN
  TPOSE=.TRUE.
  IF (IFWRD==2) TPOSE=.FALSE.
  
!                               TRANSFORM Y
  
  
  IF (DATARY) THEN
    CALL FSH26S(IGRID,IFWRD,MP,L,M,N,LDIMFT,F,FT,CFY,WSAVEY)
    DATARY=OUTARY
  ELSE
    CALL FSH26S(IGRID,IFWRD,MP,L,M,N,LDIMFT,FT,F,CFY,WSAVEY)
    DATARY=.NOT.OUTARY
  END IF
  
END IF

SELECT CASE(IFWRD)
CASE(1)   ; GO TO 285
CASE(2)   ; GO TO 100
END SELECT
285 CONTINUE

IF (L>1) THEN
  
!                               SOLVE TRIDIAGONAL SYSTEMS IN X THAT WERE
!                               PREVIOUSLY FACTORED IN FSH01S
  
!                               CALL VECTORIZED TRIDIAGONAL SOLVER
  
  DATASW=.FALSE.
  IF (NP==1) DATASW=.NOT.DATASW
  IF (MP==1) DATASW=.NOT.DATASW
  IF (DATARY) THEN
    IF (DATASW) THEN
      CALL FSH06S(L,LP,M*N,LDIMFC,LDIMFT,SCALE,A,C,F,FT,FCTRD)
      DATARY=.FALSE.
    ELSE
      CALL FSH06S(L,LP,M*N,LDIMFC,LDIMFT,SCALE,A,C,F,F,FCTRD)
    END IF
  ELSE
    IF (DATASW) THEN
      CALL FSH06S(L,LP,M*N,LDIMFC,LDIMFT,SCALE,A,C,FT,F,FCTRD)
      DATARY=.TRUE.
    ELSE
      CALL FSH06S(L,LP,M*N,LDIMFC,LDIMFT,SCALE,A,C,FT,FT,FCTRD)
    END IF
  END IF
  
END IF

IFWRD = 2

GO TO 490

510 CONTINUE

IF (.NOT.DATARY) THEN
  DO  K=1,N
    DO  J=1,M
      DO  I=1,L
        F(I,J,K)=FT(I,J,K)
      END DO
    END DO
  END DO
END IF
RETURN
END SUBROUTINE FSH03S


SUBROUTINE FSH04S(L,M,N,LDIMF,MDIMF,LDIMG,F,G)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+



INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMF
INTEGER                   :: MDIMF
INTEGER                   :: LDIMG
REAL(EB)   F(LDIMF,MDIMF, N), G(LDIMG,M,N)

!                              THIS SUBROUTINE PACKS THE SUB-ARRAY
!                              F(I,J,K), I=1,...,L, J=1,...,M, K=1,...,N
!                              INTO THE ARRAY G.

DO  K = 1,N
  DO  J = 1,M
    DO  I = 1,L
      G(I,J,K) = F(I,J,K)
    END DO
  END DO
END DO
IF (LDIMG>L) THEN
  DO  K=1,N
    DO  J=1,M
      G(LDIMG,J,K)=0._EB 
    END DO
  END DO
END IF
RETURN
END SUBROUTINE FSH04S


SUBROUTINE FSH05S(L,M,N,LDIMF,MDIMF,LDIMG,F,G)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+



INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMF
INTEGER                   :: MDIMF
INTEGER                   :: LDIMG
REAL(EB)   G(LDIMG,M,N)
REAL(EB)   F(LDIMF,MDIMF,N)

!                               THIS SUBROUTINE EXPANDS THE ARRAY G OF
!                               DIMENSION L X M X N INTO THE ARRAY F OF
!                               DIMENSION LDIMF X MDIMF X N.

DO  K = 1,N
  DO  J = 1,M
    DO  I = 1,L
      F(I,J,K) = G(I,J,K)
    END DO
  END DO
END DO
RETURN
END SUBROUTINE FSH05S


SUBROUTINE FSH06S(L,LP,M,LDIMFC,LDIMFT,SCALE,A,C,F,FT,FCTRD)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+



INTEGER                       :: L
INTEGER                   :: LP
INTEGER                       :: M
INTEGER                   :: LDIMFC
INTEGER                   :: LDIMFT
REAL(EB) :: SCALE
REAL(EB)   FT(LDIMFT,M)
REAL(EB)  FCTRD(LDIMFC,M)
REAL(EB)   A(L),C(L),F(LDIMFT,M)

!                               THIS SUBROUTINE SOLVES M TRIDIAGONAL
!                               SYSTEMS OF ORDER L THAT HAVE
!                               FACTORIZATION STORED IN ARRAY FCTRD AND
!                               RIGHT SIDES IN FT
SCALE=SCALE**2
IF (LP==1) THEN
  LH=(L+1)/2
  LQ=(LH-1)/2
  DO  I=1,LQ
    DO  J=1,M
      F(I,J)=F(I,J)-F(2*LH-I,J)
      F(LH-I,J)=F(LH-I,J)-F(LH+I,J)
      F(LH+I,J)=F(LH-I,J)+2._EB*F(LH+I,J)
      F(2*LH-I,J)=F(I,J)+2._EB*F(2*LH-I,J)
    END DO
    CALL SSWAP(M,F(I,1),LDIMFT,F(LH-I,1),LDIMFT)
  END DO
  DO  J=1,M
    F(LH,J)=2._EB*F(LH,J)
  END DO
IF (MOD(L,2)==0)  THEN
  DO  J=1,M
    F(L,J)=2._EB*F(L,J)
  END DO
END IF
IF (MOD(LH,2)==0) THEN
  DO  J=1,M
    F(  LH/2,J)=F(LH/2,J)-F(3*LH/2,J)
    F(3*LH/2,J)=F(LH/2,J)+2._EB*F(3*LH/2,J)
  END DO
END IF
SCALE=0.5_EB*SCALE
END IF
!                               FORWARD SUBSTITUTION

DO  J = 1,M
  FT(1,J) = SCALE*F(1,J)*FCTRD(1,J)
  DO  I = 2,L
    FT(I,J) = (SCALE*F(I,J)-A(I)*FT(I-1,J))*FCTRD(I,J)
  END DO
END DO

!                               BACKWARD SUBSTITUTION
DO  J = 1,M
  DO  I = L - 1,1,-1
    FT(I,J) = FT(I,J) - C(I)*FCTRD(I,J)*FT(I+1,J)
  END DO
END DO

IF (LP==1) THEN
  
  DO  I=1,LQ
    DO  J=1,M
      FT(I,J)=FT(LH+I,J)+FT(I,J)
      FT(LH-I,J)=FT(2*LH-I,J)+FT(LH-I,J)
      FT(2*LH-I,J)=2._EB*FT(2*LH-I,J)-FT(LH-I,J)
      FT(  LH+I,J)=2._EB*FT(LH+I,J)-FT(I,J)
    END DO
    CALL SSWAP(M,FT(I,1),LDIMFT,FT(LH-I,1),LDIMFT)
  END DO
  IF (MOD(LH,2)==0) THEN
    DO  J=1,M
      FT(  LH/2,J)=FT(3*LH/2,J)+FT(LH/2,J)
      FT(3*LH/2,J)=2._EB*FT(3*LH/2,J)-FT(LH/2,J)
    END DO
  END IF
  
END IF

RETURN
END SUBROUTINE FSH06S


FUNCTION FSH19S()

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+


!                               THIS FUNCTION FURNISHES PI TO MACHINE
!                               PREC

!                               PI=3.14159265358979323846264338327950288
REAL(EB) FSH19S
FSH19S = 4._EB* ATAN(1._EB) !3.14159265358979323846264338327950288_EB
RETURN
END FUNCTION FSH19S


FUNCTION FSH20S()

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+

!                               THIS FUNCTION FURNISHES THE MACHINE
!                               PREC EPS--THE SMALLEST POSITIVE
!                               MACHINE NUMBER SATISFYING

!                               FL(1._EB+EPS) > 1.

!     FSH20S = 2.4E-7_EB
REAL(EB) FSH20S
FSH20S = SPACING (1._EB)!4.E-15_EB

RETURN
END FUNCTION FSH20S


SUBROUTINE FSH26S(IGRID,IFWRD,MP,L,M,N,LDIMFT,F,FT,CFY,WSAVEY)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+


INTEGER                   :: IGRID
INTEGER                   :: IFWRD
INTEGER                   :: MP
INTEGER                   :: L
INTEGER                   :: M
INTEGER                   :: N
INTEGER                   :: LDIMFT
REAL(EB) FT(LDIMFT,M,N)
REAL(EB) CFY(4*M)
REAL(EB) WSAVEY(M+16)
REAL(EB) F(LDIMFT,M,N)

SELECT CASE(IGRID)
CASE(1)   ; GO TO 110
CASE(2)   ; GO TO 120
END SELECT

110     CONTINUE
SELECT CASE(MP)
CASE(1)   ; GO TO 130
CASE(2)   ; GO TO 160
CASE(3)   ; GO TO 170
CASE(4)   ; GO TO 220
CASE(5)   ; GO TO 230
END SELECT
120     CONTINUE
SELECT CASE(MP)
CASE(1)   ; GO TO 130
CASE(2)   ; GO TO 180
CASE(3)   ; GO TO 210
CASE(4)   ; GO TO 240
CASE(5)   ; GO TO 270
END SELECT
130     CONTINUE
SELECT CASE(IFWRD)
CASE(1)   ; GO TO 140
CASE(2)   ; GO TO 150
END SELECT
140     CONTINUE

CALL VSRFTF(F,L,M,N,LDIMFT,FT,WSAVEY)
GO TO 280
150     CONTINUE
CALL VSRFTB(F,L,M,N,LDIMFT,FT,WSAVEY)
GO TO 280
160     CONTINUE
CALL VSINT(F,L,M,N,LDIMFT,FT,CFY,WSAVEY)
GO TO 280
170     CONTINUE
SELECT CASE(IFWRD)
CASE(1)   ; GO TO 200
CASE(2)   ; GO TO 190
END SELECT
180     CONTINUE
SELECT CASE(IFWRD)
CASE(1)   ; GO TO 190
CASE(2)   ; GO TO 200
END SELECT
190     CONTINUE
CALL VSSINF(F,L,M,N,LDIMFT,FT,CFY(1),CFY(M+1),WSAVEY)
GO TO 280
200     CONTINUE
CALL VSSINB(F,L,M,N,LDIMFT,FT,CFY(1),CFY(M+1),WSAVEY)
GO TO 280
210     CONTINUE
CALL VSSINQ(F,L,M,N,LDIMFT,FT,CFY(1),CFY(M+1), CFY(2*M+1),CFY(3*M+1),WSAVEY)
GO TO 280
220     CONTINUE
CALL VCOST(F,L,M,N,LDIMFT,FT,CFY,WSAVEY)
GO TO 280
230     CONTINUE
SELECT CASE(IFWRD)
CASE(1)   ; GO TO 260
CASE(2)   ; GO TO 250
END SELECT
240     CONTINUE
SELECT CASE(IFWRD)
CASE(1)   ; GO TO 250
CASE(2)   ; GO TO 260
END SELECT
250     CONTINUE
CALL VSCOSF(F,L,M,N,LDIMFT,FT,CFY(1),CFY(M+1),WSAVEY)
GO TO 280
260     CONTINUE
CALL VSCOSB(F,L,M,N,LDIMFT,FT,CFY(1),CFY(M+1),WSAVEY)
GO TO 280
270     CONTINUE
CALL VSCOSQ(F,L,M,N,LDIMFT,FT,CFY(1),CFY(M+1), CFY(2*M+1),CFY(3*M+1),WSAVEY)
280     CONTINUE

RETURN
END SUBROUTINE FSH26S


SUBROUTINE VCOST(X,L,M,N,LDIMX,XT,C,WSAVE)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989

INTEGER                   :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                       :: LDIMX
REAL(EB) XT(LDIMX*N,M)
REAL(EB) C(M)
REAL(EB) WSAVE(M+15)
REAL(EB) X(LDIMX*N,M)

!     TPOSE = .TRUE. IF TRANSFORMING SECOND INDEX (SO TRANSPOSE)
!           = .FALSE. IF TRANSFORMING THIRD INDEX

MM1 = M-1
MS2 = M/2
SCALE=SCALE*SQRT(0.5_EB)
IF (TPOSE) THEN
  CALL VCOST1(L,M,N,LDIMX,MS2,C,XT,X)
ELSE
  DO  I=1,LDIMX*N
    XT(I,M)=  X(I,1)-X(I,M)
    XT(I,1) = X(I,1)+X(I,M)
  END DO
DO  J=2,MS2
  JC = M+1-J
  DO  I=1,LDIMX*N
    XT(I,M) = XT(I,M)+C(JC)*(X(I,J)-X(I,JC))
    XT(I,J) =  (X(I,J)+X(I,JC))-(C(J)*(X(I,J)-X(I,JC)))
    XT(I,JC) = (X(I,J)+X(I,JC))+(C(J)*(X(I,J)-X(I,JC)))
  END DO
END DO
IF (MOD(M,2) == 0) GO TO 404
DO  I=1,LDIMX*N
  XT(I,MS2+1) = 2._EB*X(I,MS2+1)
END DO
404    CONTINUE
END IF
IF (M>3) THEN
  450    CALL VRFFTF (LDIMX*N,MM1,XT,LDIMX*N,X,WSAVE)
  IF (OUTARY) THEN
    CALL VCOSTA(M,N*LDIMX,MM1,XT(1,M),X,XT)
  ELSE
    CALL VCOSTA(M,N*LDIMX,MM1,XT(1,M),XT,X)
  END IF
ELSE IF (M==2) THEN
  OUTARY=.FALSE.
ELSE
  DO  I=1,LDIMX*N
    X(I,2)=XT(I,3)
    X(I,1)=XT(I,1)+XT(I,2)
    X(I,3)=XT(I,1)-XT(I,2)
  END DO
OUTARY=.TRUE.
SCALE=SCALE*SQRT(0.5_EB)
END IF
RETURN
END SUBROUTINE VCOST


SUBROUTINE VCOST1(L,M,N,LDIMX,MS2,C,XT,X)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMX
INTEGER                   :: MS2
REAL(EB)   C(*)
REAL(EB)  XT(LDIMX,N,M)
REAL(EB)   X(LDIMX,M,N)
DO  K=1,N
  DO  I=1,L
    XT(I,K,M)=X(I,1,K)-X(I,M,K)
    XT(I,K,1) = X(I,1,K)+X(I,M,K)
  END DO
END DO
DO  J=2,MS2
  JC = M+1-J
  DO  K=1,N
    DO  I=1,L
      XT(I,K,M) = XT(I,K,M)+C(JC)*(X(I,J,K)-X(I,JC,K))
      XT(I,K,J) =  (X(I,J,K)+X(I,JC,K))-(C(J)*(X(I,J,K)-X(I,JC,K)))
      XT(I,K,JC) = (X(I,J,K)+X(I,JC,K))+(C(J)*(X(I,J,K)-X(I,JC,K)))
    END DO
  END DO
END DO
IF (MOD(M,2) == 0) GO TO 404
DO  K=1,N