FDS Source: update Z_LS calculation to properly account for terrain c…

Source/vege.f90

@@ -18,6 +18,7 @@ MODULE VEGE
 TYPE(MESH_TYPE),    POINTER :: M
 TYPE(WALL_TYPE),    POINTER :: WC
 TYPE(CFACE_TYPE),   POINTER :: CFA
+TYPE(CC_CUTFACE_TYPE),   POINTER :: CF
 TYPE(SURFACE_TYPE), POINTER :: SF
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
 TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2
@@ -37,7 +38,7 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_1(NM)
 USE COMPLEX_GEOMETRY, ONLY : CC_IDCF
 INTEGER, INTENT(IN) :: NM
 INTEGER :: ICF,IW,I,J,SURF_INDEX
-REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: SUM_AREA
+REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: MAX_AREA,SUM_AREA
 
 T_NOW = CURRENT_TIME()
 
@@ -78,38 +79,43 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_1(NM)
 
 IF (CC_IBM) THEN
 
+   ALLOCATE(MAX_AREA(0:IBP1,0:JBP1)) ; MAX_AREA = 0._EB
    ALLOCATE(SUM_AREA(0:IBP1,0:JBP1)) ; SUM_AREA = 0._EB
    ALLOCATE(M%LS_KLO_TERRAIN(0:IBP1,0:JBP1),STAT=IZERO) ; CALL ChkMemErr('READ','LS_KLO_TERRAIN',IZERO)
    LS_KLO_TERRAIN => M%LS_KLO_TERRAIN ; LS_KLO_TERRAIN = 2*KBP1+1 ! Number larger that KBP1.
    ALLOCATE(M%LS_KHI_TERRAIN(0:IBP1,0:JBP1),STAT=IZERO) ; CALL ChkMemErr('READ','LS_KHI_TERRAIN',IZERO)
    LS_KHI_TERRAIN => M%LS_KHI_TERRAIN ; LS_KHI_TERRAIN = -1 ! Number smaller than 0.
    DO ICF=1,M%N_CUTFACE_MESH
-      IF (CUT_FACE(ICF)%STATUS/=2 .OR. CUT_FACE(ICF)%NFACE<1) CYCLE ! CC_INBOUNDARY == 2
+      CF  => CUT_FACE(ICF)
+      IF (CF%STATUS/=2 .OR. CF%NFACE<1) CYCLE ! CC_INBOUNDARY == 2
       ! Location of CFACE with largest AREA, to define SURF_INDEX:
-      IW  = MAXLOC(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE),DIM=1)
-      CFA => CFACE( CUT_FACE(ICF)%CFACE_INDEX(IW) )
-      BC  => BOUNDARY_COORD(CFA%BC_INDEX)
-      IF (BC%KKG < LS_KLO_TERRAIN(BC%IIG,BC%JJG)) LS_KLO_TERRAIN(BC%IIG,BC%JJG) = BC%KKG
-      IF (BC%KKG > LS_KHI_TERRAIN(BC%IIG,BC%JJG)) LS_KHI_TERRAIN(BC%IIG,BC%JJG) = BC%KKG
+      IW  = MAXLOC(CF%AREA(1:CF%NFACE),DIM=1)
+      CFA => CFACE(CF%CFACE_INDEX(IW) )
+      BC  => BOUNDARY_COORD(CFA%BC_INDEX)      
       IF (BC%NVEC(KAXIS)>-TWO_EPSILON_EB .AND. CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
-         IF (SUM(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE))<SUM_AREA(BC%IIG,BC%JJG)) THEN
+         IF (BC%KKG < LS_KLO_TERRAIN(BC%IIG,BC%JJG)) LS_KLO_TERRAIN(BC%IIG,BC%JJG) = BC%KKG
+         IF (BC%KKG > LS_KHI_TERRAIN(BC%IIG,BC%JJG)) LS_KHI_TERRAIN(BC%IIG,BC%JJG) = BC%KKG
+         SUM_AREA(BC%IIG,BC%JJG) = SUM_AREA(BC%IIG,BC%JJG) + SUM(CF%AREA(1:CF%NFACE))
+         Z_LS(BC%IIG,BC%JJG) = Z_LS(BC%IIG,BC%JJG) + DOT_PRODUCT(CF%XYZCEN(KAXIS,1:CF%NFACE),CF%AREA(1:CF%NFACE))
+         IF (SUM(CF%AREA(1:CF%NFACE))<MAX_AREA(BC%IIG,BC%JJG)) THEN
             CYCLE  ! This CUT_FACE does not contain the majority of the area corresponding to the FDS cell area, DX*DY
          ELSE
-            SUM_AREA(BC%IIG,BC%JJG) = SUM(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE))
+            MAX_AREA(BC%IIG,BC%JJG) = SUM(CF%AREA(1:CF%NFACE))
          ENDIF
-         ! Area averaged Z height of CFACES within this cut-cell (containing CC_INBOUNDARY CFACES):
-         Z_LS(BC%IIG,BC%JJG) = DOT_PRODUCT(CUT_FACE(ICF)%XYZCEN(KAXIS,1:CUT_FACE(ICF)%NFACE), &
-                                           CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE)       ) / SUM_AREA(BC%IIG,BC%JJG)
          IF (BC%KKG > K_LS(BC%IIG,BC%JJG)) K_LS(BC%IIG,BC%JJG) = BC%KKG
          LS_SURF_INDEX(BC%IIG,BC%JJG) = CFA%SURF_INDEX
       ENDIF
    ENDDO
+   WHERE (SUM_AREA>TWO_EPSILON_EB)  Z_LS = Z_LS/SUM_AREA
+   DEALLOCATE(MAX_AREA)
+   DEALLOCATE(SUM_AREA)
+
    DO J=1,JBAR
       DO I=1,IBAR
          IF (K_LS(I,J)==KBAR .AND. FCVAR(I,J,K_LS(I,J),CC_IDCF,KAXIS)>0) LS_SURF_INDEX(I,J) = 0
       ENDDO
    ENDDO
-   DEALLOCATE(SUM_AREA)
+
 ENDIF
 
 ! Set up level set on cartesian faces only where they are not under a GEOM
@@ -155,7 +161,7 @@ END SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_1
 SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_2(NM,MODE)
 
 INTEGER, INTENT(IN) :: NM,MODE
-INTEGER :: I,IM1,IM2,IIG,IP1,IP2,J,JJG,JM1,JP1
+INTEGER :: I,IM1,IIG,IP1,J,JJG,JM1,JP1
 REAL(EB) :: DZT_DUM,DZTDX_DUM,DZTDY_DUM,G_EAST,G_WEST,G_SOUTH,G_NORTH
 
 T_NOW = CURRENT_TIME()
@@ -236,17 +242,13 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_2(NM,MODE)
 
 GRADIENT_ILOOP: DO I = 1,IBAR
 
-   IM1=I-1 ; IM2=I-2
-   IP1=I+1 ; IP2=I+2
-   IF (I==1) IM1 = I
-   IF (I==IBAR) IP1 = I
+   IM1=I-1
+   IP1=I+1
 
    DO J = 1,JBAR
 
       JM1=J-1
       JP1=J+1
-      IF (J==1) JM1 = J
-      IF (J==IBAR) JP1 = J
 
       G_EAST  = 0.5_EB*( Z_LS(I,J) + Z_LS(IP1,J) )
       G_WEST  = 0.5_EB*( Z_LS(I,J) + Z_LS(IM1,J) )
@@ -489,7 +491,6 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
             ! and then calculating effective wind speed, phi_s is converted to an effected wind speed and added
             ! to UMF calculated from the wind. Effective U has units of m/min in Wilson formula.
             ! 0.3048 ~= 1/3.281
-            ! if phi_s < 0 then a complex value (NaN) results. Using abs(phi_s) and sign function to correct.
 
             UMF_TMP = &
             0.3048_EB/PHI_S*(((SF%VEG_LSET_BETA / BETA_OP_ROTH)**E_ROTH)*PHI_S/C_ROTH)**(1._EB/B_ROTH)