pois.f90

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

LDIMFC=L
IF (LDIMF>L .AND. MOD(L,2)==0) LDIMFC=L+1

IGRID = 2

CALL FSH00S(IGRID,LPEROD,L,MPEROD,M,NPEROD,N,LDIMFC, SCAL,A,B,C,D,SAVE)

RETURN
END SUBROUTINE S3CFIS


SUBROUTINE S3CFSS(LDIMF,MDIMF,F,SAVE,W)

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

INTEGER:: LDIMF, MDIMF
REAL(EB)   F(LDIMF,MDIMF,*), SAVE(-3:*), W(*)

!                               CHECK VALUE OF IERROR (=SAVE(1)).
!                               IF NON-ZERO, RETURN.

IF (SAVE(1)/=0._EB) RETURN

CALL FSH02S(LDIMF,MDIMF,F,SAVE,W)

RETURN
END SUBROUTINE S3CFSS


SUBROUTINE FSH00S(IGRID,LPEROD,L,MPEROD,M,NPEROD,N,LDIMFC,C2, A,B,C,D,SAVE)

! +--------------------------------------------------------------------+
! |                                                                    |
! |                       COPYRIGHT (C) 1989 BY                        |
! |                          ROLAND A. SWEET                           |
! |                        ALL RIGHTS RESERVED                         |
! |                                                                    |
! +--------------------------------------------------------------------+
REAL(EB)   a(l),b(l),c(l),d(l),save(-3:*), C2
INTEGER:: IGRID, LPEROD, L,  MPEROD,  M, NPEROD, N, LDIMFC
!                               THIS SUBROUTINE INITIALIZES FFT SOLVER

!                               ALLOCATE SAVE ARRAY

IA = 12
IB = IA + L
IC = IB + L
ICFY = IC + L

IF (IGRID==1) THEN
  ICFZ = ICFY + 2*M
  IFCTRD = ICFZ + 2*N
ELSE
  ICFZ = ICFY + 4*M
  IFCTRD = ICFZ + 4*N
END IF

IWSY = IFCTRD + LDIMFC*M*N
IWSZ = IWSY + M + 16
IZRT = IWSZ + N + 16

!                               COPY COEFFICIENT ARRAYS A,B, AND C INTO
!                               SAVE ARRAY.  A COPY OF B IS MADE BECAUSE
!                               IN THE NEXT LEVEL ROUTINE, BOUNDARY
!                               ELEMENTS OF B MAY BE CHANGED.

DO  I = 0,L - 1
  SAVE(IA+I) = A(I+1)
  SAVE(IB+I) = B(I+1)
  SAVE(IC+I) = C(I+1)
END DO
LP = LPEROD + 1
MP = MPEROD + 1
NP = NPEROD + 1

!                               CALL LOWER LEVEL INITIALIZATION ROUTINE

CALL FSH01S(IGRID,L,LP,M,MP,D,N,NP,LDIMFC,C2,SAVE(IA),SAVE(IB),  &
    SAVE(IC),SAVE(ICFY),SAVE(ICFZ),SAVE(IFCTRD), SAVE(IWSY),SAVE(IWSZ),SAVE(IZRT))

!                               SAVE PARAMETERS FOR SUBROUTINE SOLVER

SAVE(2) = L
SAVE(3) = LP
SAVE(4) = M
SAVE(5) = MP
SAVE(6) = N
SAVE(7) = NP
SAVE(8) = ICFZ
SAVE(9) = IWSZ
SAVE(10) = IZRT
SAVE(11) = IGRID

RETURN
END SUBROUTINE FSH00S


SUBROUTINE FSH01S(IGRID,L,LP,M,MP,D,N,NP,LDIMFC,C2,A,B,C,CFY,CFZ,  &
    FCTRD,WSAVEY,WSAVEZ,ZRT)

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

REAL(EB) a(l),b(l),c(l),d(l),cfy(4*m),cfz(4*n),fctrd(ldimfc,n,m),wsavey(m+16),wsavez(n+16),zrt(n)
INTEGER:: IGRID, L, LP, M, MP, N, NP, LDIMFC, K, I
REAL(EB):: C2


!                               INITIALIZATION ROUTINE FOR FFT SOLVERS

PI = FSH19S()
EPS = FSH20S()

IF (L>1 .AND. LP==1) THEN
  LH = (L+1)/2
  LODD = 1
  
  IF (2*LH==L) LODD = 2
  C(LH-1) = 0._EB 
  A(LH) = 0._EB 
  C(LH) = 2._EB*C(LH)
  
  IF (LODD==1) THEN
    B(LH-1) = B(LH-1) - A(LH-1)
    B(L) = B(L) + A(L)
  END IF
  
  IF (LODD==2) A(L) = C(LH)
END IF

!                               COMPUTE TRANSFORM ROOTS FOR J-DIRECTION

IF (M==1) THEN
  CFY(1) = 0._EB 
ELSE
  
  IF (IGRID==1) THEN
    MRDEL = ((MP-1)* (MP-3)* (MP-5))/3
    DEL = PI/ (2._EB* (M+MRDEL))
  ELSE
    DEL = PI/ (2*M)
  END IF
  
  IF (MP==1) THEN
    CFY(1) = 0._EB 
    CFY(M) = -4._EB 
    DO J = 2,M - 1,2
      CFY(J) = -4._EB*SIN(J*DEL)**2
      CFY(J+1) = CFY(J)
    END DO
  END IF
  
  IF (MP==2) THEN
    DO J = 1,M
      CFY(J) = -4._EB*SIN(J*DEL)**2
    END DO
  END IF
  
  IF (MP==3 .OR. MP==5) THEN
    DO J = 1,M
      CFY(J) = -4._EB*SIN((J-.5_EB)*DEL)**2
    END DO
  END IF
  
  IF (MP==4) THEN
    DO J = 1,M
      CFY(J) = -4._EB*SIN((J-1)*DEL)**2
    END DO
  END IF
  
END IF

!                               COMPUTE TRANSFORM ROOTS IN K-DIRECTION

IF (N==1) THEN
  ZRT(1) = 0._EB 
ELSE
  
  IF (IGRID==1) THEN
    NRDEL = ((NP-1)* (NP-3)* (NP-5))/3
    DEL = PI/ (2._EB* (N+NRDEL))
  ELSE
    DEL = PI/ (2*N)
  END IF
  
  IF (NP==1) THEN
    ZRT(1) = 0._EB 
    ZRT(N) = -4._EB*C2
    DO K = 2,N - 1,2
      ZRT(K) = -4._EB*C2*SIN(K*DEL)**2
      ZRT(K+1) = ZRT(K)
    END DO
  END IF
  
  IF (NP==2) THEN
    DO K = 1,N
      ZRT(K) = -4._EB*C2*SIN(K*DEL)**2
    END DO
  END IF
  
  IF (NP==3 .OR. NP==5) THEN
    DO K = 1,N
      ZRT(K) = -4._EB*C2*SIN((K-.5_EB)*DEL)**2
    END DO
  END IF
  
  IF (NP==4) THEN
    DO K = 1,N
      ZRT(K) = -4._EB*C2*SIN((K-1)*DEL)**2
    END DO
  END IF
  
END IF

IF (L>1) THEN
  
!                               FACTOR M*N TRIDIAGONAL SYSTEMS.
!                               FIRST, DO THE POSSIBLY SINGULAR
!                               CASE CORRESPONDING TO J = K = 1.
  
  FCTRD(1,1,1) = 1._EB/ (B(1)+D(1)*CFY(1)+ZRT(1))
  DO  I = 2,L - 1
    FCTRD(I,1,1) = 1._EB/ (B(I)+D(I)*CFY(1)+ZRT(1)- A(I)*C(I-1)*FCTRD(I-1,1,1))
  END DO
  
!                               IF TRIDIAGONAL SYSTEM
!                               (...,A(I),B(I),C(I),...) IS SINGULAR
!                               THEN FCTRD(1,1,L) IS 1._EB/0.  IF
!                               DENOMINATOR IS WITHIN ROUND-OFF OF 0,
!                               SET FCTRD(1,1,L) ARBITRARILY
  
  DEN = B(L) + D(L)*CFY(1) + ZRT(1) - A(L)*C(L-1)*FCTRD(L-1,1,1)
  BMAX = ABS(B(1))
  DO  I = 2,L
    BMAX = MAX(BMAX,ABS(B(I)))
  END DO
  
  IF (ABS(DEN/BMAX)<=10._EB *EPS) DEN = BMAX
  FCTRD(L,1,1) = 1._EB/DEN
  
!                               FACTOR CASES J=1, K=2,...,N.
  
  DO  K = 2,N
    FCTRD(1,K,1) = 1._EB/ (B(1)+D(1)*CFY(1)+ZRT(K))
  END DO
  DO  I = 2,L
    DO  K = 2,N
      FCTRD(I,K,1) = 1._EB/ (B(I)+D(I)*CFY(1)+ZRT(K)- A(I)*C(I-1)*FCTRD(I-1,K,1))
    END DO
  END DO
  
!                               FACTOR CASES K=1, J=2,...,M.
  
  DO  J = 2,M
    FCTRD(1,1,J) = 1._EB/ (B(1)+D(1)*CFY(J)+ZRT(1))
  END DO
  DO  I = 2,L
    DO  J = 2,M
      FCTRD(I,1,J) = 1._EB/ (B(I)+D(I)*CFY(J)+ZRT(1)- A(I)*C(I-1)*FCTRD(I-1,1,J))
    END DO
  END DO
  
!                               FACTOR REMAINING CASES.
  
  DO  K = 2,N
    DO  J = 2,M
      FCTRD(1,K,J) = 1._EB/ (B(1)+D(1)*CFY(J)+ZRT(K))
    END DO
  END DO
  DO  I = 2,L
    DO  K = 2,N
      DO  J = 2,M
        FCTRD(I,K,J) = 1._EB/ (B(I)+D(I)*CFY(J)+ZRT(K)-  &
            A(I)*C(I-1)*FCTRD(I-1,K,J))
      END DO
    END DO
  END DO
END IF

!                               INITIALIZE FFT TRANSFORMS AND
!                               PRE-PROCESSING COEFFICIENTS IN J

IF (M/=1) THEN
  
  SELECT CASE(IGRID)
  CASE(1)   ; GO TO 440
  CASE(2)   ; GO TO 450
END SELECT

440     CONTINUE
SELECT CASE(MP)
CASE(1)   ; GO TO 460
CASE(2)   ; GO TO 470
CASE(3)   ; GO TO 480
CASE(4)   ; GO TO 500
CASE(5)   ; GO TO 510
END SELECT

450     CONTINUE
SELECT CASE(MP)
CASE(1)   ; GO TO 460
CASE(2)   ; GO TO 480
CASE(3)   ; GO TO 490
CASE(4)   ; GO TO 510
CASE(5)   ; GO TO 520
END SELECT

460     CONTINUE
CALL VSRFTI(M,WSAVEY)

GO TO 530

470     CONTINUE
CALL VSINTI(M,CFY,WSAVEY)

GO TO 530

480     CONTINUE
!EB          CALL VSSINI(M,CFY(1),CFY(M+1),WSAVEY)
CALL VSCOSI(M,CFY(1),CFY(M+1),WSAVEY)

GO TO 530

490     CONTINUE
!EB          CALL VSSNQI(M,CFY(1),CFY(M+1),CFY(2*M+1),CFY(3*M+1),WSAVEY)
CALL VSCSQI(M,CFY(1),CFY(M+1),CFY(2*M+1),CFY(3*M+1),WSAVEY)

GO TO 530

500     CONTINUE
CALL VCOSTI(M,CFY,WSAVEY)

GO TO 530

510     CONTINUE
CALL VSCOSI(M,CFY(1),CFY(M+1),WSAVEY)

GO TO 530

520     CONTINUE
CALL VSCSQI(M,CFY(1),CFY(M+1),CFY(2*M+1),CFY(3*M+1),WSAVEY)
530     CONTINUE
END IF

!                               INITIALIZE FFT TRANSFORMS AND
!                               PRE-PROCESSING COEFFICIENTS IN K

IF (N/=1) THEN
  
  SELECT CASE(IGRID)
  CASE(1)   ; GO TO 540
  CASE(2)   ; GO TO 550
END SELECT

540     CONTINUE
SELECT CASE(NP)
CASE(1)   ; GO TO 560
CASE(2)   ; GO TO 570
CASE(3)   ; GO TO 580
CASE(4)   ; GO TO 600
CASE(5)   ; GO TO 610
END SELECT

550     CONTINUE
SELECT CASE(NP)
CASE(1)   ; GO TO 560
CASE(2)   ; GO TO 580
CASE(3)   ; GO TO 590
CASE(4)   ; GO TO 610
CASE(5)   ; GO TO 620
END SELECT

560     CONTINUE
CALL VSRFTI(N,WSAVEZ)

GO TO 630

570     CONTINUE
CALL VSINTI(N,CFZ,WSAVEZ)

GO TO 630

580     CONTINUE
!EB          CALL VSSINI(N,CFZ(1),CFZ(N+1),WSAVEZ)
CALL VSCOSI(N,CFZ(1),CFZ(N+1),WSAVEZ)

GO TO 630

590     CONTINUE
!EB          CALL VSSNQI(N,CFZ(1),CFZ(N+1),CFZ(2*N+1),CFZ(3*N+1),WSAVEZ)
CALL VSCSQI(N,CFZ(1),CFZ(N+1),CFZ(2*N+1),CFZ(3*N+1),WSAVEZ)

GO TO 630

600     CONTINUE
CALL VCOSTI(N,CFZ,WSAVEZ)

GO TO 630

610     CONTINUE
CALL VSCOSI(N,CFZ(1),CFZ(N+1),WSAVEZ)

GO TO 630

620     CONTINUE
CALL VSCSQI(N,CFZ(1),CFZ(N+1),CFZ(2*N+1),CFZ(3*N+1),WSAVEZ)
630     CONTINUE
END IF

RETURN
END SUBROUTINE FSH01S


SUBROUTINE FSH02S(LDIMF,MDIMF,F,SAVE,W)

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



INTEGER:: LDIMF, MDIMF
REAL(EB)   F(LDIMF,MDIMF,*), W(*), SAVE(-3:*)

!                               RETRIEVE CONSTANTS FROM SAVE ARRAY

L = SAVE(2)
LP = SAVE(3)
M = SAVE(4)
MP = SAVE(5)
N = SAVE(6)
NP = SAVE(7)
IGRID = SAVE(11)

!                               ALLOCATION OF SAVE ARRAY
IA = 12
IC = IA + 2*L
ICFY = IC + L

IF (IGRID==1) THEN
  ICFZ = ICFY + 2*M
  IFCTRD = ICFZ + 2*N
ELSE
  ICFZ = ICFY + 4*M
  IFCTRD = ICFZ + 4*N
END IF
LDIMFC=L
IF (LDIMF>L .AND. MOD(L,2)==0) LDIMFC=L+1
IWSY = IFCTRD + LDIMFC*M*N
IWSZ = IWSY + M + 16
LDIMFT=L

IF (LDIMF==L .AND. MDIMF==M) THEN
  
!                               NO HOLES IN DATA ARRAY, SO CALL SOLVER
  
  CALL FSH03S(IGRID,L,LP,M,MP,N,NP,LDIMFC,LDIMFT,F,SAVE(ICFY), SAVE(ICFZ),  &
      W,SAVE(IA),SAVE(IC),SAVE(IFCTRD),SAVE(IWSY), SAVE(IWSZ))
ELSE
  IF (LDIMF>L .AND. MOD(L,2)==0) LDIMFT=L+1
  
!                               PACK DATA ARRAY, CALL SOLVER,
!                               AND THEN UNPACK SOLUTION ARRAY
  
  CALL FSH04S(L,M,N,LDIMF,MDIMF,LDIMFT,F,W)
  CALL FSH03S(IGRID,L,LP,M,MP,N,NP,LDIMFC,LDIMFT,W,SAVE(ICFY),  &
      SAVE(ICFZ),F, SAVE(IA),SAVE(IC),SAVE(IFCTRD),SAVE(IWSY),  &
      SAVE(IWSZ))
  CALL FSH05S(L,M,N,LDIMF,MDIMF,LDIMFT,F,W)
END IF

RETURN
END SUBROUTINE FSH02S


SUBROUTINE FSH03S(IGRID,L,LP,M,MP,N,NP,LDIMFC,LDIMFT,F,CFY,CFZ,  &
    FT,A,C,FCTRD,WSAVEY,WSAVEZ)

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


INTEGER                   :: IGRID
INTEGER                       :: L
INTEGER                   :: LP
INTEGER                       :: M
INTEGER                   :: MP
INTEGER                       :: N
INTEGER                   :: NP
INTEGER                   :: LDIMFC
INTEGER                   :: LDIMFT
REAL(EB)  CFY(4*M)
REAL(EB)  CFZ(4*N)
REAL(EB)  FT(LDIMFT,M,N)
REAL(EB)   FCTRD(LDIMFC,M,N)
REAL(EB)   WSAVEY(M+16)
REAL(EB)   WSAVEZ(N+16)
REAL(EB)   A(L),C(L),F(LDIMFT,M,N)
    
LOGICAL :: DATARY,DATASW
!                               ZERO OUT BOTTOM PLANE OF ARRAY FT

DO  K=1,N
  DO  J=1,M
    FT(LDIMFT,J,K)=0._EB 
  END DO
END DO

NOCOPY=.TRUE.