Merge pull request #14722 from ericvmueller/master

Correct the vector implementation of Rothermel slope factor for level set spread

Author: ericvmueller

Source/mesh.f90

@@ -333,7 +333,7 @@ MODULE MESH_VARIABLES
    REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: FLUX0_LS,FLUX1_LS,PHI_LS,PHI1_LS,ROS_BACKU, &
                                             ROS_HEAD,ROS_FLANK,WIND_EXP, &
                                             SR_X_LS,SR_Y_LS,U_LS,V_LS,Z_LS,DZTDX,DZTDY,MAG_ZT, &
-                                            PHI_WS,UMF,THETA_ELPS,PHI_S,PHI_S_X,PHI_S_Y,PHI_W,LS_WORK1,LS_WORK2
+                                            PHI_WS,UMF,THETA_ELPS,PHI_S_X,PHI_S_Y,PHI_W,LS_WORK1,LS_WORK2
 
 END TYPE MESH_TYPE
 
@@ -475,7 +475,7 @@ MODULE MESH_POINTERS
 REAL(EB), POINTER, DIMENSION(:,:) :: FLUX0_LS,FLUX1_LS,PHI_LS,PHI1_LS,ROS_BACKU, &
                                      ROS_HEAD,ROS_FLANK,WIND_EXP, &
                                      SR_X_LS,SR_Y_LS,U_LS,V_LS,Z_LS,DZTDX,DZTDY,MAG_ZT, &
-                                     PHI_WS,UMF,THETA_ELPS,PHI_S,PHI_S_X,PHI_S_Y,PHI_W,LS_WORK1,LS_WORK2
+                                     PHI_WS,UMF,THETA_ELPS,PHI_S_X,PHI_S_Y,PHI_W,LS_WORK1,LS_WORK2
 
 
 CONTAINS
@@ -875,7 +875,6 @@ SUBROUTINE POINT_TO_MESH(NM)
 PHI_WS => M%PHI_WS
 UMF => M%UMF
 THETA_ELPS => M%THETA_ELPS
-PHI_S => M%PHI_S
 PHI_S_X => M%PHI_S_X
 PHI_S_Y => M%PHI_S_Y
 PHI_W => M%PHI_W

Source/vege.f90

@@ -156,7 +156,7 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_2(NM,MODE)
 
 INTEGER, INTENT(IN) :: NM,MODE
 INTEGER :: I,IM1,IM2,IIG,IP1,IP2,J,JJG,JM1,JP1
-REAL(EB) :: DZT_DUM,G_EAST,G_WEST,G_SOUTH,G_NORTH
+REAL(EB) :: DZT_DUM,DZTDX_DUM,DZTDY_DUM,G_EAST,G_WEST,G_SOUTH,G_NORTH
 
 T_NOW = CURRENT_TIME()
 
@@ -218,7 +218,6 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_2(NM,MODE)
 ! Rothermel 'Phi' factors for effects of Wind and Slope on ROS
 
 ALLOCATE(M%PHI_WS(IBAR,JBAR))   ; CALL ChkMemErr('VEGE:LEVEL SET','PHI_W',IZERO)   ; PHI_WS => M%PHI_WS    ; PHI_WS = 0.0_EB
-ALLOCATE(M%PHI_S(IBAR,JBAR))    ; CALL ChkMemErr('VEGE:LEVEL SET','PHI_S',IZERO)   ; PHI_S => M%PHI_S
 ALLOCATE(M%PHI_S_X(IBAR,JBAR))  ; CALL ChkMemErr('VEGE:LEVEL SET','PHI_S_X',IZERO) ; PHI_S_X => M%PHI_S_X
 ALLOCATE(M%PHI_S_Y(IBAR,JBAR))  ; CALL ChkMemErr('VEGE:LEVEL SET','PHI_S_Y',IZERO) ; PHI_S_Y => M%PHI_S_Y
 
@@ -296,16 +295,19 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_2(NM,MODE)
          E_ROTH = 0.715_EB * EXP(-0.01094_EB * SF%VEG_LSET_SIGMA)
          BETA_OP_ROTH = 0.20395_EB * (SF%VEG_LSET_SIGMA**(-0.8189_EB))! Optimum packing ratio
 
-         ! Limit effect to slope lte 80 degrees. Phi_s_x,y are slope factors (Rothermel model)
-
-         DZT_DUM = MIN(5.67_EB,ABS(DZTDX(IIG,JJG))) ! 5.67 ~ tan 80 deg
-         PHI_S_X(IIG,JJG) = 5.275_EB * ((SF%VEG_LSET_BETA)**(-0.3_EB)) * DZT_DUM**2
-         PHI_S_X(IIG,JJG) = SIGN(PHI_S_X(IIG,JJG),DZTDX(IIG,JJG))
-         DZT_DUM = MIN(1.73_EB,ABS(DZTDY(IIG,JJG))) ! 1.73 ~ tan 60 deg
-         PHI_S_Y(IIG,JJG) = 5.275_EB * ((SF%VEG_LSET_BETA)**(-0.3_EB)) * DZT_DUM**2
-         PHI_S_Y(IIG,JJG) = SIGN(PHI_S_Y(IIG,JJG),DZTDY(IIG,JJG))
+         ! Slope vector
+         DZTDX_DUM = DZTDX(IIG,JJG)
+         DZTDY_DUM = DZTDY(IIG,JJG)
+         DZT_DUM = SQRT(DZTDX_DUM**2._EB+DZTDY_DUM**2._EB)
+         ! Limit effect to slope lte 80 degrees
+         IF (DZT_DUM>5.67_EB) THEN
+            DZTDX_DUM = DZTDX_DUM * 5.67_EB/DZT_DUM
+            DZTDY_DUM = DZTDY_DUM * 5.67_EB/DZT_DUM
+            DZT_DUM = 5.67_EB
+         ENDIF
 
-         PHI_S(IIG,JJG) = SQRT(PHI_S_X(IIG,JJG)**2 + PHI_S_Y(IIG,JJG)**2)
+         PHI_S_X(IIG,JJG) = 5.275_EB * ((SF%VEG_LSET_BETA)**(-0.3_EB)) * DZTDX_DUM * DZT_DUM
+         PHI_S_Y(IIG,JJG) = 5.275_EB * ((SF%VEG_LSET_BETA)**(-0.3_EB)) * DZTDY_DUM * DZT_DUM
 
       ENDIF IF_ELLIPSE
 
@@ -332,8 +334,8 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
 REAL(EB), INTENT(IN) :: T,DT
 INTEGER :: IIG,IW,JJG,IC,OUTPUT_INDEX
 INTEGER :: KDUM,KWIND,ICF,IKT
-REAL(EB) :: UMF_TMP,PHX,PHY,MAG_PHI,PHI_W_X,PHI_W_Y,UMF_X,UMF_Y,UMAG,ROS_MAG,UMF_MAG,WIND_FACTOR,&
-            SIN_THETA,COS_THETA,THETA,ZWIND(2),U_Z(2),V_Z(2),REF_WIND_HEIGHT
+REAL(EB) :: UMF_TMP,PHX,PHY,MAG_PHI,PHI_S,PHI_W_X,PHI_W_Y,UMF_X,UMF_Y,UMAG,ROS_MAG,UMF_MAG,&
+            WIND_FACTOR,SIN_THETA,COS_THETA,THETA,ZWIND(2),U_Z(2),V_Z(2),REF_WIND_HEIGHT
 
 T_NOW = CURRENT_TIME()
 
@@ -464,7 +466,9 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
 
          ! Include Rothermel slope factor
 
-         IF (PHI_S(IIG,JJG) > 0.0_EB) THEN
+         PHI_S = SQRT(PHI_S_X(IIG,JJG)+PHI_S_Y(IIG,JJG))
+
+         IF (PHI_S > 0.0_EB) THEN
 
             PHX = PHI_W_X + PHI_S_X(IIG,JJG)
             PHY = PHI_W_Y + PHI_S_Y(IIG,JJG)
@@ -487,13 +491,11 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
             ! 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 = (((ABS(PHI_S_X(IIG,JJG)) * (SF%VEG_LSET_BETA / BETA_OP_ROTH)**E_ROTH)/C_ROTH)**(1/B_ROTH))*0.3048
-            UMF_TMP = SIGN(UMF_TMP,PHI_S_X(IIG,JJG))
-            UMF_X = UMF_X + UMF_TMP
+            UMF_TMP = &
+            0.3048_EB/PHI_S*(((SF%VEG_LSET_BETA / BETA_OP_ROTH)**E_ROTH)*PHI_S/C_ROTH)**(1._EB/B_ROTH)
 
-            UMF_TMP = (((ABS(PHI_S_Y(IIG,JJG)) * (SF%VEG_LSET_BETA / BETA_OP_ROTH)**E_ROTH)/C_ROTH)**(1/B_ROTH))*0.3048
-            UMF_TMP = SIGN(UMF_TMP,PHI_S_Y(IIG,JJG))
-            UMF_Y = UMF_Y + UMF_TMP
+            UMF_X = UMF_X + UMF_TMP*PHI_S_X(IIG,JJG)
+            UMF_Y = UMF_Y + UMF_TMP*PHI_S_Y(IIG,JJG)
 
          ELSE
 

Utilities/Matlab/scripts/level_set_ellipse.m

@@ -163,7 +163,7 @@
 end
 
 max_err = max(error_table{:,2:end}(:));
-if max_err>.15
+if max_err>.2
     display(['Matlab Warning: LS_ellipse is out of tolerance. Max error = ',num2str(max_err)])
 end
 
@@ -218,8 +218,9 @@
 function [phi_s_x,phi_s_y] = slope_adj(dzdx,dzdy)
     beta=0.01; % packing ratio (-), FDS default
     % Rothermel optimum packing ratio
-    phi_s_x = sign(dzdx)*5.275*beta^-0.3*dzdx^2;
-    phi_s_y = sign(dzdy)*5.275*beta^-0.3*dzdy^2;
+    dzds = sqrt(dzdx^2+dzdy^2);
+    phi_s_x = 5.275*beta^-0.3*dzdx*dzds;
+    phi_s_y = 5.275*beta^-0.3*dzdy*dzds;
 end
 
 % calculate virtual wind vector for slope
@@ -233,7 +234,13 @@
     C = 7.47*exp(-0.8711*sigma.^0.55);
     B = 0.15988*sigma.^0.54;
     E = 0.715*exp(-0.01094*sigma);
+    phi_s=sqrt(phi_s_x^2+phi_s_y^2);
+    if (phi_s>0)
+        uv_tmp = 0.3048/phi_s*(phi_s/C*beta_ratio^E)^(1/B);
+    else
+        uv_tmp = 0;
+    end
     % divide by 60 to maintain units of m/s elsewhere
-    u_virtual = 0.3048*(phi_s_x/C*beta_ratio^E)^(1/B)/60;
-    v_virtual = 0.3048*(phi_s_y/C*beta_ratio^E)^(1/B)/60;
+    u_virtual = uv_tmp*phi_s_x/60;
+    v_virtual = uv_tmp*phi_s_y/60;
 end