diff --git a/Source/mesh.f90 b/Source/mesh.f90 index c2219812691..85d97a67d4e 100644 --- a/Source/mesh.f90 +++ b/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 diff --git a/Source/vege.f90 b/Source/vege.f90 index af0d1e47602..9e76fe42c8d 100644 --- a/Source/vege.f90 +++ b/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 diff --git a/Utilities/Matlab/scripts/level_set_ellipse.m b/Utilities/Matlab/scripts/level_set_ellipse.m index 43adfcc482c..f33f883ca67 100644 --- a/Utilities/Matlab/scripts/level_set_ellipse.m +++ b/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