pres.f90

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

      IF (LBC==5 .OR. LBC==6)             H(0,J,K) = H(1,J,K)
   ENDDO
ENDDO
 
DO K=1,KBAR
   DO I=1,IBAR
      IF (MBC==3 .OR. MBC==4) H(I,0,K)    = H(I,1,K)    - DETA*BYS(I,K)
      IF (MBC==3 .OR. MBC==2) H(I,JBP1,K) = H(I,JBAR,K) + DETA*BYF(I,K)
      IF (MBC==1 .OR. MBC==2) H(I,0,K)    =-H(I,1,K)    + 2._EB*BYS(I,K)
      IF (MBC==1 .OR. MBC==4) H(I,JBP1,K) =-H(I,JBAR,K) + 2._EB*BYF(I,K)
   ENDDO
ENDDO
 
DO J=1,JBAR
   DO I=1,IBAR
      IF (NBC==3 .OR. NBC==4)  H(I,J,0)    = H(I,J,1)    - DZETA*BZS(I,J)
      IF (NBC==3 .OR. NBC==2)  H(I,J,KBP1) = H(I,J,KBAR) + DZETA*BZF(I,J)
      IF (NBC==1 .OR. NBC==2)  H(I,J,0)    =-H(I,J,1)    + 2._EB*BZS(I,J)
      IF (NBC==1 .OR. NBC==4)  H(I,J,KBP1) =-H(I,J,KBAR) + 2._EB*BZF(I,J)
   ENDDO
ENDDO
 
! ************************* Check the Solution *************************
 
IF (CHECK_POISSON) THEN     
   POIS_ERR = 0.
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            RHSS = (R(I-1)*FVX(I-1,J,K)-R(I)*FVX(I,J,K))*RDX(I)*RRN(I) + (FVY(I,J-1,K)-FVY(I,J,K))*RDY(J) + &
                   (FVZ(I,J,K-1)-FVZ(I,J,K))*RDZ(K) - DDDT(I,J,K)
            LHSS = ((H(I+1,J,K)-H(I,J,K))*RDXN(I)*R(I) - (H(I,J,K)-H(I-1,J,K))*RDXN(I-1)*R(I-1) )*RDX(I)*RRN(I) &
                   +((H(I,J+1,K)-H(I,J,K))*RDYN(J) -  (H(I,J,K)-H(I,J-1,K))*RDYN(J-1) )*RDY(J) &
                   +((H(I,J,K+1)-H(I,J,K))*RDZN(K) - (H(I,J,K)-H(I,J,K-1))*RDZN(K-1) )*RDZ(K)
            RES = ABS(RHSS-LHSS)
            POIS_ERR = MAX(RES,POIS_ERR)
         ENDDO
      ENDDO
   ENDDO
ENDIF
!
! **********************************************************************
!
TUSED(5,NM)=TUSED(5,NM)+SECOND()-TNOW
END SUBROUTINE PRESSURE_SOLVER
 
 
SUBROUTINE COMPUTE_A_B(A,B,NM)
USE GLOBAL_CONSTANTS, ONLY: NCGC, SOLID_BOUNDARY, OPEN_BOUNDARY, PREDICTOR
REAL(EB) :: A(NCGC,NCGC),B(NCGC),FVX_AVG,FVY_AVG,FVZ_AVG
TYPE (MESH_TYPE), POINTER :: M2
TYPE (OMESH_TYPE), POINTER :: OM
INTEGER :: I,J,K,NM,II,JJ,KK,IIO,JJO,KKO,NOM,IW,IC,JC,KC,N,NO
 
CALL POINT_TO_MESH(NM)
 
K_COARSE: DO KC=1,KBAR2
   J_COARSE: DO JC=1,JBAR2
      I_COARSE: DO IC=1,IBAR2
 
         N = CGI2(IC,JC,KC)
         II = I_LO(IC)-1
         DO KK=K_LO(KC),K_HI(KC)
            DO JJ=J_LO(JC),J_HI(JC)
               IF (II>0) THEN
                  NO = CGI(II,JJ,KK)
                  A(N,N)  = A(N,N)  - DY(JJ)*DZ(KK)*RDXN(II)
                  A(N,NO) = A(N,NO) + DY(JJ)*DZ(KK)*RDXN(II)
                  B(N)    = B(N)    + DY(JJ)*DZ(KK)*FVX(II,JJ,KK)
               ELSE 
                  IW = WALL_INDEX(CELL_INDEX(II+1,JJ,KK),-1)
                  NOM = IJKW(9,IW)
                  IF (NOM>0) THEN
                     OM  => OMESH(NOM)
                     M2  => MESHES(NOM)
                     IIO = IJKW(10,IW)
                     JJO = IJKW(11,IW)
                     KKO = IJKW(12,IW)
                     NO  = M2%CGI(IIO,JJO,KKO)
                     A(N,N)  = A(N,N)  - DY(JJ)*DZ(KK)*RDXN(II)
                     A(N,NO) = A(N,NO) + DY(JJ)*DZ(KK)*RDXN(II)
                     FVX_AVG = 0.5_EB*(OM%FVX(M2%IBAR,JJO,KKO)+FVX(II,JJ,KK))
                     B(N)    = B(N)    + DY(JJ)*DZ(KK)*FVX_AVG 
                  ENDIF
                  IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) B(N) = B(N) + DUWDT(IW)*DY(JJ)*DZ(KK)
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
                     A(N,N) = A(N,N) - 2._EB*RDXN(II)   *DY(JJ)*DZ(KK)
                     B(N)   = B(N)   + FVX(II,JJ,KK) *DY(JJ)*DZ(KK)- RDXN(II)*DY(JJ)*DZ(KK)*(H(II,JJ,KK)+H(II+1,JJ,KK))
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
 
         II = I_HI(IC)
         DO KK=K_LO(KC),K_HI(KC)
            DO JJ=J_LO(JC),J_HI(JC)
               IF (II<IBAR) THEN
                  NO = CGI(II+1,JJ,KK)
                  A(N,N)  = A(N,N)  - DY(JJ)*DZ(KK)*RDXN(II)
                  A(N,NO) = A(N,NO) + DY(JJ)*DZ(KK)*RDXN(II)
                  B(N)    = B(N)    - DY(JJ)*DZ(KK)*FVX(II,JJ,KK)
               ELSE
                  IW = WALL_INDEX(CELL_INDEX(II,JJ,KK),1)
                  NOM = IJKW(9,IW)
                  IF (NOM>0) THEN
                     OM  => OMESH(NOM)
                     M2  => MESHES(NOM)
                     IIO = IJKW(10,IW)
                     JJO = IJKW(11,IW)
                     KKO = IJKW(12,IW)
                     NO  = M2%CGI(IIO,JJO,KKO)
                     A(N,N)  = A(N,N)  - DY(JJ)*DZ(KK)*RDXN(II)
                     A(N,NO) = A(N,NO) + DY(JJ)*DZ(KK)*RDXN(II)
                     FVX_AVG = 0.5_EB*(OM%FVX(0,JJO,KKO)+FVX(II,JJ,KK))
                     B(N)    = B(N)    - DY(JJ)*DZ(KK)*FVX_AVG
                  ENDIF
                  IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) B(N) = B(N) + DUWDT(IW)*DY(JJ)*DZ(KK)
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
                     A(N,N) = A(N,N) - 2._EB*RDXN(II)  *DY(JJ)*DZ(KK)
                     B(N)   = B(N)   - FVX(II,JJ,KK)*DY(JJ)*DZ(KK) - RDXN(II)*DY(JJ)*DZ(KK)*(H(II,JJ,KK)+H(II+1,JJ,KK))
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
 
         JJ = J_LO(JC)-1
         DO KK=K_LO(KC),K_HI(KC)
            DO II=I_LO(IC),I_HI(IC)
               IF (JJ>0) THEN
                  NO = CGI(II,JJ,KK)
                  A(N,N)  = A(N,N)  - DX(II)*DZ(KK)*RDYN(JJ)
                  A(N,NO) = A(N,NO) + DX(II)*DZ(KK)*RDYN(JJ)
                  B(N)    = B(N)    + DX(II)*DZ(KK)*FVY(II,JJ,KK)
               ELSE
                  IW = WALL_INDEX(CELL_INDEX(II,JJ+1,KK),-2)
                  NOM = IJKW(9,IW)
                  IF (NOM>0) THEN
                     OM  => OMESH(NOM)
                     M2  => MESHES(NOM)
                     IIO = IJKW(10,IW)
                     JJO = IJKW(11,IW)
                     KKO = IJKW(12,IW)
                     NO  = M2%CGI(IIO,JJO,KKO)
                     A(N,N)  = A(N,N)  - DX(II)*DZ(KK)*RDYN(JJ)
                     A(N,NO) = A(N,NO) + DX(II)*DZ(KK)*RDYN(JJ)
                     FVY_AVG = 0.5_EB*(OM%FVY(IIO,M2%JBAR,KKO)+FVY(II,JJ,KK))
                     B(N)    = B(N)    + DX(II)*DZ(KK)*FVY_AVG
                  ENDIF
                  IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) B(N) = B(N) + DUWDT(IW)*DX(II)*DZ(KK)
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
                     A(N,N) = A(N,N) - 2._EB*RDYN(JJ)  *DX(II)*DZ(KK)
                     B(N)   = B(N)   + FVY(II,JJ,KK)*DX(II)*DZ(KK) - RDYN(JJ)*DX(II)*DZ(KK)*(H(II,JJ,KK)+H(II,JJ+1,KK))
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
 
         JJ = J_HI(JC)
         DO KK=K_LO(KC),K_HI(KC)
            DO II=I_LO(IC),I_HI(IC)
               IF (JJ<JBAR) THEN
                  NO = CGI(II,JJ+1,KK)
                  A(N,N)  = A(N,N)  - DX(II)*DZ(KK)*RDYN(JJ)
                  A(N,NO) = A(N,NO) + DX(II)*DZ(KK)*RDYN(JJ)
                  B(N)    = B(N)    - DX(II)*DZ(KK)*FVY(II,JJ,KK)
               ELSE
                  IW = WALL_INDEX(CELL_INDEX(II,JJ,KK),2)
                  NOM = IJKW(9,IW)
                  IF (NOM>0) THEN
                     OM  => OMESH(NOM)
                     M2  => MESHES(NOM)
                     IIO = IJKW(10,IW)
                     JJO = IJKW(11,IW)
                     KKO = IJKW(12,IW)
                     NO  = M2%CGI(IIO,JJO,KKO)
                     A(N,N)  = A(N,N)  - DX(II)*DZ(KK)*RDYN(JJ)
                     A(N,NO) = A(N,NO) + DX(II)*DZ(KK)*RDYN(JJ)
                     FVY_AVG  = 0.5_EB*(OM%FVY(IIO,0,KKO)+FVY(II,JJ,KK))
                     B(N)     = B(N)   - DX(II)*DZ(KK)*FVY_AVG
                  ENDIF
                  IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) B(N) = B(N) + DUWDT(IW)*DX(II)*DZ(KK)
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
                     A(N,N) = A(N,N) - 2._EB*RDYN(JJ)  *DX(II)*DZ(KK)
                     B(N)   = B(N)   - FVY(II,JJ,KK)*DX(II)*DZ(KK) - RDYN(JJ)*DX(II)*DZ(KK)*(H(II,JJ,KK)+H(II,JJ+1,KK))
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
 
         KK = K_LO(KC)-1
         DO JJ=J_LO(JC),J_HI(JC)
            DO II=I_LO(IC),I_HI(IC)
               IF (KK>0) THEN
                  NO = CGI(II,JJ,KK)
                  A(N,N)  = A(N,N)  - DX(II)*DY(JJ)*RDZN(KK)
                  A(N,NO) = A(N,NO) + DX(II)*DY(JJ)*RDZN(KK)
                  B(N)    = B(N)    + DX(II)*DY(JJ)*FVZ(II,JJ,KK)
               ELSE
                  IW = WALL_INDEX(CELL_INDEX(II,JJ,KK+1),-3)
                  NOM = IJKW(9,IW)
                  IF (NOM>0) THEN
                     OM  => OMESH(NOM)
                     M2  => MESHES(NOM)
                     IIO = IJKW(10,IW)
                     JJO = IJKW(11,IW)
                     KKO = IJKW(12,IW)
                     NO  = M2%CGI(IIO,JJO,KKO)
                     A(N,N)  = A(N,N)  - DX(II)*DY(JJ)*RDZN(KK)
                     A(N,NO) = A(N,NO) + DX(II)*DY(JJ)*RDZN(KK)
                     FVZ_AVG = 0.5_EB*(OM%FVZ(IIO,JJO,M2%KBAR)+FVZ(II,JJ,KK))
                     B(N)    = B(N)    + DX(II)*DY(JJ)*FVZ_AVG
                  ENDIF
                  IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) B(N) = B(N) + DUWDT(IW)*DX(II)*DY(JJ)
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
                     A(N,N) = A(N,N) - 2._EB*RDZN(KK)  *DX(II)*DY(JJ)
                     B(N)   = B(N)   + FVZ(II,JJ,KK)*DX(II)*DY(JJ) - RDZN(KK)*DX(II)*DY(JJ)*(H(II,JJ,KK)+H(II,JJ,KK+1))
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
 
         KK = K_HI(KC)
         DO JJ=J_LO(JC),J_HI(JC)
            DO II=I_LO(IC),I_HI(IC)
               IF (KK<KBAR) THEN
                  NO = CGI(II,JJ,KK+1)
                  A(N,N)  = A(N,N)  - DX(II)*DY(JJ)*RDZN(KK)
                  A(N,NO) = A(N,NO) + DX(II)*DY(JJ)*RDZN(KK)
                  B(N)    = B(N)    - DX(II)*DY(JJ)*FVZ(II,JJ,KK)
               ELSE
                  IW = WALL_INDEX(CELL_INDEX(II,JJ,KK),3)
                  NOM = IJKW(9,IW)
                  IF (NOM>0) THEN
                     OM  => OMESH(NOM)
                     M2  => MESHES(NOM)
                     IIO = IJKW(10,IW)
                     JJO = IJKW(11,IW)
                     KKO = IJKW(12,IW)
                     NO  = M2%CGI(IIO,JJO,KKO)
                     A(N,N)  = A(N,N)  - DX(II)*DY(JJ)*RDZN(KK)
                     A(N,NO) = A(N,NO) + DX(II)*DY(JJ)*RDZN(KK)
                     FVZ_AVG = 0.5_EB*(OM%FVZ(IIO,JJO,0)+FVZ(II,JJ,KK))
                     B(N)    = B(N)    - DX(II)*DY(JJ)*FVZ_AVG
                  ENDIF
                  IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) B(N) = B(N) + DUWDT(IW)*DX(II)*DY(JJ)
                  IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
                     A(N,N) = A(N,N) - 2._EB*RDZN(KK)  *DX(II)*DY(JJ)
                     B(N)   = B(N)   - FVZ(II,JJ,KK)*DX(II)*DY(JJ) - RDZN(KK)*DX(II)*DY(JJ)*(H(II,JJ,KK)+H(II,JJ,KK+1))
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
 
         DO K=K_LO(KC),K_HI(KC)
            DO J=J_LO(JC),J_HI(JC)
               DO I=I_LO(IC),I_HI(IC)
                  B(N) = B(N) - DDDT(I,J,K)*DX(I)*RC(I)*DY(J)*DZ(K)
               ENDDO
            ENDDO
         ENDDO
 
      ENDDO I_COARSE
   ENDDO J_COARSE
ENDDO K_COARSE
 
END SUBROUTINE COMPUTE_A_B
 

SUBROUTINE COMPUTE_CORRECTION_PRESSURE(B,NM)
USE GLOBAL_CONSTANTS, ONLY: NCGC, SOLID_BOUNDARY, OPEN_BOUNDARY, PREDICTOR, TWO_D, CYLINDRICAL
USE POIS, ONLY: H3CZSS, H2CZSS
 
REAL(EB) :: B(NCGC),DHDX_F,DHDX_S,DHDY_F, &
            DHDY_S,DHDZ_F,DHDZ_S,DA, &
            AREA_XS,AREA_XF,AREA_YS,AREA_YF,AREA_ZS,AREA_ZF, &
            DUDT,DUDTO,DVDT,DVDTO,DWDT,DWDTO,RDT, &
            U_NEXT,U_NEXT_O,V_NEXT,V_NEXT_O,W_NEXT,W_NEXT_O, &
            B_OTHER,FVX_AVG,FVY_AVG,FVZ_AVG,D_FAC, &
            AREA_XS_CLOSED,AREA_XF_CLOSED, &
            AREA_YS_CLOSED,AREA_YF_CLOSED, &
            AREA_ZS_CLOSED,AREA_ZF_CLOSED
INTEGER :: NM,II,JJ,KK,IW,IOR,I,J,K,IIO,JJO,KKO,NOM, &
            IOR_PATCH,II_LOW,II_HIGH,JJ_LOW,JJ_HIGH, KK_LOW,KK_HIGH,IC,JC,KC,N,NO
TYPE (MESH_TYPE), POINTER :: M2
TYPE (OMESH_TYPE), POINTER :: OM
 
CALL POINT_TO_MESH(NM)
 
RDT = 1._EB/DT
!     D_FAC = 0.5
D_FAC = 0.0_EB
 
BXS = 0._EB
BXF = 0._EB
BYS = 0._EB
BYF = 0._EB
BZS = 0._EB
BZF = 0._EB
 
ORIENT_LOOP:  DO IOR_PATCH=-3,3
   IF (IOR_PATCH==0) CYCLE ORIENT_LOOP
   KC_LOOP: DO KC=1,KBAR2
      IF (IOR_PATCH== 3 .AND. KC/=1)     CYCLE KC_LOOP
      IF (IOR_PATCH==-3 .AND. KC/=KBAR2) CYCLE KC_LOOP
      JC_LOOP: DO JC=1,JBAR2
         IF (IOR_PATCH== 2 .AND. JC/=1)     CYCLE JC_LOOP
         IF (IOR_PATCH==-2 .AND. JC/=JBAR2) CYCLE JC_LOOP
         IC_LOOP: DO IC=1,IBAR2
            IF (IOR_PATCH== 1 .AND. IC/=1)     CYCLE IC_LOOP
            IF (IOR_PATCH==-1 .AND. IC/=IBAR2) CYCLE IC_LOOP
!
            N = CGI2(IC,JC,KC)
!
            SELECT CASE(IOR_PATCH)
!
               CASE(1)
                  II_LOW    = I_LO(IC)-1
                  II_HIGH   = I_LO(IC)-1
                  JJ_LOW    = J_LO(JC)
                  JJ_HIGH   = J_HI(JC)
                  KK_LOW    = K_LO(KC)
                  KK_HIGH   = K_HI(KC)
               CASE(-1)
                  II_LOW    = I_HI(IC)+1
                  II_HIGH   = I_HI(IC)+1
                  JJ_LOW    = J_LO(JC)
                  JJ_HIGH   = J_HI(JC)
                  KK_LOW    = K_LO(KC)
                  KK_HIGH   = K_HI(KC)
               CASE(2)
                  II_LOW    = I_LO(IC)
                  II_HIGH   = I_HI(IC)
                  JJ_LOW    = J_LO(JC)-1
                  JJ_HIGH   = J_LO(JC)-1
                  KK_LOW    = K_LO(KC)
                  KK_HIGH   = K_HI(KC)
               CASE(-2)
                  II_LOW    = I_LO(IC)
                  II_HIGH   = I_HI(IC)
                  JJ_LOW    = J_HI(JC)+1
                  JJ_HIGH   = J_HI(JC)+1
                  KK_LOW    = K_LO(KC)
                  KK_HIGH   = K_HI(KC)
               CASE(3)
                  II_LOW    = I_LO(IC)
                  II_HIGH   = I_HI(IC)
                  JJ_LOW    = J_LO(JC)
                  JJ_HIGH   = J_HI(JC)
                  KK_LOW    = K_LO(KC)-1
                  KK_HIGH   = K_LO(KC)-1
               CASE(-3)
                  II_LOW    = I_LO(IC)
                  II_HIGH   = I_HI(IC)
                  JJ_LOW    = J_LO(JC)
                  JJ_HIGH   = J_HI(JC)
                  KK_LOW    = K_HI(KC)+1
                  KK_HIGH   = K_HI(KC)+1
!
            END SELECT
!
            DHDX_S = 0._EB
            DHDX_F = 0._EB
            DHDY_S = 0._EB 
            DHDY_F = 0._EB 
            DHDZ_S = 0._EB 
            DHDZ_F = 0._EB 
!
            AREA_XS_CLOSED=0._EB
            AREA_XF_CLOSED=0._EB
            AREA_YS_CLOSED=0._EB
            AREA_YF_CLOSED=0._EB
            AREA_ZS_CLOSED=0._EB
            AREA_ZF_CLOSED=0._EB
!
            AREA_XS = 0._EB 
            AREA_XF = 0._EB 
            AREA_YS = 0._EB 
            AREA_YF = 0._EB 
            AREA_ZS = 0._EB 
            AREA_ZF = 0._EB 
!
            WALL_CELL_LOOP: DO IW=1,NEWC
!
               II  = IJKW(1,IW)
               JJ  = IJKW(2,IW)
               KK  = IJKW(3,IW)
               IOR = IJKW(4,IW)
               IF (IOR/=IOR_PATCH) CYCLE WALL_CELL_LOOP
               IF (II<II_LOW .OR. II>II_HIGH) CYCLE WALL_CELL_LOOP
               IF (JJ<JJ_LOW .OR. JJ>JJ_HIGH) CYCLE WALL_CELL_LOOP
               IF (KK<KK_LOW .OR. KK>KK_HIGH) CYCLE WALL_CELL_LOOP
               NOM = IJKW(9,IW)
               IIO = IJKW(10,IW)
               JJO = IJKW(11,IW)
               KKO = IJKW(12,IW)
               OM  => OMESH(NOM)
               M2  => MESHES(NOM)

               IF_INTERPOLATED_BOUNDARY: IF (NOM>0) THEN
                  NO = M2%CGI(IIO,JJO,KKO)