Skip to content

Commit e93c161

Browse files
authored
Merge pull request #14666 from ericvmueller/ls_wind
Add VEG_LSET_WIND_HEIGHT and VEG_LSET_WIND_RAMP for user-defined wind speed functions in level set spread
2 parents 42e0d6d + 93e62c5 commit e93c161

File tree

11 files changed

+246
-26
lines changed

11 files changed

+246
-26
lines changed

Manuals/FDS_User_Guide/FDS_User_Guide.tex

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7741,6 +7741,26 @@ \section{Level Set Model for Wildland Fire Spread}
77417741
For all models, the total mineral content $S_{\rm t}=0.056$, the effective mineral content $S_{\rm e}=0.01$, the heat of combustion $\Delta H=18600$~kW/kg, density $\rho_{\rm p}=510$~kg/m$^3$, moisture content of fine, medium, and large dead vegetation, $M_{\rm d,1}=0.03$, $M_{\rm d,2}=0.04$, $M_{\rm d,3}=0.05$, moisture content of live woody and herbaceous vegetation $M_{\rm l,w}=M_{\rm l,h}= 0.70$.
77427742
\end{sidewaystable}
77437743

7744+
\subsection{Custom Wind Speed Function}
7745+
The default Rothermel-Albini formulation~\cite{Rothermel:1972,Albini:1976} represents the influence of wind on the head fire spread rate as
7746+
\be
7747+
R = R_{\rm 0} \, (1+\phi_{\rm W})
7748+
\ee
7749+
where the no-wind, no-slope spread rate, $R_{\rm 0}$, is increased by a factor $\phi_W$, which is a function of the wind speed. In the default model the wind speed is taken as the value at mid-flame height, $U_{\rm mf}$.
7750+
7751+
This introduces two complications. The first is determining the functional form of $\phi_W$. The default formulation is an empirical function which depends on fuel bed properties. However, this has not been rigorously validated for all possible fuel configurations - the original fit was determined with just three fuel bed types~\cite{Rothermel:1972}. In addition, the relationship assumes a reference wind which is unaffected by the fire, which can be contradictory to the approach of \ct{LEVEL_SET_MODE=4}. Therefore, it may be preferable to specify a custom function $\phi_W$ to better match observational data for a given fuel type. This can be done by specifying a \ct{RAMP} using \ct{VEG_LSET_WIND_RAMP}:
7752+
\begin{lstlisting}
7753+
&SURF ..., VEG_LSET_WIND_RAMP='WR' /
7754+
7755+
&RAMP ID='WR', Z = 0, F=0 /
7756+
&RAMP ID='WR', Z = 10, F=4 /
7757+
\end{lstlisting}
7758+
Here, the variable \ct{Z} represents the reference wind speed in m/s and \ct{F} is the value of $\phi_W$. In this example the spread rate increases linearly with wind speed between 0~m/s and 10~m/s, but any functional form can be represented in a piecwise manner using the \ct{RAMP}.
7759+
7760+
The second complication is determining the appropriate mid-flame wind speed. In the default formulation the wind is taken from a reference height of 6.1~m (based on a standard weather station measurement height of 20~ft). This is then re-scaled to an assumed mid-flame height (one fuel-bed height above the fuel) using an assumed logarithmic profile of the wind. This is the so-called wind adjustment factor~\cite{Andrews:2012,Bova:IJWF2015}. However, these assumptions may begin to break down, for example, when the wind profile is not logarithmic. Therefore, it is possible to specify a reference wind height via \ct{VEG_LSET_WIND_HEIGHT} on the \ct{SURF} line. The wind sampled at this height will be directly used in the $\phi_{\rm W}$ function, without any adjustment, whether using the default model or a custom \ct{RAMP}. Note that it is possible to include an adjustment factor when using \ct{VEG_LSET_WIND_HEIGHT} by implicitly including this factor when specifying a \ct{VEG_LSET_WIND_RAMP}.
7761+
7762+
7763+
77447764
\subsection{Level Set Vegetation Drag}
77457765
When a surface is defined using level set fuel parameters, a drag force is applied in the gas phase grid cell adjacent to the boundary following the same approach described for the Boundary Fuel Model in Section~\ref{info:boundary_fuel_model}. In this case, the fuel surface-area to volume ratio (\ct{VEG_LSET_SIGMA}), packing ratio (\ct{VEG_LSET_BETA}), and depth (\ct{VEG_LSET_HT}) are either obtained from the default fuel models or specified by the user for a custom fuel. The shape factor (\ct{SHAPE_FACTOR}) and drag coefficient (\ct{DRAG_COEFFICIENT}) can also be modified from their default values of 0.25 and 2.8, respectively. A comparison of the consistency of the level set drag with the particle and boundary models is included in Figure~\ref{ground_vegetation_drag_fig}.
77467766

@@ -11480,6 +11500,7 @@ \section{Gas Phase Output Quantities}
1148011500
\ct{TOTAL MASS FLUX X}$^1$ & Section~\ref{info:mass_flow} & kg/s/m$^2$ & D,I,P,S \\ \hline
1148111501
\ct{TOTAL MASS FLUX Y}$^1$ & Section~\ref{info:mass_flow} & kg/s/m$^2$ & D,I,P,S \\ \hline
1148211502
\ct{TOTAL MASS FLUX Z}$^1$ & Section~\ref{info:mass_flow} & kg/s/m$^2$ & D,I,P,S \\ \hline
11503+
\ct{U_LS, V_LS} & Velocity for level set spread~\ref{info:level_set}& m/s & D,I,P,S \\ \hline
1148311504
\ct{U-VELOCITY} & Gas velocity component, $u$ & m/s & D,I,P,S \\ \hline
1148411505
\ct{V-VELOCITY} & Gas velocity component, $v$ & m/s & D,I,P,S \\ \hline
1148511506
\ct{W-VELOCITY} & Gas velocity component, $w$ & m/s & D,I,P,S \\ \hline
@@ -13603,6 +13624,8 @@ \section{\texorpdfstring{{\tt SURF}}{SURF} (Surface Properties)}
1360313624
\ct{VEG_LSET_SURF_LOAD} & Real & Section~\ref{info:level_set} & kg/m$^2$ & 1.0 \\ \hline
1360413625
\ct{VEG_LSET_TAN2} & Real & Section~\ref{info:level_set} & & \\ \hline
1360513626
\ct{VEG_LSET_WIND_EXP} & Real & Section~\ref{info:level_set} & & 1. \\ \hline
13627+
\ct{VEG_LSET_WIND_HEIGHT} & Real & Section~\ref{info:level_set} & m & 6.1 \\ \hline
13628+
\ct{VEG_LSET_WIND_RAMP} & Character & Section~\ref{info:level_set} & & \\ \hline
1360613629
\ct{VEL} & Real & Section~\ref{info:Velocity_BC} & m/s & \\ \hline
1360713630
\ct{VEL_BULK} & Real & Section~\ref{info:profiles} & m/s & \\ \hline
1360813631
\ct{VEL_GRAD} & Real & Section~\ref{info:vel_grad} & 1/s & \\ \hline

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7116,6 +7116,35 @@ \subsection{Elliptical Fire Spread}
71167116
\label{fig:level_set_ellipse_error}
71177117
\end{figure}
71187118

7119+
\subsection{Custom Wind Fuction (\texorpdfstring{\ct{LS\_wind\_ramp}}{LS\_wind\_ramp})}
7120+
\label{LS_wind_ramp}
7121+
7122+
The level set model determines a wind-aided head fire spread rate based on the function
7123+
\be
7124+
R = R_{\rm 0} \, (1+\phi_{\rm W})
7125+
\ee
7126+
where the no-wind, no-slope spread rate, $R_{\rm 0}$, is increased by a factor $\phi_W$, which is a function of the wind speed. Note that there may be an additional term accounting for slope, $\phi_{\rm S}$, but this is excluded for simplicity in this example. The ability to customize the form of $\phi_W$ as a function of the reference wind speed, using \ct{VEG_LSET_WIND_RAMP}, is tested here. A line fire is ignited on a level set surface with a no-wind, no-slope spread rate, \ct{VEG_LSET_ROS_00}, of 0.03~m/s. The wind speed, using a flat wind profile, is progressively increased from 0~m/s to 5~m/s. In the first case,\ct{LS_wind_ramp_lin.fds}, the wind ramp is set to reproduce a linear function
7127+
\be
7128+
R = R_{\rm 0} \, (1+0.4 \, U)
7129+
\ee
7130+
and in the second case, \ct{LS_wind_ramp_quad.fds}, it is set to match a quadratic function
7131+
\be
7132+
R = R_{\rm 0} \, (1+ 0.1 \, U^2)
7133+
\ee
7134+
The maximum spread rate as a function of wind speed is plotted in Figure~\ref{fig:LS_wind_ramp} to confirm that the functions are properly applied.
7135+
7136+
7137+
\begin{figure}[h]
7138+
\begin{center}
7139+
\begin{tabular}{cc}
7140+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/LS_wind_ramp.pdf}
7141+
\end{tabular}
7142+
\end{center}
7143+
\caption[Level set wind ramp]{Expected versus predicted spread rate as a function of linear and quadratic wind speed functions.}
7144+
\label{fig:LS_wind_ramp}
7145+
\end{figure}
7146+
7147+
71197148
\subsection{Ember Generation (\texorpdfstring{\ct{LS4\_ember\_yield}}{LS4\_ember\_yield})}
71207149
\label{LS4_ember_yield}
71217150

Source/data.f90

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1304,6 +1304,14 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
13041304
OUTPUT_QUANTITY(551)%UNITS = '1/s'
13051305
OUTPUT_QUANTITY(551)%SHORT_NAME = 'cart_div'
13061306

1307+
OUTPUT_QUANTITY(552)%NAME = 'U_LS'
1308+
OUTPUT_QUANTITY(552)%UNITS = 'm/s'
1309+
OUTPUT_QUANTITY(552)%SHORT_NAME = 'U_LS'
1310+
1311+
OUTPUT_QUANTITY(553)%NAME = 'V_LS'
1312+
OUTPUT_QUANTITY(553)%UNITS = 'm/s'
1313+
OUTPUT_QUANTITY(553)%SHORT_NAME = 'V_LS'
1314+
13071315
! Boundary Quantities (Negative indices)
13081316

13091317
OUTPUT_QUANTITY(-1)%NAME = 'RADIATIVE HEAT FLUX'

Source/dump.f90

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -8510,11 +8510,18 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
85108510
(ADV_FZ(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FZ(II,JJ,KK,1:N_TRACKED_SPECIES)) )
85118511
ENDIF
85128512
ENDIF
8513-
CASE(550) ! CUTCELL VELOCITY DIVERGENCE
8514-
GAS_PHASE_OUTPUT_RES = CARTVELDIV(II,JJ,KK)
85158513

8516-
CASE(551) ! CARTESIAN VELOCITY DIVERGENCE
8517-
GAS_PHASE_OUTPUT_RES = CARTVELDIV(II,JJ,KK)
8514+
CASE(550) ! CUTCELL VELOCITY DIVERGENCE
8515+
GAS_PHASE_OUTPUT_RES = CARTVELDIV(II,JJ,KK)
8516+
8517+
CASE(551) ! CARTESIAN VELOCITY DIVERGENCE
8518+
GAS_PHASE_OUTPUT_RES = CARTVELDIV(II,JJ,KK)
8519+
8520+
CASE(552) ! U_LS
8521+
GAS_PHASE_OUTPUT_RES = U_LS(II,JJ)
8522+
8523+
CASE(553) ! V_LS
8524+
GAS_PHASE_OUTPUT_RES = V_LS(II,JJ)
85188525

85198526
END SELECT IND_SELECT
85208527

Source/read.f90

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7775,7 +7775,8 @@ SUBROUTINE READ_SURF(QUICK_READ)
77757775
MATL_ID(MAX_LAYERS,MAX_MATERIALS),PROFILE,BACKING,GEOMETRY,RAMP_EF,RAMP_PART,NAME_LIST(MAX_MATERIALS*MAX_LAYERS),&
77767776
SPEC_ID(MAX_SPECIES),RAMP_TMP_BACK,RAMP_TMP_GAS_BACK,RAMP_TMP_GAS_FRONT,&
77777777
RAMP_V_X,RAMP_V_Y,RAMP_V_Z,NEAR_WALL_TURBULENCE_MODEL,SURF_DEFAULT,&
7778-
RAMP_HEAT_TRANSFER_COEFFICIENT,RAMP_HEAT_TRANSFER_COEFFICIENT_BACK,RAMP_IHS(MAX_LAYERS)
7778+
RAMP_HEAT_TRANSFER_COEFFICIENT,RAMP_HEAT_TRANSFER_COEFFICIENT_BACK,RAMP_IHS(MAX_LAYERS),&
7779+
VEG_LSET_WIND_RAMP
77797780
CHARACTER(LABEL_LENGTH), DIMENSION(10) :: INIT_IDS
77807781
LOGICAL :: ADIABATIC,BURN_AWAY,FREE_SLIP,NO_SLIP,CONVERT_VOLUME_TO_MASS,HORIZONTAL,DIRICHLET_FRONT,DIRICHLET_BACK, BLOWING, &
77817782
INERT_Q_REF,ALLOW_UNDERSIDE_PARTICLES,ALLOW_SURFACE_PARTICLES
@@ -7806,7 +7807,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
78067807
REAL(EB) :: VEG_LSET_IGNITE_TIME,VEG_LSET_QCON,VEG_LSET_ROS_HEAD,VEG_LSET_ROS_FLANK,VEG_LSET_ROS_BACK, &
78077808
VEG_LSET_WIND_EXP,VEG_LSET_BETA,VEG_LSET_HT,VEG_LSET_SIGMA,VEG_LSET_ROS_00, &
78087809
VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_SURF_LOAD,VEG_LSET_FIREBASE_TIME,&
7809-
VEG_LSET_CHAR_FRACTION,VEL_PART,INIT_PER_AREA,TIME_STEP_FACTOR
7810+
VEG_LSET_CHAR_FRACTION,VEG_LSET_WIND_HEIGHT,VEL_PART,INIT_PER_AREA,TIME_STEP_FACTOR
78107811
LOGICAL :: DEFAULT,VEG_LSET_SPREAD,VEG_LSET_TAN2,TGA_ANALYSIS,VARIABLE_THICKNESS,HT3D,THERM_THICK,VEG_LSET_ROS_FIXED
78117812
LOGICAL, ALLOCATABLE, DIMENSION(:) :: DUPLICATE
78127813
! Ember generating variables
@@ -7843,7 +7844,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
78437844
VEG_LSET_BETA,VEG_LSET_CHAR_FRACTION,VEG_LSET_FIREBASE_TIME,VEG_LSET_FUEL_INDEX,VEG_LSET_HT,VEG_LSET_IGNITE_TIME,&
78447845
VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_QCON,&
78457846
VEG_LSET_ROS_00,VEG_LSET_ROS_BACK,VEG_LSET_ROS_FLANK,VEG_LSET_ROS_HEAD,VEG_LSET_ROS_FIXED,VEG_LSET_SIGMA,&
7846-
VEG_LSET_SURF_LOAD,VEG_LSET_TAN2,VEG_LSET_WIND_EXP,&
7847+
VEG_LSET_SURF_LOAD,VEG_LSET_TAN2,VEG_LSET_WIND_EXP,VEG_LSET_WIND_RAMP,VEG_LSET_WIND_HEIGHT,&
78477848
VEL,VEL_BULK,VEL_GRAD,VEL_PART,VEL_T,VOLUME_FLOW,WIDTH,XYZ,Z0,Z_0
78487849

78497850
! Count the SURF lines in the input file
@@ -8147,6 +8148,9 @@ SUBROUTINE READ_SURF(QUICK_READ)
81478148
SF%VEG_LSET_CHAR_FRACTION= VEG_LSET_CHAR_FRACTION
81488149
SF%VEG_LSET_FIREBASE_TIME= VEG_LSET_FIREBASE_TIME
81498150
IF (VEG_LSET_FIREBASE_TIME<0._EB .AND. VEG_LSET_SIGMA>0._EB) SF%VEG_LSET_FIREBASE_TIME = 75600._EB/VEG_LSET_SIGMA
8151+
IF (VEG_LSET_WIND_RAMP/='null') CALL GET_RAMP_INDEX(VEG_LSET_WIND_RAMP,'PROFILE',SF%I_RAMP_LS_WIND)
8152+
SF%VEG_LSET_WIND_HEIGHT = VEG_LSET_WIND_HEIGHT
8153+
81508154

81518155
IF (SF%VEG_LSET_FUEL_INDEX>0 .AND. COLOR=='null' .AND. ANY(RGB<0)) THEN
81528156
SELECT CASE(SF%VEG_LSET_FUEL_INDEX)
@@ -9080,6 +9084,8 @@ SUBROUTINE SET_SURF_DEFAULTS
90809084
VEG_LSET_SURF_LOAD = -1.0_EB !kg/m^2
90819085
VEG_LSET_FIREBASE_TIME = -1.0_EB
90829086
VEG_LSET_CHAR_FRACTION = 0.20_EB
9087+
VEG_LSET_WIND_RAMP = 'null'
9088+
VEG_LSET_WIND_HEIGHT = -1.0_EB
90839089

90849090
END SUBROUTINE SET_SURF_DEFAULTS
90859091

Source/type.f90

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1001,8 +1001,8 @@ MODULE TYPES
10011001
REAL(EB) :: VEG_LSET_IGNITE_T,VEG_LSET_ROS_HEAD,VEG_LSET_ROS_00,VEG_LSET_QCON,VEG_LSET_ROS_FLANK,VEG_LSET_ROS_BACK, &
10021002
VEG_LSET_WIND_EXP,VEG_LSET_SIGMA,VEG_LSET_HT,VEG_LSET_BETA,&
10031003
VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_SURF_LOAD,VEG_LSET_FIREBASE_TIME, &
1004-
VEG_LSET_CHAR_FRACTION
1005-
INTEGER :: VEG_LSET_FUEL_INDEX
1004+
VEG_LSET_CHAR_FRACTION,VEG_LSET_WIND_HEIGHT
1005+
INTEGER :: VEG_LSET_FUEL_INDEX,I_RAMP_LS_WIND=-1
10061006

10071007
END TYPE SURFACE_TYPE
10081008

Source/vege.f90

Lines changed: 30 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -332,8 +332,8 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
332332
REAL(EB), INTENT(IN) :: T,DT
333333
INTEGER :: IIG,IW,JJG,IC,OUTPUT_INDEX
334334
INTEGER :: KDUM,KWIND,ICF,IKT
335-
REAL(EB) :: UMF_TMP,PHX,PHY,MAG_PHI,PHI_W_X,PHI_W_Y,UMF_X,UMF_Y,UMAG,ROS_MAG,UMF_MAG,ROTH_FACTOR,&
336-
SIN_THETA,COS_THETA,THETA,ZWIND(2),U_Z(2),V_Z(2)
335+
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,&
336+
SIN_THETA,COS_THETA,THETA,ZWIND(2),U_Z(2),V_Z(2),REF_WIND_HEIGHT
337337

338338
T_NOW = CURRENT_TIME()
339339

@@ -362,6 +362,10 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
362362

363363
! Establish the wind field
364364

365+
REF_WIND_HEIGHT = SF%VEG_LSET_WIND_HEIGHT
366+
! If not set, use Behave/Andrews approach
367+
IF (SF%VEG_LSET_WIND_HEIGHT<0._EB) REF_WIND_HEIGHT = 6.1_EB
368+
365369
IF_CFD_COUPLED: IF (LEVEL_SET_COUPLED_WIND .AND. .NOT. LEVEL_SET_MODE==5) THEN ! The wind speed is derived from the CFD
366370

367371
U_LS(IIG,JJG) = 0.5_EB*(U(IIG-1,JJG,K_LS(IIG,JJG))+U(IIG,JJG,K_LS(IIG,JJG)))
@@ -372,20 +376,23 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
372376
! Evaluate time and height varying profiles, using 6.1 m reference height
373377
IF (I_RAMP_DIRECTION_T/=0 .OR. I_RAMP_DIRECTION_Z/=0) THEN
374378
IF (I_RAMP_DIRECTION_T==0) THEN
375-
THETA=EVALUATE_RAMP(6.1_EB,I_RAMP_DIRECTION_Z)*DEG2RAD
379+
THETA=EVALUATE_RAMP(REF_WIND_HEIGHT,I_RAMP_DIRECTION_Z)*DEG2RAD
376380
ELSEIF (I_RAMP_DIRECTION_Z==0) THEN
377381
THETA=EVALUATE_RAMP(T,I_RAMP_DIRECTION_T)*DEG2RAD
378382
ELSE
379-
THETA=(EVALUATE_RAMP(6.1_EB,I_RAMP_DIRECTION_Z)+EVALUATE_RAMP(T,I_RAMP_DIRECTION_T))*DEG2RAD
383+
THETA=(EVALUATE_RAMP(REF_WIND_HEIGHT,I_RAMP_DIRECTION_Z)+&
384+
EVALUATE_RAMP(T,I_RAMP_DIRECTION_T))*DEG2RAD
380385
ENDIF
381386
SIN_THETA = -SIN(THETA)
382387
COS_THETA = -COS(THETA)
383388
ELSE
384389
SIN_THETA = 1._EB
385390
COS_THETA = 1._EB
386391
ENDIF
387-
U_LS(IIG,JJG) = U0*EVALUATE_RAMP(6.1_EB,I_RAMP_SPEED_Z)*EVALUATE_RAMP(T,I_RAMP_SPEED_T)*SIN_THETA
388-
V_LS(IIG,JJG) = V0*EVALUATE_RAMP(6.1_EB,I_RAMP_SPEED_Z)*EVALUATE_RAMP(T,I_RAMP_SPEED_T)*COS_THETA
392+
U_LS(IIG,JJG) = U0*EVALUATE_RAMP(REF_WIND_HEIGHT,I_RAMP_SPEED_Z)*&
393+
EVALUATE_RAMP(T,I_RAMP_SPEED_T)*SIN_THETA
394+
V_LS(IIG,JJG) = V0*EVALUATE_RAMP(REF_WIND_HEIGHT,I_RAMP_SPEED_Z)*&
395+
EVALUATE_RAMP(T,I_RAMP_SPEED_T)*COS_THETA
389396

390397
ENDIF IF_CFD_COUPLED
391398

@@ -396,14 +403,14 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
396403
IF (LEVEL_SET_COUPLED_WIND) THEN
397404

398405
! Check if sample height is in the mesh
399-
IF (ZC(KBAR)<6.1_EB) THEN
406+
IF (ZC(KBAR)<REF_WIND_HEIGHT) THEN
400407
U_LS(IIG,JJG) = 0.5_EB*(U(IIG-1,JJG,KBAR)+U(IIG,JJG,KBAR))
401408
V_LS(IIG,JJG) = 0.5_EB*(V(IIG,JJG-1,KBAR)+V(IIG,JJG,KBAR))
402409
ELSE
403410

404411
KWIND = 0
405412
DO KDUM = K_LS(IIG,JJG),KBAR
406-
IF ((ZC(KDUM)-Z_LS(IIG,JJG))>=6.1_EB) THEN
413+
IF ((ZC(KDUM)-Z_LS(IIG,JJG))>=REF_WIND_HEIGHT) THEN
407414
KWIND = KDUM
408415
ZWIND(1) = ZC(KWIND-1)-Z_LS(IIG,JJG)
409416
ZWIND(2) = ZC(KWIND)-Z_LS(IIG,JJG)
@@ -419,30 +426,37 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
419426
IF (KWIND==1) THEN
420427
U_Z(1) = 0._EB; V_Z(1) = 0._EB; ZWIND(1) = 0._EB
421428
ENDIF
422-
U_LS(IIG,JJG) = U_Z(1) + (6.1_EB-ZWIND(1))/(ZWIND(2)-ZWIND(1))*(U_Z(2)-U_Z(1))
423-
V_LS(IIG,JJG) = V_Z(1) + (6.1_EB-ZWIND(1))/(ZWIND(2)-ZWIND(1))*(V_Z(2)-V_Z(1))
429+
U_LS(IIG,JJG) = U_Z(1) + (REF_WIND_HEIGHT-ZWIND(1))/(ZWIND(2)-ZWIND(1))*(U_Z(2)-U_Z(1))
430+
V_LS(IIG,JJG) = V_Z(1) + (REF_WIND_HEIGHT-ZWIND(1))/(ZWIND(2)-ZWIND(1))*(V_Z(2)-V_Z(1))
424431

425432
ENDIF
426433

427434
ENDIF
428435

429-
! Wind at midflame height (UMF). From Andrews 2012, USDA FS Gen Tech Rep. RMRS-GTR-266 (with added SI conversion)
436+
UMF_TMP = 1._EB
430437

431-
UMF_TMP = 1.83_EB / LOG((20.0_EB + 1.18_EB * SF%VEG_LSET_HT) /(0.43_EB * SF%VEG_LSET_HT)) ! Bova et al., Eq. A2
438+
! Wind at midflame height (UMF). From Andrews 2012, USDA FS Gen Tech Rep. RMRS-GTR-266 (with added SI conversion)
439+
IF (SF%VEG_LSET_WIND_HEIGHT<0._EB) &
440+
UMF_TMP = 1.83_EB / LOG((20.0_EB + 1.18_EB * SF%VEG_LSET_HT) /(0.43_EB * SF%VEG_LSET_HT)) ! Bova et al., Eq. A2
432441

433442
! Factor 60 converts U from m/s to m/min which is used in elliptical model.
434443

435444
UMF_X = UMF_TMP * U_LS(IIG,JJG) * 60.0_EB
436445
UMF_Y = UMF_TMP * V_LS(IIG,JJG) * 60.0_EB
437446
UMF_MAG = SQRT(UMF_X**2 + UMF_Y**2)
438447

439-
! Components of Rothermel wind factor - affects spread rate
448+
! Compute wind factor affecting spread rate R(U) = R_0*(1+WIND_FACTOR)
449+
450+
IF (SF%I_RAMP_LS_WIND>0) THEN
451+
WIND_FACTOR = EVALUATE_RAMP(UMF_MAG/60._EB,SF%I_RAMP_LS_WIND)
452+
ELSE
453+
WIND_FACTOR = C_ROTH * ((3.281_EB * UMF_MAG)**B_ROTH) * (SF%VEG_LSET_BETA/BETA_OP_ROTH)**(-E_ROTH) ! Bova et al., Eq. A1
454+
ENDIF
440455

441-
ROTH_FACTOR = C_ROTH * ((3.281_EB * UMF_MAG)**B_ROTH) * (SF%VEG_LSET_BETA / BETA_OP_ROTH)**(-E_ROTH) ! Bova et al., Eq. A1
442456

443457
IF (UMF_MAG>TWO_EPSILON_EB) THEN
444-
PHI_W_X = ROTH_FACTOR*UMF_X/UMF_MAG
445-
PHI_W_Y = ROTH_FACTOR*UMF_Y/UMF_MAG
458+
PHI_W_X = WIND_FACTOR*UMF_X/UMF_MAG
459+
PHI_W_Y = WIND_FACTOR*UMF_Y/UMF_MAG
446460
ELSE
447461
PHI_W_X = 0.0_EB
448462
PHI_W_Y = 0.0_EB

0 commit comments

Comments
 (0)