pois.f90

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

   END IF
      DO  J=2,M-1,2
         JBY2=J/2
         DO  I=1,LDIMF*N
            FT(I,J)   =  F(I,J+1)*C1(JBY2)+F(I,J)*C2(JBY2)
            FT(I,J+1) = -F(I,J+1)*C2(JBY2)+F(I,J)*C1(JBY2)
         END DO
      END DO
END IF

!     REAL(EB),PERIODIC SYNTHESIS

CALL VRFFTB(LDIMF*N,M,FT,LDIMF*N,F,WORK)

!     POSTPROCESSING

IF (OUTARY) THEN
  DO  J=1,M
    DO  I=1,LDIMF*N
      F(I,J)=C4(J)*FT(I,J)+C3(J)*FT(I,M+1-J)
    END DO
  END DO
ELSE
   DO  J=1,M
      DO  I=1,LDIMF*N
         FT(I,J)=C4(J)*F(I,J)+C3(J)*F(I,M+1-J)
      END DO
   END DO
END IF
RETURN
END SUBROUTINE VSCOSQ


SUBROUTINE VSCSB1(L,M,N,LDIMF,F,FT,C1,C2)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMF
REAL(EB)  FT(LDIMF,N,M)
REAL(EB)  C1(M)
REAL(EB)  C2(M)
REAL(EB)   F(LDIMF,M,N)
DO  K=1,N
  DO  I=1,L
    FT(I,K,1)=.5_EB*F(I,1,K)
  END DO
END DO
DO  K=1,N
  DO  J=2,M
    DO  I=1,L
      FT(I,K,J) =(C1(J)*F(I,J,K)+C2(J)*F(I,M-J+2,K))
    END DO
  END DO
END DO
RETURN
END SUBROUTINE VSCSB1


SUBROUTINE VSCSBA(M,LDIMF,FT,F)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: M
INTEGER                       :: LDIMF
REAL(EB)   F(LDIMF,M), FT(LDIMF,M)

DO  I=1,LDIMF
  F(I,1) = FT(I,1)
END DO
IF (MOD(M,2)==0) THEN
  DO  I=1,LDIMF
    F(I,M) = FT(I,M)
  END DO
END IF
DO  J=2,M-1,2
  DO  I=1,LDIMF
    F(I,J)=FT(I,J)-FT(I,J+1)
    F(I,J+1)=FT(I,J)+FT(I,J+1)
  END DO
END DO
RETURN
END SUBROUTINE VSCSBA


SUBROUTINE VSCSF1(L,M,N,LDIMF,F,FT)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMF
REAL(EB)   FT(LDIMF,N,M)
REAL(EB)   F(LDIMF,M,N)
DO  K=1,N
  DO  I=1,L
    FT(I,K,1)=F(I,1,K)
  END DO
END DO
IF (MOD(M,2)==0) THEN
  DO  K=1,N
    DO  I=1,L
      FT(I,K,M)=F(I,M,K)
    END DO
  END DO
END IF
DO  K=1,N
  DO  J=2,M-1,2
    DO  I=1,L
      FT(I,K,J)   =  0.5_EB*(F(I,J,K)+F(I,J+1,K))
      FT(I,K,J+1) = -0.5_EB*(F(I,J,K)-F(I,J+1,K))
    END DO
  END DO
END DO
RETURN
END SUBROUTINE VSCSF1


SUBROUTINE VSCSFA(M,LDIMF,F,FT,C1,C2)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: M
INTEGER                       :: LDIMF
REAL(EB)    FT(LDIMF,M)
REAL(EB)    C1(M)
REAL(EB)   F(LDIMF,M),  C2(M)

DO  J=2,M
  DO  I=1,LDIMF
    F(I,J)=C1(J)*FT(I,J)-C2(J)*FT(I,M+2-J)
  END DO
END DO
DO  I=1,LDIMF
  F(I,1)=FT(I,1)
END DO
RETURN
END SUBROUTINE VSCSFA


SUBROUTINE VSCSQ1(L,M,N,LDIMF,F,FT,C1,C2)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMF
REAL(EB)  FT(LDIMF,N,M)
REAL(EB)  C1(M)
REAL(EB)    C2(M)
REAL(EB)   F(LDIMF,M,N)
DO  K=1,N
  DO  I=1,L
    FT(I,K,1)=F(I,1,K)
  END DO
END DO
IF (MOD(M,2)==0) THEN
  DO  K=1,N
    DO  I=1,L
      FT(I,K,M)=-F(I,M,K)
    END DO
  END DO
END IF
DO  J=2,M-1,2
  JBY2=J/2
  DO  K=1,N
    DO  I=1,L
      FT(I,K,J)   =  F(I,J+1,K)*C1(JBY2)+F(I,J,K)*C2(JBY2)
      FT(I,K,J+1) = -F(I,J+1,K)*C2(JBY2)+F(I,J,K)*C1(JBY2)
    END DO
  END DO
END DO
RETURN
END SUBROUTINE VSCSQ1


SUBROUTINE VSCSQI(N,C1,C2,C3,C4,WSAVE)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989

!      ENTRY VSSNQI(N,C1,C2,C3,C4,WSAVE)


INTEGER                       :: N
REAL(EB)  C3(N)
REAL(EB)   C4(N)
REAL(EB)   WSAVE(N+15)
REAL(EB)   C1(N),C2(N)

PI=FSH19S()
DX=PI/N
SCALE=SQRT(.5_EB)

!     GENERATE A(I)+-B(I)

DO  I=1,(N-1)/2
  C=COS(I*DX)
  S=SIN(I*DX)
  C1(I)=.5_EB*(S+C)
  C2(I)=.5_EB*(C-S)
END DO

DX=PI/(2*N)
DO  I=1,N
  C=COS((I-.5_EB)*DX)
  S=SIN((I-.5_EB)*DX)
  C3(I)=SCALE*(C+S)
  C4(I)=SCALE*(C-S)
END DO

!     INITIALIZE VRFFTPK ROUTINE

CALL VRFFTI(N,WSAVE)
RETURN
END SUBROUTINE VSCSQI


SUBROUTINE VSINT(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+1)
REAL(EB)  C(M/2)
REAL(EB)  WSAVE(M+16)
REAL(EB)  X(LDIMX*N,M)

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

SQRT2I=SQRT(.5_EB)
MODM = MOD(M,2)
103 MP1 = M+1
MS2 = M/2
DO  I=1,LDIMX*N
  XT(I,1) = 0._EB 
  XT(I,M+1)=0._EB 
END DO
!                         ZERO OUT LAST PLANE BECAUSE FOR SINE
!                         IT DOESN'T GET DONE IN FSH02
IF (TPOSE) THEN
  CALL VSINT1(L,M,N,LDIMX,MODM,MS2,C,XT,X)
ELSE
  DO  J=1,MS2
    JC=M+1-J
    DO  I=1,LDIMX*N
      XT(I,J+1) = (X(I,J)-X(I,JC))+C(J)*(X(I,J)+X(I,JC))
      XT(I,JC+1) = C(J)*(X(I,J)+X(I,JC))-(X(I,J)-X(I,JC))
    END DO
  END DO
  IF (MODM == 0) GO TO 405
  DO  I=1,LDIMX*N
    XT(I,MS2+2) = 4._EB*X(I,MS2+1)
  END DO
405    CONTINUE
END IF
CALL VRFFTF (LDIMX*N,MP1,XT,LDIMX*N,X,WSAVE)
SCALE=SCALE*SQRT2I
IF (OUTARY) THEN
  CALL VSINTA(M,LDIMX*N,MODM,X,XT)
ELSE
  CALL VSINTA(M,LDIMX*N,MODM,XT,X)
END IF
RETURN
END SUBROUTINE VSINT


SUBROUTINE VSINT1(L,M,N,LDIMX,MODM,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                   :: MODM
INTEGER                   :: MS2
REAL(EB)   C(*)
REAL(EB)   X(LDIMX,M,N)
REAL(EB)   XT(LDIMX,N,M+1)
DO  J=1,MS2
  JC = M+1-J
  DO  K=1,N
    DO  I=1,L
      XT(I,K,J+1) = (X(I,J,K)-X(I,JC,K))+C(J)*(X(I,J,K)+X(I,JC,K))
      XT(I,K,JC+1) = C(J)*(X(I,J,K)+X(I,JC,K))-(X(I,J,K)-X(I,JC,K))
    END DO
  END DO
END DO
IF (MODM == 0) GO TO 405
DO  K=1,N
  DO  I=1,L
    XT(I,K,MS2+2) = 4._EB*X(I,MS2+1,K)
  END DO
END DO
405 CONTINUE
RETURN
END SUBROUTINE VSINT1


SUBROUTINE VSINTA(M,LDIMX,MODM,X,XT)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989



INTEGER                       :: M
INTEGER                       :: LDIMX
INTEGER                   :: MODM
REAL(EB)   X(LDIMX,M),  XT(LDIMX,M+1)

DO  I=1,LDIMX
  X(I,1) = 0.5_EB*XT(I,1)
END DO
DO  J=2,M-1,2
  DO  I=1,LDIMX
    X(I,J) = -XT(I,J+1)
    X(I,J+1) = X(I,J-1)+XT(I,J)
  END DO
END DO
IF (MODM /= 0) GO TO 900
DO  I=1,LDIMX
  X(I,M) = -XT(I,M+1)
END DO
900 CONTINUE
RETURN
END SUBROUTINE VSINTA


SUBROUTINE VSINTI(N,C,WSAVE)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989

!     *******   NOTE   ******

!     SPECIAL DOCUMENTATION FOR USE OF INTERNAL PARAMETER TPOSE.  IF
!     TPOSE = .TRUE., THEN ALL PREPROCESSING ROUTINES TRANSPOSE SECOND
!                     AND THIRD INDICES IN PREPARATION FOR TRANSFORM OF
!                     ORIGINAL SECOND INDEX.  IF
!     TPOSE = .FALSE., THEN ALL PREPROCESSING ROUTINES PREPARE ONLY FOR
!                      TRANSFORM OF THIRD INDEX.  TO USE IN THIS CASE,
!                      CALL ROUTINE WITH M AND N INTERCHANGED.



INTEGER                       :: N
REAL(EB)   C(N/2),WSAVE (N+16)

!                               INITIALIZE NOCOPY AND TPOSE TO DEFAULT
!                               VALUES

NOCOPY = .FALSE.
TPOSE = .TRUE.

IF (N <= 1) RETURN
PI=FSH19S()
NP1 = N+1
NS2 = N/2
DT = PI/REAL(NP1,EB)
DO  K=1,NS2
  C(K) = 2._EB*SIN(K*DT)
END DO
CALL VRFFTI (NP1,WSAVE)
RETURN
END SUBROUTINE VSINTI


SUBROUTINE VSRFTB(F,L,M,N,LDIMF,FT,WSAVE)

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

!     PACKAGE VFFTPAK, VERSION 1, JUNE 1989


!     *******   NOTE   ******

!     SPECIAL DOCUMENTATION FOR USE OF INTERNAL PARAMETER TPOSE.  IF
!     TPOSE = .TRUE., THEN ALL PREPROCESSING ROUTINES TRANSPOSE SECOND
!                     AND THIRD INDICES IN PREPARATION FOR TRANSFORM OF
!                     ORIGINAL SECOND INDEX.  IF
!     TPOSE = .FALSE., THEN ALL PREPROCESSING ROUTINES PREPARE ONLY FOR
!                      TRANSFORM OF THIRD INDEX.  TO USE IN THIS CASE,
!                      CALL ROUTINE WITH M AND N INTERCHANGED.


!     VSFFTPK, VERSION 2, JUNE 1988


INTEGER                       :: L
INTEGER                       :: M
INTEGER                       :: N
INTEGER                   :: LDIMF
REAL(EB)   FT(LDIMF,N,M)
REAL(EB)   WSAVE(M+15)
REAL(EB)   F(LDIMF,M,N)

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

IF (TPOSE) THEN
  
!        RE-ORDER INPUT
  
  DO  K=1,N
    DO  J=1,M
      DO  I=1,L
        FT(I,K,J)=F(I,J,K)
      END DO
    END DO
  END DO
  
!        REAL(EB), PERIODIC TRANSFORM
  
  CALL VRFFTB(LDIMF*N,M,FT,LDIMF*N,F,WSAVE)
  OUTARY=.NOT.OUTARY
  IF (.NOT.NOCOPY) THEN
    CALL VSRTB1(M,N*LDIMF,F,FT)
  END IF
ELSE
  CALL VRFFTB(LDIMF*N,M,F,LD