pois.f90

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

MODULE POIS
 
! CODE CONVERTED USING TO_F90 BY ALAN MILLER
! DATE: 2006-01-18  TIME: 13:48:46

! POISSON SOLVER ROUTINES

USE PRECISION_PARAMETERS
IMPLICIT REAL(EB) (A-H,O-Z)
IMPLICIT INTEGER (I-N)
PRIVATE
REAL(EB) SCALE
INTEGER :: KAPPA,NMAX,IKPWR
LOGICAL :: TPOSE,NOCOPY,OUTARY
CHARACTER(255), PARAMETER :: poisid='$Id: pois.f90 607 2007-09-17 17:53:03Z mcgratta $'
CHARACTER(255), PARAMETER :: poisrev='$Revision: 607 $'
CHARACTER(255), PARAMETER :: poisdate='$Date: 2007-09-17 13:53:03 -0400 (Mon, 17 Sep 2007) $'

PUBLIC H3CZIS,H3CZSS,H2CZSS,H2CYSS,H3CSSS,H2CZIS,H3CSIS,H2CYIS,GET_REV_pois

CONTAINS

SUBROUTINE H3CZIS(XS,XF,L,LBDCND,YS,YF,M,MBDCND,ZS,ZF,N,NBDCND,  &
    H,ELMBDA,LDIMF,MDIMF,IERROR,SAVE)

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

REAL(EB)::  XS, XF, YS, YF, ZS, ZF, ELMBDA, DLZSQR
INTEGER:: L, LBDCND, M, MBDCND, N, NBDCND, LDIMF, MDIMF, IERROR, LPEROD
REAL(EB) SAVE(-3:*),H(0:*)

!                               CHECK FOR INVALID INPUT

IERROR = 0

IF (XS>XF) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 1._EB
END IF

IF (L<3) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 2._EB
END IF

IF (LBDCND<0 .OR. LBDCND>4) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 3._EB 
END IF

IF (YS>YF) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 4._EB 
END IF

IF (M<3) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 5._EB 
END IF

IF (MBDCND<0 .OR. MBDCND>4) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 6._EB 
END IF

IF (ZS>ZF) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 7._EB 
END IF

IF (N<3) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 8._EB 
END IF

IF (NBDCND<0 .OR. NBDCND>4) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 9._EB 
END IF

IF (LDIMF<L) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 10._EB 
END IF

IF (MDIMF<M) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 11._EB 
END IF

IF (IERROR/=0) THEN
  RETURN
ELSE
  SAVE(1) = 0._EB 
END IF

!                               DEFINE GRID PARAMETERS

DX = (XF-XS)/L
DY = (YF-YS)/M
DZ = (ZF-ZS)/N
DLXSQR = 1._EB/DX**2
DLYSQR = 1._EB/DY**2
DLZSQR = 1._EB/DZ**2
LP = LBDCND + 1
MP = MBDCND + 1
NP = NBDCND + 1

!                               ALLOCATE SAVE ARRAY

IA = 12
IB = IA + L
IC = IB + L
ID = IC + L
IS = ID + L

!                               DEFINE THE A,B,C COEFFICIENTS
!                               IN SAVE ARRAY

DO  I = 0,L - 1
  HM = .5_EB*(H(I)+H(I+1))
  HP = .5_EB*(H(I+1)+H(I+2))
  SAVE(IA+I) = DLXSQR/(H(I+1)*HM)
  SAVE(IC+I) = DLXSQR/(H(I+1)*HP)
  SAVE(IB+I) = - (SAVE(IA+I)+SAVE(IC+I)) + ELMBDA
  SAVE(ID+I) = DLYSQR
END DO

SELECT CASE(LP)
CASE(2)
SAVE(IB) = SAVE(IB) - SAVE(IA)
SAVE(IC-1) = SAVE(IC-1) - SAVE(ID-1)
CASE(3)
SAVE(IB) = SAVE(IB) - SAVE(IA)
SAVE(IC-1) = SAVE(IC-1) + SAVE(ID-1)
CASE(4)
SAVE(IB) = SAVE(IB) + SAVE(IA)
SAVE(IC-1) = SAVE(IC-1) + SAVE(ID-1)
CASE(5)
SAVE(IB) = SAVE(IB) + SAVE(IA)
SAVE(IC-1) = SAVE(IC-1) - SAVE(ID-1)
END SELECT


!                               DETERMINE WHETHER OR NOT BOUNDARY
!                               CONDITION IS PERIODIC IN X

IF (LBDCND==0) THEN
  LPEROD = 0
ELSE
  LPEROD = 1
END IF

!                               INITIALIZE SOLVER ROUTINE S3CFIS.
CALL S3CFIS(LPEROD,L,MBDCND,M,NBDCND,N,DLZSQR,SAVE(IA),SAVE(IB),  &
    SAVE(IC),SAVE(ID),LDIMF,MDIMF,IR,SAVE(IS))

!                               TEST ERROR FLAG FROM S3CFIS FOR
!                               INTERNAL ERROR

IF (IR/=0) THEN
  SAVE(1) = 99._EB 
  IERROR = 1
  
  RETURN
END IF

!                               SAVE PARAMETERS FOR H3CZSS IN SAVE ARRAY

SAVE(2) = DX
SAVE(3) = L
SAVE(4) = LP
SAVE(5) = DY
SAVE(6) = M
SAVE(7) = MP
SAVE(8) = DZ
SAVE(9) = N
SAVE(10) = NP
SAVE(11) = ELMBDA

RETURN
END SUBROUTINE H3CZIS

SUBROUTINE H3CZSS(BDXS,BDXF,BDYS,BDYF,BDZS,BDZF,LDIMF,MDIMF,F,  &
    PERTRB,SAVE,W,H)

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


REAL(EB):: BDXS(MDIMF,*), BDXF(MDIMF,*), BDYS(LDIMF,*), BDYF(LDIMF,*), BDZS(LDIMF,*), BDZF(LDIMF,*), &
           F(LDIMF,MDIMF,*),SAVE(-3:*),W(*),H(0:*), PERTRB
INTEGER:: LDIMF, MDIMF
 

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

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

!                               GET PARAMETERS FOR H3CZSS FROM SAVE
!                               ARRAY WHERE THEY WERE STORED IN
!                               INITIALIZATION SUBROUTINE H3CZIS.

DX = SAVE(2)
L = SAVE(3)
LP = SAVE(4)
DY = SAVE(5)
M = SAVE(6)
MP = SAVE(7)
DZ = SAVE(8)
N = SAVE(9)
NP = SAVE(10)
ELMBDA = SAVE(11)

DLYRCP = 1._EB/DY
TWDYSQ = 2._EB/DY**2
DLZRCP = 1._EB/DZ
TWDZSQ = 2._EB/DZ**2

!                               ALLOCATE SAVE ARRAY

IA = 12
IB = IA + L
IC = IB + L
ID = IC + L
IS = 12 + 4*L


!                               ENTER BOUNDARY DATA FOR X-BOUNDARIES

IF (LP==2 .OR. LP==3) THEN
  DO K = 1,N
    DO J = 1,M
      F(1,J,K) = F(1,J,K) - 2._EB*BDXS(J,K)*SAVE(IA)
    END DO
  END DO
END IF

IF (LP==4 .OR. LP==5) THEN
  DO K = 1,N
    DO J = 1,M
      F(1,J,K) = F(1,J,K) + SAVE(IA)*DX*BDXS(J,K)
    END DO
  END DO
END IF

IF (LP==2 .OR. LP==5) THEN
  DO K = 1,N
    DO J = 1,M
      F(L,J,K) = F(L,J,K) - 2._EB*BDXF(J,K)*SAVE(ID-1)
    END DO
  END DO
END IF

IF (LP==3 .OR. LP==4) THEN
  DO K = 1,N
    DO J = 1,M
      F(L,J,K) = F(L,J,K) - SAVE(ID-1)*DX*BDXF(J,K)
    END DO
  END DO
END IF

!                               ENTER BOUNDARY DATA FOR Y-BOUNDARIES

IF (MP==2 .OR. MP==3) THEN
  DO K = 1,N
    DO I = 1,L
      F(I,1,K) = F(I,1,K) - BDYS(I,K)*TWDYSQ
    END DO
  END DO
END IF

IF (MP==4 .OR. MP==5) THEN
  DO K = 1,N
    DO I = 1,L
      F(I,1,K) = F(I,1,K) + BDYS(I,K)*DLYRCP
    END DO
  END DO
END IF

IF (MP==2 .OR. MP==5) THEN
  DO K = 1,N
    DO I = 1,L
      F(I,M,K) = F(I,M,K) - BDYF(I,K)*TWDYSQ
    END DO
  END DO
END IF

IF (MP==3 .OR. MP==4) THEN
  DO K = 1,N
    DO I = 1,L
      F(I,M,K) = F(I,M,K) - BDYF(I,K)*DLYRCP
    END DO
  END DO
END IF

!                               ENTER BOUNDARY DATA FOR Z-BOUNDARIES

IF (NP==2 .OR. NP==3) THEN
  DO J = 1,M
    DO I = 1,L
      F(I,J,1) = F(I,J,1) - BDZS(I,J)*TWDZSQ
    END DO
  END DO
END IF

IF (NP==4 .OR. NP==5) THEN
  DO J = 1,M
    DO I = 1,L
      F(I,J,1) = F(I,J,1) + BDZS(I,J)*DLZRCP
    END DO
  END DO
END IF

IF (NP==2 .OR. NP==5) THEN
  DO J = 1,M
    DO I = 1,L
      F(I,J,N) = F(I,J,N) - BDZF(I,J)*TWDZSQ
    END DO
  END DO
END IF

IF (NP==3 .OR. NP==4) THEN
  DO J = 1,M
    DO I = 1,L
      F(I,J,N) = F(I,J,N) - BDZF(I,J)*DLZRCP
    END DO
  END DO
END IF

PERTRB = 0._EB 
PERT   = 0._EB 
ISING = 0

!                               FOR SINGULAR PROBLEMS ADJUST DATA TO
!                               INSURE A SOLUTION WILL EXIST.  GO THRU
!                               THIS CODE TWICE: ISING=1 FOR CALCULATING
!                               PERTRB; ISING=2 FOR NORMALIZING SOLUTION
!                               AFTER IT IS COMPUTED.

SELECT CASE(LP)
CASE(1)   ; GO TO 630
CASE(2:3) ; GO TO 750
CASE(4)   ; GO TO 630
CASE(5)   ; GO TO 750
END SELECT

630 CONTINUE
SELECT CASE(MP)
CASE(1)   ; GO TO 640
CASE(2:3) ; GO TO 750
CASE(4)   ; GO TO 640
CASE(5)   ; GO TO 750
END SELECT

640 CONTINUE
SELECT CASE(NP)
CASE(1)   ; GO TO 650
CASE(2:3) ; GO TO 750
CASE(4)   ; GO TO 650
CASE(5)   ; GO TO 750
END SELECT

650 CONTINUE
IF (ELMBDA/=0._EB) GO TO 750
ISING = 1
660 CONTINUE
PERT = 0._EB 
DO  I = 1,L
  W(I) = 0._EB 
END DO
DO  K = 1,N
  DO  J = 1,M
    DO  I = 1,L
      W(I) = W(I) + F(I,J,K)
    END DO
  END DO
END DO
S1 = 0._EB 
S3 = 0._EB 
DO  I = 1,L
  S3 = S3 + H(I)
  S1 = S1 + H(I)*W(I)
END DO

S3   = S3*M*N
PERT = S1/S3


!                               ADJUST F ARRAY BY PERT

DO  K = 1,N
  DO  J = 1,M
    DO  I = 1,L
      F(I,J,K) = F(I,J,K) - PERT
    END DO
  END DO
END DO
750 CONTINUE

!                               IF NORMALIZING SOLUTION, RESTORE PERTRB
!                               AND JUMP TO END

IF (ISING==2) THEN
  PERTRB = PRTSAV
  
  GO TO 800
  
END IF

PRTSAV = PERT

!                               SOLVE THE EQUATION

CALL S3CFSS(LDIMF,MDIMF,F,SAVE(IS),W)

!                               IF A SINGULAR PROBLEM,
!                               RE-NORMALIZE SOLUTION (ISING=2)

IF (ISING==1) THEN
  ISING = 2
  
  GO TO 660
  
END IF

800 CONTINUE
RETURN
END SUBROUTINE H3CZSS


SUBROUTINE S3CFIS(LPEROD,L,MPEROD,M,NPEROD,N,SCAL,A,B,C,D,LDIMF,  &
    MDIMF,IERROR,SAVE)

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



INTEGER:: LPEROD, L, MPEROD, M, NPEROD, N
REAL(EB):: SCAL
REAL(EB):: A(L), B(L), C(L), D(L), SAVE(-3:*)
INTEGER:: LDIMF, MDIMF, IERROR
!                               CHECK FOR INVALID INPUT

IERROR = 0

IF (LPEROD/=0 .AND. LPEROD/=1) THEN
  IERROR = 1
  SAVE(1) = 1._EB 
END IF

IF (L<3) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 2._EB 
END IF

IF (MPEROD<0 .AND. MPEROD>4) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 3._EB 
END IF

IF (M<3) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 4._EB 
END IF

IF (NPEROD<0 .AND. NPEROD>4) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 5._EB 
END IF

IF (N<3) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 6._EB 
END IF

IF (LDIMF<L) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 7._EB 
END IF

IF (MDIMF<M) THEN
  IERROR = IERROR + 1
  SAVE(IERROR) = 8._EB 
END IF

IF (LPEROD==0) THEN
  DO  I = 1,L
    
    IF (A(I)/=A(1)) GO TO 110
    IF (B(I)/=B(1)) GO TO 110
    IF (C(I)/=A(1)) GO TO 110
    IF (D(I)/=D(1)) GO TO 110
    
  END DO
  
  GO TO 120
  
  110     CONTINUE
  IERROR = IERROR + 1
  SAVE(IERROR) = 9._EB 
END IF

120 CONTINUE

IF (IERROR/=0) THEN
  RETURN
ELSE
  SAVE(1) = IERROR
END IF