diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 537e74824d..93520c4c56 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -2124,6 +2124,12 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft) new_patch%fabi_sha_z(:,:,:) = 0._r8 new_patch%scorch_ht(:) = 0._r8 new_patch%frac_burnt = 0._r8 + new_patch%canopy_bulk_density = 0._r8 + new_patch%canopy_fuel_load = 0._r8 + new_patch%passive_crown_FI = 0._r8 + new_patch%ros_torch = 0._r8 + new_patch%heat_per_area = 0._r8 + new_patch%lb = 0._r8 new_patch%litter_moisture(:) = 0._r8 new_patch%fuel_eff_moist = 0._r8 new_patch%livegrass = 0._r8 @@ -2235,11 +2241,18 @@ subroutine zero_patch(cp_p) currentPatch%fi = nan ! average fire intensity of flaming front during day. ! backward ros plays no role. kj/m/s or kw/m. currentPatch%fire = 999 ! sr decide_fire.1=fire hot enough to proceed. 0=stop everything- no fires today + currentPatch%active_crown_fire_flg = 9999 ! flag to indicate active crown fire ignition currentPatch%fd = nan ! fire duration (mins) currentPatch%ros_back = nan ! backward ros (m/min) currentPatch%scorch_ht(:) = nan ! scorch height of flames on a given PFT - currentPatch%frac_burnt = nan ! fraction burnt daily - currentPatch%burnt_frac_litter(:) = nan + currentPatch%frac_burnt = nan ! fraction burnt daily + currentPatch%canopy_bulk_density = nan ! available canopy fuel bulk density in patch (kg biomass/m3) + currentPatch%canopy_fuel_load = nan ! available canopy fuel load in patch (kg biomass) + currentPatch%passive_crown_FI = nan ! fire intensity for ignition of passive canopy fuel (kW/m) + currentPatch%ros_torch = nan ! rate of spread (ROS) for crown torch initiation (m/min) + currentPatch%heat_per_area = nan ! heat release per unit area (kJ/m2) for surface fuel + currentPatch%lb = nan ! length to breadth ratio of fire ellipse (unitless) + currentPatch%burnt_frac_litter(:) = nan currentPatch%btran_ft(:) = 0.0_r8 currentPatch%canopy_layer_tlai(:) = 0.0_r8 @@ -2658,6 +2671,12 @@ subroutine fuse_2_patches(csite, dp, rp) rp%ros_back = (dp%ros_back*dp%area + rp%ros_back*rp%area) * inv_sum_area rp%scorch_ht(:) = (dp%scorch_ht(:)*dp%area + rp%scorch_ht(:)*rp%area) * inv_sum_area rp%frac_burnt = (dp%frac_burnt*dp%area + rp%frac_burnt*rp%area) * inv_sum_area + rp%canopy_bulk_density = (dp%canopy_bulk_density*dp%area + rp%canopy_bulk_density*rp%area) * inv_sum_area + rp%canopy_fuel_load = (dp%canopy_fuel_load*dp%area + rp%canopy_fuel_load*rp%area) * inv_sum_area + rp%passive_crown_FI = (dp%passive_crown_FI*dp%area + rp%passive_crown_FI*rp%area) * inv_sum_area + rp%ros_torch = (dp%ros_torch*dp%area + rp%ros_torch*rp%area) * inv_sum_area + rp%heat_per_area = (dp%heat_per_area*dp%area + rp%heat_per_area*rp%area) * inv_sum_area + rp%lb = (dp%lb*dp%area + rp%lb*rp%area) * inv_sum_area rp%burnt_frac_litter(:) = (dp%burnt_frac_litter(:)*dp%area + rp%burnt_frac_litter(:)*rp%area) * inv_sum_area rp%btran_ft(:) = (dp%btran_ft(:)*dp%area + rp%btran_ft(:)*rp%area) * inv_sum_area rp%zstar = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index aedcb4aa7c..fd1d630a29 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -27,6 +27,7 @@ module SFMainMod use EDtypesMod , only : ed_cohort_type use EDtypesMod , only : AREA use EDtypesMod , only : DL_SF +! use EDtypesMod , only : crown_fire_threshold ! TODO slevis: NOT used; if end up using, look at SF_val_fire_threshold as template use EDTypesMod , only : TW_SF use EDtypesMod , only : LB_SF use EDtypesMod , only : LG_SF @@ -35,7 +36,6 @@ module SFMainMod use EDtypesMod , only : TR_SF use FatesLitterMod , only : litter_type - use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : all_carbon_elements use PRTGenericMod, only : leaf_organ @@ -53,11 +53,13 @@ module SFMainMod public :: fire_model public :: fire_danger_index - public :: charecteristics_of_fuel + public :: characteristics_of_fuel + public :: characteristics_of_crown public :: rate_of_spread public :: ground_fuel_consumption public :: wind_effect public :: area_burnt_intensity + public :: active_crown_fire public :: crown_scorching public :: crown_damage public :: cambial_damage_kill @@ -87,11 +89,15 @@ subroutine fire_model( currentSite, bc_in) type (ed_patch_type), pointer :: currentPatch + !zero fire things currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) currentPatch%frac_burnt = 0.0_r8 + currentPatch%FI = 0.0_r8 + currentPatch%FD = 0.0_r8 currentPatch%fire = 0 + currentPatch%active_crown_fire_flg = 0 currentPatch => currentPatch%older enddo @@ -102,10 +108,12 @@ subroutine fire_model( currentSite, bc_in) if( hlm_spitfire_mode > hlm_sf_nofire_def )then call fire_danger_index(currentSite, bc_in) call wind_effect(currentSite, bc_in) - call charecteristics_of_fuel(currentSite) + call characteristics_of_fuel(currentSite) + call characteristics_of_crown(currentSite) call rate_of_spread(currentSite) call ground_fuel_consumption(currentSite) call area_burnt_intensity(currentSite, bc_in) + call active_crown_fire (currentSite) call crown_scorching(currentSite) call crown_damage(currentSite) call cambial_damage_kill(currentSite) @@ -153,21 +161,22 @@ subroutine fire_danger_index ( currentSite, bc_in) else yipsolon = (SF_val_fdi_a* temp_in_C)/(SF_val_fdi_b+ temp_in_C)+log(max(1.0_r8,rh)/100.0_r8) dewpoint = (SF_val_fdi_b*yipsolon)/(SF_val_fdi_a-yipsolon) !Standard met. formula - d_NI = ( temp_in_C-dewpoint)* temp_in_C !follows Nesterov 1968. Equation 5. Thonicke et al. 2010. + d_NI = ( temp_in_C-dewpoint)* temp_in_C !follows Nesterov 1968. Eq 5, Thonicke et al. 2010. if (d_NI < 0.0_r8) then !Change in NI cannot be negative. d_NI = 0.0_r8 !check endif endif - currentSite%acc_NI = currentSite%acc_NI + d_NI !Accumulate Nesterov index over the fire season. + currentSite%acc_NI = currentSite%acc_NI + d_NI !Accumulate Nesterov index over fire season. end subroutine fire_danger_index !***************************************************************** - subroutine charecteristics_of_fuel ( currentSite ) + subroutine characteristics_of_fuel ( currentSite) !***************************************************************** - use SFParamsMod, only : SF_val_drying_ratio, SF_val_SAV, SF_val_FBD + use SFParamsMod, only: SF_val_drying_ratio, SF_val_SAV, SF_val_FBD, & + SF_val_miner_total type(ed_site_type), intent(in), target :: currentSite @@ -176,8 +185,8 @@ subroutine charecteristics_of_fuel ( currentSite ) type(litter_type), pointer :: litt_c real(r8) alpha_FMC(nfsc) ! Relative fuel moisture adjusted per drying ratio - real(r8) fuel_moisture(nfsc) ! Scaled moisture content of small litter fuels. - real(r8) MEF(nfsc) ! Moisture extinction factor of fuels integer n + real(r8) fuel_moisture(nfsc) ! Scaled moisture content of small litter fuels + real(r8) MEF(nfsc) ! Moisture extinction factor of fuels, integer n fuel_moisture(:) = 0.0_r8 @@ -237,9 +246,9 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_frac(lg_sf) = currentPatch%livegrass / currentPatch%sum_fuel ! MEF (moisure of extinction) depends on compactness of fuel, depth, particle size, wind, slope - ! Eqn here is eqn 27 from Peterson and Ryan (1986) "Modeling Postfire Conifer Mortality for Long-Range Planning" + ! Eq here is Eq 27 from Peterson and Ryan (1986) "Modeling Postfire Conifer Mortality for Long-Range Planning" ! but lots of other approaches in use out there... - ! MEF: pine needles=0.30 (text near EQ 28 Rothermal 1972) + ! MEF: pine needles=0.30 (text near Eq 28 Rothermal 1972) ! Table II-1 NFFL mixed fuels models from Rothermal 1983 Gen. Tech. Rep. INT-143 ! MEF: short grass=0.12,tall grass=0.25,chaparral=0.20,closed timber litter=0.30,hardwood litter=0.25 ! Thonicke 2010 SAV values propagated thru P&R86 eqn below gives MEF:tw=0.355, sb=0.44, lb=0.525, tr=0.63, dg=0.248, lg=0.248 @@ -322,19 +331,172 @@ subroutine charecteristics_of_fuel ( currentSite ) if ( hlm_masterproc == itrue ) write(fates_log(),*) 'problem with spitfire fuel averaging' endif + ! remove mineral content from net fuel load per Thonicke 2010 + ! for ir calculation in subr. rate_of_spread + ! slevis moved here because rate_of_spread is now called twice/timestep + currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals + currentPatch => currentPatch%younger enddo !end patch loop - end subroutine charecteristics_of_fuel + end subroutine characteristics_of_fuel + + !**************************************************************** + subroutine characteristics_of_crown ( currentSite ) + !****************************************************************. + + !returns the live crown fuel characteristics within each patch. + ! passive_crown_FI is minimum fire intensity to ignite canopy crown fuel + + use SFParamsMod, only : SF_VAL_CWD_FRAC + + type(ed_site_type), intent(in), target :: currentSite + + type(ed_patch_type) , pointer :: currentPatch + type(ed_cohort_type), pointer :: currentCohort + + ! ARGUMENTS + + ! LOCAL + real(r8) :: crown_depth ! depth of crown (m) + real(r8) :: height_cbb ! clear branch bole height or crown base height (m) + real(r8) :: max_height ! max cohort on patch (m) + real(r8) :: crown_ignite_energy ! heat yield for crown (kJ/kg) + real(r8) :: tree_sapw_struct_c ! above-ground tree struct and sap biomass in cohort (kgC) + real(r8) :: leaf_c ! leaf carbon (kgC) + real(r8) :: sapw_c ! sapwood carbon (kgC) + real(r8) :: struct_c ! structure carbon (kgC) + real(r8) :: twig_sapw_struct_c ! above-ground twig sap and struct in cohort (kgC) + real(r8) :: crown_fuel_c ! biomass of 1 hr fuels (leaves,twigs) in cohort (kg C) + real(r8) :: crown_fuel_biomass ! biomass of crown fuel in cohort (kg biomass) + real(r8) :: crown_fuel_per_m ! crown fuel per 1m section in cohort + real(r8) :: height_base_canopy ! lowest height of fuels in patch to carry fire in crown + + integer :: ih ! counter + + real, dimension(70):: biom_matrix ! matrix to track biomass from bottom to 70m + real(r8),parameter :: min_density_canopy_fuel = 0.011_r8 !min canopy fuel density (kg/m3) sufficient to + !propogate fire vertically through canopy + !Scott and Reinhardt 2001 RMRS-RP-29 + real(r8),parameter :: foliar_moist_content = 1.0_r8 !foliar moisture content default 100% Scott & Reinhardt 2001 + + + !returns the live crown fuel characteristics within each patch. + ! passive_crown_FI is the required minimum fire intensity to ignite canopy crown fuel + + currentPatch => currentSite%oldest_patch + + do while(associated(currentPatch)) + !zero Patch level variables + height_base_canopy = 0.0_r8 + currentPatch%canopy_fuel_load = 0.0_r8 + currentPatch%passive_crown_FI = 0.0_r8 + currentPatch%canopy_bulk_density = 0.0_r8 + max_height = 0._r8 + biom_matrix(:) = 0._r8 + + currentCohort=>currentPatch%tallest + do while(associated(currentCohort)) + + !zero cohort level variables + tree_sapw_struct_c = 0.0_r8 + leaf_c = 0.0_r8 + sapw_c = 0.0_r8 + struct_c = 0.0_r8 + twig_sapw_struct_c = 0.0_r8 + crown_fuel_c = 0.0_r8 + crown_fuel_biomass = 0.0_r8 + crown_fuel_per_m = 0.0_r8 + + ! Calculate crown 1hr fuel biomass (leaf, twig sapwood, twig structural biomass) + if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees + + call CrownDepth(currentCohort%hite,currentCohort%pft,crown_depth) + height_cbb = currentCohort%hite - crown_depth + + !find patch max height for stand canopy fuel + if (currentCohort%hite > max_height) then + max_height = currentCohort%hite + endif + + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + + tree_sapw_struct_c = currentCohort%n * & + (prt_params%allom_agb_frac(currentCohort%pft)*(sapw_c + struct_c)) + twig_sapw_struct_c = tree_sapw_struct_c * SF_VAL_CWD_frac(1) !only 1hr fuel + + crown_fuel_c = (currentCohort%n * leaf_c) + twig_sapw_struct_c !crown fuel (kgC) + + crown_fuel_biomass = crown_fuel_c / 0.45_r8 ! crown fuel (kg biomass) + + crown_fuel_per_m = crown_fuel_biomass / crown_depth ! kg biomass per m + + !sort crown fuel into bins from bottom to top of crown + !accumulate across cohorts to find density within canopy 1m sections + do ih = max(1, nint(height_cbb)), max(1, nint(currentCohort%hite)) + biom_matrix(ih) = biom_matrix(ih) + crown_fuel_per_m + end do + + !accumulate available canopy fuel for patch (kg biomass) + ! use this in CFB (crown fraction burn) calculation and FI final + currentPatch%canopy_fuel_load = currentPatch%canopy_fuel_load + crown_fuel_biomass !canopy fuel in patch + + endif !trees only + + currentCohort => currentCohort%shorter; + + enddo !end cohort loop + + biom_matrix(:) = biom_matrix(:) / currentPatch%area !kg biomass/m3 + + !loop from 1m to 70m to find bin with total density = 0.011 kg/m3 + !min canopy fuel density to propogate fire vertically in canopy across patch + do ih=1,70 + if (biom_matrix(ih) > min_density_canopy_fuel) then + height_base_canopy = float(ih) + exit + end if + end do + + !canopy_bulk_denisty (kg/m3) for Patch + if (max_height - height_base_canopy > 0._r8) then + currentPatch%canopy_bulk_density = sum(biom_matrix) / (max_height - height_base_canopy) + else + currentPatch%canopy_bulk_density = 0._r8 + end if + + ! Note: crown_ignition_energy to be calculated based on PFT foliar moisture content from FATES-Hydro + ! or create foliar moisture % based on BTRAN + ! Use foliar_moisture(currentCohort%pft) and compute weighted PFT average with Eq 3 Van Wagner 1977 + ! in place of foliar_moist_content parameter + + ! Eq 3 Van Wagner 1977, Eq 11 Scott & Reinhardt 2001 + ! h = 460.0 + 25.9*m + ! h = crown_ignite_energy (kJ/kg), m = foliar moisture content based on dry fuel (%) + crown_ignite_energy = 460.0 + 25.9 * foliar_moist_content + + ! Crown fuel ignition potential (kW/m), Eq 4 Van Wagner 1977, Eq 11 Scott & Reinhardt 2001 + ! FI = (Czh)**3/2 where z=canopy base height,h=heat of crown ignite energy, FI=fire intensity + ! 0.01 = C, empirical constant Van Wagner 1977 Eq 4 for 6m canopy base height, 100% FMC, FI 2500kW/m + ! passive_crown_FI = min fire intensity to ignite canopy fuel (kW/m or kJ/m/s) + currentPatch%passive_crown_FI = (0.01_r8 * height_base_canopy * crown_ignite_energy)**1.5_r8 + + currentPatch => currentPatch%younger; + + enddo !end patch loop + + end subroutine characteristics_of_crown !***************************************************************** subroutine wind_effect ( currentSite, bc_in) !*****************************************************************. ! Routine called daily from within ED within a site loop. - ! Calculates the effective windspeed based on vegetation charecteristics. + ! Calculates the effective windspeed based on vegetation characteristics. ! currentSite%wind is daily wind converted to m/min for Spitfire units use FatesConstantsMod, only : sec_per_min @@ -414,13 +576,13 @@ subroutine wind_effect ( currentSite, bc_in) end subroutine wind_effect - !***************************************************************** - subroutine rate_of_spread ( currentSite ) + !******************************************************************* + subroutine rate_of_spread ( currentSite) !*****************************************************************. !Routine called daily from within ED within a site loop. !Returns the updated currentPatch%ROS_front value for each patch. - use SFParamsMod, only : SF_val_miner_total, & + use SFParamsMod, only : SF_val_miner_total, & SF_val_part_dens, & SF_val_miner_damp, & SF_val_fuel_energy @@ -431,6 +593,10 @@ subroutine rate_of_spread ( currentSite ) type(ed_patch_type), pointer :: currentPatch + ! ARGUMENTS + + ! LOCAL VARIABLES + ! Rothermal fire spread model parameters. real(r8) beta,beta_op ! weighted average of packing ratio (unitless) real(r8) ir ! reaction intensity (kJ/m2/min) @@ -441,16 +607,24 @@ subroutine rate_of_spread ( currentSite ) real(r8) beta_ratio ! ratio of beta/beta_op real(r8) a_beta ! dummy variable for product of a* beta_ratio for react_v_opt equation real(r8) a,b,c,e ! function of fuel sav + real(r8) time_r ! residence time (min) + real(r8) temp1, temp2 + + real(r8),parameter :: q_dry = 581.0_r8 !heat of pre-ignition of dry fuels (kJ/kg) + real(r8),parameter :: wind_reduce = 0.2_r8 !wind reduction factor (%) - logical, parameter :: debug_windspeed = .false. !for debugging - real(r8),parameter :: q_dry = 581.0_r8 !heat of pre-ignition of dry fuels (kJ/kg) currentPatch=>currentSite%oldest_patch; do while(associated(currentPatch)) + +!! clean up this initialise to zero section?? - ! remove mineral content from net fuel load per Thonicke 2010 for ir calculation - currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals + ! ---initialise parameters to zero.--- + beta_ratio = 0.0_r8; q_ig = 0.0_r8; eps = 0.0_r8; a = 0.0_r8; b = 0.0_r8; c = 0.0_r8; e = 0.0_r8 + phi_wind = 0.0_r8; xi = 0.0_r8; reaction_v_max = 0.0_r8; reaction_v_opt = 0.0_r8; mw_weight = 0.0_r8 + moist_damp = 0.0_r8; ir = 0.0_r8; a_beta = 0.0_r8; + currentPatch%ROS_front = 0.0_r8 ! ----start spreading--- @@ -463,7 +637,7 @@ subroutine rate_of_spread ( currentSite ) ! fraction of fuel array volume occupied by fuel or compactness of fuel bed beta = currentPatch%fuel_bulkd / SF_val_part_dens - ! Equation A6 in Thonicke et al. 2010 + ! Eq A6 in Thonicke et al. 2010 ! packing ratio (unitless) beta_op = 0.200395_r8 *(currentPatch%fuel_sav**(-0.8189_r8)) @@ -476,20 +650,20 @@ subroutine rate_of_spread ( currentSite ) endif ! ---heat of pre-ignition--- - ! Equation A4 in Thonicke et al. 2010 - ! Rothermal EQ12= 250 Btu/lb + 1116 Btu/lb * fuel_eff_moist - ! conversion of Rothermal (1972) EQ12 in BTU/lb to current kJ/kg + ! Eq A4 in Thonicke et al. 2010, Eq 12 Rothermel 1972 + ! 50 Btu/lb + 1116 Btu/lb * fuel_eff_moist + ! conversion of Rothermel (1972) Eq 12 in BTU/lb to current kJ/kg ! q_ig in kJ/kg - q_ig = q_dry +2594.0_r8 * currentPatch%fuel_eff_moist + q_ig = q_dry + 2594.0_r8 * currentPatch%fuel_eff_moist ! ---effective heating number--- - ! Equation A3 in Thonicke et al. 2010. + ! Eq A3 in Thonicke et al. 2010. eps = exp(-4.528_r8 / currentPatch%fuel_sav) - ! Equation A7 in Thonicke et al. 2010 per eqn 49 from Rothermel 1972 + ! Eq A7 in Thonicke et al. 2010 per Eq 49 Rothermel 1972 b = 0.15988_r8 * (currentPatch%fuel_sav**0.54_r8) - ! Equation A8 in Thonicke et al. 2010 per eqn 48 from Rothermel 1972 + ! Eq A8 in Thonicke et al. 2010 per Eq 48 Rothermel 1972 c = 7.47_r8 * (exp(-0.8711_r8 * (currentPatch%fuel_sav**0.55_r8))) - ! Equation A9 in Thonicke et al. 2010. (appears to have typo, using coefficient eqn.50 Rothermel 1972) + ! Eq A9 in Thonicke et al. 2010. (has typo, using coefficient Eq 50 Rothermel 1972) e = 0.715_r8 * (exp(-0.01094_r8 * currentPatch%fuel_sav)) if (debug) then @@ -501,35 +675,35 @@ subroutine rate_of_spread ( currentSite ) if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - e ',e endif - ! Equation A5 in Thonicke et al. 2010 + ! Eq A5 in Thonicke et al. 2010 ! phi_wind (unitless) ! convert current_wspeed (wind at elev relevant to fire) from m/min to ft/min for Rothermel ROS eqn phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) ! ---propagating flux---- - ! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermal 1972 + ! Eq A2 in Thonicke et al.2010 and Eq 42 Rothermel 1972 ! xi (unitless) xi = (exp((0.792_r8 + 3.7597_r8 * (currentPatch%fuel_sav**0.5_r8)) * (beta+0.1_r8))) / & (192_r8+7.9095_r8 * currentPatch%fuel_sav) ! ---reaction intensity---- - ! Equation in table A1 Thonicke et al. 2010. + ! Eq in table A1 Thonicke et al. 2010. a = 8.9033_r8 * (currentPatch%fuel_sav**(-0.7913_r8)) a_beta = exp(a*(1.0_r8-beta_ratio)) !dummy variable for reaction_v_opt equation - ! Equation in table A1 Thonicke et al. 2010. + ! Eq in table A1 Thonicke et al. 2010. ! reaction_v_max and reaction_v_opt = reaction velocity in units of per min - ! reaction_v_max = Equation 36 in Rothermal 1972 and Fig 12 + ! reaction_v_max = Eq 36 in Rothermel 1972 and Fig 12 reaction_v_max = 1.0_r8 / (0.0591_r8 + 2.926_r8* (currentPatch%fuel_sav**(-1.5_r8))) - ! reaction_v_opt = Equation 38 in Rothermal 1972 and Fig 11 + ! reaction_v_opt = Eq 38 in Rothermel 1972 and Fig 11 reaction_v_opt = reaction_v_max*(beta_ratio**a)*a_beta ! mw_weight = relative fuel moisture/fuel moisture of extinction ! average values for litter pools (dead leaves, twigs, small and large branches) plus grass mw_weight = currentPatch%fuel_eff_moist/currentPatch%fuel_mef - ! Equation in table A1 Thonicke et al. 2010. + ! Eq in table A1 Thonicke et al. 2010. ! moist_damp is unitless moist_damp = max(0.0_r8,(1.0_r8 - (2.59_r8 * mw_weight) + (5.11_r8 * (mw_weight**2.0_r8)) - & (3.52_r8*(mw_weight**3.0_r8)))) @@ -540,15 +714,37 @@ subroutine rate_of_spread ( currentSite ) ! write(fates_log(),*) 'ir',gamma_aptr,moist_damp,SF_val_fuel_energy,SF_val_miner_damp - if (((currentPatch%fuel_bulkd) <= 0.0_r8).or.(eps <= 0.0_r8).or.(q_ig <= 0.0_r8)) then + + if (currentPatch%fuel_bulkd <= 0.0_r8 .or. eps <= 0.0_r8 .or. & + q_ig <= 0.0_r8 .or. ir <= 0.0_r8 .or. & + currentPatch%fuel_sav <= 0.0_r8 .or. xi <= 0.0_r8 .or. & + c <= 0.0_r8 .or. beta_ratio <= 0.0_r8 .or. b <= 0.0_r8) then currentPatch%ROS_front = 0.0_r8 - else ! Equation 9. Thonicke et al. 2010. +! currentPatch%ROS_torch = 0.0_r8 ! potentially useful diagnostic + else ! Eq 9. Thonicke et al. 2010. ! forward ROS in m/min currentPatch%ROS_front = (ir*xi*(1.0_r8+phi_wind)) / (currentPatch%fuel_bulkd*eps*q_ig) ! write(fates_log(),*) 'ROS',currentPatch%ROS_front,phi_wind,currentPatch%effect_wspeed ! write(fates_log(),*) 'ros calcs',currentPatch%fuel_bulkd,ir,xi,eps,q_ig + + ! calculate heat release per unit area (HPA)(kJ/m2), Eq 2 Scott & Reinhardt 2001 + ! and residence time (min), Eq 3 Scott & Reinhardt 2001 + time_r = 12.595_r8 / currentPatch%fuel_sav + currentPatch%heat_per_area = ir * time_r + + ! calculate torching index based on wind speed and crown fuels + ! ROS for crown torch initation (m/min), Eq 18 Scott & Reinhardt 2001 + ! temp1 requires max(0... to avoid raising a negative value to a + ! fractional temp2 power. + ! TODO slevis: omitting "- phi_s" after the -1 bc phi_s = 0 for now. + ! @jkshuman recommended naming it slope_factor. +! temp1 = max(0._r8, (60._r8 * currentPatch%passive_crown_FI * & +! currentPatch%fuel_bulkd * eps * q_ig / (currentPatch%heat_per_area * ir * xi) - 1._r8) / & +! (c * beta_ratio**-e)) +! temp2 = 1._r8 / b +! currentPatch%ROS_torch = temp1**temp2 / (54.683_r8 * wind_reduce) endif - ! Equation 10 in Thonicke et al. 2010 + ! Eq 10 in Thonicke et al. 2010 ! backward ROS from Can FBP System (1992) in m/min ! backward ROS wind not changed by vegetation currentPatch%ROS_back = currentPatch%ROS_front*exp(-0.012_r8*currentSite%wind) @@ -583,7 +779,7 @@ subroutine ground_fuel_consumption ( currentSite ) do while(associated(currentPatch)) currentPatch%burnt_frac_litter(:) = 1.0_r8 ! Calculate fraction of litter is burnt for all classes. - ! Equation B1 in Thonicke et al. 2010--- + ! Eq B1 in Thonicke et al. 2010--- do c = 1, nfsc !work out the burnt fraction for all pools, even if those pools dont exist. moist = currentPatch%litter_moisture(c) ! 1. Very dry litter @@ -647,7 +843,7 @@ end subroutine ground_fuel_consumption !***************************************************************** - subroutine area_burnt_intensity ( currentSite, bc_in ) + subroutine area_burnt_intensity ( currentSite, bc_in) !***************************************************************** !returns the updated currentPatch%FI value for each patch. @@ -658,42 +854,44 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) !currentPatch%ROS_front forward ROS (m/min) !currentPatch%TFC_ROS total fuel consumed by flaming front (kgC/m2 of burned area) - use FatesInterfaceTypesMod, only : hlm_spitfire_mode use EDParamsMod, only : ED_val_nignitions use EDParamsMod, only : cg_strikes ! fraction of cloud-to-ground ligtning strikes use FatesConstantsMod, only : years_per_day use SFParamsMod, only : SF_val_fdi_alpha,SF_val_fuel_energy, & - SF_val_max_durat, SF_val_durat_slope, SF_val_fire_threshold + SF_val_max_durat, SF_val_durat_slope, SF_val_fire_threshold type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), pointer :: currentPatch type(bc_in_type), intent(in) :: bc_in - real(r8) ROS !m/s - real(r8) W !kgBiomass/m2 - real(r8) :: tree_fraction_patch ! patch level. no units - real(r8) lb !length to breadth ratio of fire ellipse (unitless) - real(r8) df !distance fire has travelled forward in m - real(r8) db !distance fire has travelled backward in m - real(r8) AB !daily area burnt in m2 per km2 - - real(r8) size_of_fire !in m2 - real(r8) cloud_to_ground_strikes ! [fraction] depends on hlm_spitfire_mode - real(r8) anthro_ign_count ! anthropogenic ignition count/km2/day - integer :: iofp ! index of oldest fates patch + ! ARGUMENTS + + ! LOCAL VARIABLES + real(r8) ROS !rate of spread (m/s) + real(r8) W !available fuel (kgBiomass/m2) + real(r8) :: tree_fraction_patch !patch level. no units + real(r8) df !distance fire has travelled forward (m) + real(r8) db !distance fire has travelled backward (m) + real(r8) AB !daily area burnt (m2 per km2) + real(r8) size_of_fire !in m2 + real(r8) cloud_to_ground_strikes ! [fraction] depends on hlm_spitfire_mode + real(r8) anthro_ign_count ! anthropogenic ignition count/km2/day + integer :: iofp ! index of oldest fates patch real(r8), parameter :: pot_hmn_ign_counts_alpha = 0.0035_r8 ! Potential human ignition counts (alpha in Li et al. 2012) (#/person/month) - real(r8), parameter :: km2_to_m2 = 1000000.0_r8 !area conversion for square km to square m + real(r8), parameter :: km2_to_m2 = 1000000.0_r8 ! area conversion for square km to square m real(r8), parameter :: m_per_min__to__km_per_hour = 0.06_r8 ! convert wind speed from m/min to km/hr - real(r8), parameter :: forest_grassland_lengthtobreadth_threshold = 0.55_r8 ! tree canopy cover below which to use grassland length-to-breadth eqn + real(r8), parameter :: forest_grassland_lengthtobreadth_threshold = 0.55_r8 ! tree canopy cover below which to use + ! grassland length-to-breadth eqn + ! 0.55 = benchmark forest cover, Staver 2010 ! ---initialize site parameters to zero--- currentSite%NF_successful = 0._r8 - ! Equation 7 from Venevsky et al GCB 2002 (modification of equation 8 in Thonicke et al. 2010) + ! Eq 7 from Venevsky et al GCB 2002 (modification of Eqn 8, Thonicke et al. 2010) ! FDI 0.1 = low, 0.3 moderate, 0.75 high, and 1 = extreme ignition potential for alpha 0.000337 if (hlm_spitfire_mode == hlm_sf_successful_ignitions_def) then - currentSite%FDI = 1.0_r8 ! READING "SUCCESSFUL IGNITION" DATA - ! force ignition potential to be extreme + currentSite%FDI = 1.0_r8 ! READING "SUCCESSFUL IGNITION" DATA + ! force ignition potential to be extreme cloud_to_ground_strikes = 1.0_r8 ! cloud_to_ground = 1 = use 100% incoming observed ignitions else ! USING LIGHTNING DATA currentSite%FDI = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI) @@ -732,14 +930,14 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) if (currentSite%NF > 0.0_r8) then - ! Equation 14 in Thonicke et al. 2010 + ! Eq 14 in Thonicke et al. 2010 ! fire duration in minutes currentPatch%FD = (SF_val_max_durat+1.0_r8) / (1.0_r8 + SF_val_max_durat * & exp(SF_val_durat_slope*currentSite%FDI)) if(write_SF == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fire duration minutes',currentPatch%fd endif - !equation 15 in Arora and Boer CTEM model.Average fire is 1 day long. + !Eq 15 in Arora and Boer CTEM model.Average fire is 1 day long. !currentPatch%FD = 60.0_r8 * 24.0_r8 !no minutes in a day tree_fraction_patch = 0.0_r8 @@ -753,15 +951,16 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) endif if ((currentPatch%effect_wspeed*m_per_min__to__km_per_hour) < 1._r8) then !16.67m/min = 1km/hr - lb = 1.0_r8 + currentPatch%lb = 1.0_r8 else - if (tree_fraction_patch > forest_grassland_lengthtobreadth_threshold) then !benchmark forest cover, Staver 2010 - ! EQ 79 forest fuels (Canadian Forest Fire Behavior Prediction System Ont.Inf.Rep. ST-X-3, 1992) - lb = (1.0_r8 + (8.729_r8 * & - ((1.0_r8 -(exp(-0.03_r8 * m_per_min__to__km_per_hour * currentPatch%effect_wspeed)))**2.155_r8))) - else ! EQ 80 grass fuels (CFFBPS Ont.Inf.Rep. ST-X-3, 1992, but with a correction from an errata published within - ! Information Report GLC-X-10 by Bottom et al., 2009 because there is a typo in CFFBPS Ont.Inf.Rep. ST-X-3, 1992) - lb = (1.1_r8*((m_per_min__to__km_per_hour * currentPatch%effect_wspeed)**0.464_r8)) + if (tree_fraction_patch > forest_grassland_lengthtobreadth_threshold) then + ! Eq 79 forest fuels (Canadian Forest Fire Behavior Prediction System Ont.Inf.Rep. ST-X-3, 1992) + currentPatch%lb = (1.0_r8 + (8.729_r8 * & + ((1.0_r8 -(exp(-0.03_r8 * m_per_min__to__km_per_hour * currentPatch%effect_wspeed)))**2.155_r8))) + else ! Eq 80 grass fuels (CFFBPS Ont.Inf.Rep. ST-X-3, 1992, with correction from errata published in + ! Inf.Rep. GLC-X-10 (Bottom et al., 2009) because of typo in CFFBPS Ont.Inf.Rep. ST-X-3, 1992) + currentPatch%lb = 1.1_r8 * & + ((m_per_min__to__km_per_hour * currentPatch%effect_wspeed)**0.464_r8) endif endif @@ -773,15 +972,15 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) df = currentPatch%ROS_front * currentPatch%FD !m ! --- calculate area burnt--- - if(lb > 0.0_r8) then - - ! Equation 1 in Thonicke et al. 2010 + if (currentPatch%lb > 0.0_r8) then + + ! Eq 1 in Thonicke et al. 2010 ! To Do: Connect here with the Li & Levis GDP fire suppression algorithm. - ! Equation 16 in arora and boer model JGR 2005 + ! Eq 16 in arora and boer model JGR 2005 ! AB = AB *3.0_r8 - !size of fire = equation 14 Arora and Boer JGR 2005 (area of an ellipse) - size_of_fire = ((pi_const/(4.0_r8*lb))*((df+db)**2.0_r8)) + !size of fire = Eq 14 Arora and Boer JGR 2005 (area of an ellipse) + size_of_fire = (df + db) * (df + db) * pi_const / (4.0_r8 * currentPatch%lb) ! AB = daily area burnt = size fires in m2 * num ignitions per day per km2 * prob ignition starts fire ! AB = m2 per km2 per day @@ -803,12 +1002,12 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) currentPatch%frac_burnt = 0._r8 endif ! lb - ROS = currentPatch%ROS_front / 60.0_r8 !m/min to m/sec - W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 of burned area to kgbiomass/m2 of burned area + ROS = currentPatch%ROS_front / 60.0_r8 !m/min to m/sec for FI calculation + W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 of burned area to kgbiomass/m2 of burned area - ! EQ 15 Thonicke et al 2010 - !units of fire intensity = (kJ/kg)*(kgBiomass/m2)*(m/min) - currentPatch%FI = SF_val_fuel_energy * W * ROS !kj/m/s, or kW/m + ! Eq 15 Thonicke et al 2010 + !units of fire intensity = (kJ/kg)*(kgBiomass/m2)*(m/sec) + currentPatch%FI = SF_val_fuel_energy * W * ROS !kj/m/s, or kW/m if(write_sf == itrue)then if( hlm_masterproc == itrue ) write(fates_log(),*) 'fire_intensity',currentPatch%fi,W,currentPatch%ROS_front @@ -835,6 +1034,286 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) end subroutine area_burnt_intensity + !***************************************************************** + subroutine active_crown_fire ( currentSite) + !***************************************************************** + + !evaluates if there will be an active crown fire based on canopy fuel and rate of spread + !returns final rate of spread and fire intensity in patch with added fuel from active crown fire. + !currentCohort%fraction_crown_burned is the proportion of crown affected by fire + + use EDParamsMod, only : active_crown_fire_switch + use SFParamsMod, only : SF_val_miner_total, SF_val_part_dens, SF_val_miner_damp, & + SF_val_fuel_energy, SF_val_drying_ratio + + + type(ed_site_type), intent(in), target :: currentSite + + type(ed_patch_type) , pointer :: currentPatch + type(ed_cohort_type), pointer :: currentCohort + + ! ARGUMENTS + + ! LOCAL VARIABLES + ! Active crown Rothermel fire spread model parameters using FM 10 + real(r8) beta,beta_op ! weighted average of packing ratio (unitless) + real(r8) ir ! reaction intensity (kJ/m2/min) + real(r8) xi,eps,phi_wind ! all are unitless + real(r8) q_ig ! heat of pre-ignition (kJ/kg) + real(r8) reaction_v_opt,reaction_v_max !reaction velocity (per min)!optimum and maximum + real(r8) moist_damp,mw_weight ! moisture dampening coefficient and ratio fuel moisture to extinction + real(r8) beta_ratio ! ratio of beta/beta_op + real(r8) a_beta ! dummy variable for product of a* beta_ratio for react_v_opt equation + real(r8) a,b,c,e ! function of fuel sav + real(r8) total_fuel ! total fuel (kg biomass/m2) + real(r8) net_fuel ! net fuel (kg biomass/m2) without minerals + real(r8) fuel_depth ! fuel depth (m) + real(r8) fuel_bd ! fuel bulk density (kg biomass/m3) + real(r8) fuel_sav ! fuels average sav + real(r8) fuel_eff_moist ! fuels effective moisture + real(r8) fuel_moist1hr ! moisture 1 hour fuels + real(r8) fuel_moist10hr ! moisture 10 hour fuels + real(r8) fuel_moist100hr ! moisture 100 hour fuels + real(r8) fuel_moistlive ! moisture live fuels + real(r8) fuel_1hr ! FM 10 1-hr fuel loading (kg/m2) + real(r8) fuel_10hr ! FM 10 10-hr fuel loading (kg/m2) + real(r8) fuel_100hr ! FM 10 100-hr fuel loading (kg/m2) + real(r8) fuel_live ! FM 10 live fuel loading (kg/m2) + real(r8) SAV_1hr ! surface area to volume 1 hour fuels (twigs) + real(r8) SAV_10hr ! surface area to volume 10 hour fuels (small branches) + real(r8) SAV_100hr ! surface area to volume 100 hour fuels (large branches) + real(r8) SAV_live ! surface area to volume live fuels + real(r8) midflame_wind ! 40% of open wind speed, Scott & Reinhardt 2001 + real(r8) db ! distance fire has traveld backward (m) + real(r8) df ! distance fire has travelled forward (m) + real(r8) AB ! daily area burnt (m2 per km2) + real(r8) size_of_fire ! in m2 + real(r8) ROS_active ! actual rate of spread (m/min) using FM 10 fuels + real(r8) ROS_active_min ! minimum rate of spread to ignite active crown fire + real(r8) CI_temp ! temporary variable to calculate wind_active_min + real(r8) wind_active_min ! open windspeed to sustain active crown fire where ROS_SA = ROS_active_min + real(r8) ROS_SA ! rate of spread for surface fire with wind_active_min + real(r8) canopy_frac_burnt ! fraction of canopy fuels consumed (0, surface fire to 1,active crown fire) + real(r8) ROS_final ! final rate of spread for combined surface and canopy spread (m/min) + real(r8) FI_final ! final fireline intensity (kW/m or kJ/m/sec) with canopy consumption + + real(r8),parameter :: q_dry = 581.0_r8 !heat of pre-ignition of dry fuels (kJ/kg) + ! fuel loading, MEF, and depth from Anderson 1982 Aids to determining fuel models for fire behavior + ! SAV values from BEHAVE model Burgan & Rothermel (1984) + real(r8),parameter :: fuel_1hr_ton = 3.01_r8 ! FM 10 1-hr fuel loading (US tons/acre) + real(r8),parameter :: fuel_10hr_ton = 2.0_r8 ! FM 10 10-hr fuel loading (US tons/acre) + real(r8),parameter :: fuel_100hr_ton = 5.01_r8 ! FM 10 100-hr fuel loading (US tons/acre) + real(r8),parameter :: fuel_live_ton = 2.0_r8 ! FM 10 live fuel loading (US tons/acre) + real(r8),parameter :: fuel_mef = 0.25_r8 ! FM 10 moisture of extinction (volumetric) + real(r8),parameter :: fuel_depth_ft= 1.0_r8 ! FM 10 fuel depth (ft) + real(r8),parameter :: sav_1hr_ft = 2000.0_r8 ! FM 10 1-hr SAV (ft2/ft3) + real(r8),parameter :: sav_10hr_ft = 109.0_r8 ! FM 10 10-hr SAV (ft2/ft3) + real(r8),parameter :: sav_100hr_ft = 30.0_r8 ! FM 10 100-hr SAV (ft2/ft3) + real(r8),parameter :: sav_live_ft = 1650.0_r8 ! FM 10 live SAV (ft2/ft3) + real(r8),parameter :: tonnes_acre_to_kg_m2 = 0.2241701 ! convert tons/acre to kg/m2 + real(r8),parameter :: sqft_cubicft_to_sqm_cubicm = 0.03280844 !convert ft2/ft3 to m2/m3 + real(r8),parameter :: canopy_ignite_energy = 18000.0_r8 ! heat yield for canopy fuels (kJ/kg) + real(r8),parameter :: critical_mass_flow_rate = 0.05_r8 ! critical mass flow rate (kg/m2/sec)for crown fire + real(r8),parameter :: km2_to_m2 = 1000000.0_r8 ! area conversion for square km to square m + + integer :: passive_canopy_fuel_flg ! flag if canopy fuel true for vertical spread + + + currentPatch => currentSite%oldest_patch + + do while(associated(currentPatch)) + + if (currentPatch%fire == 1) then + passive_canopy_fuel_flg = 0 !does patch have canopy fuels for vertical spread? + ROS_active = 0.0_r8 + + ! check initiation of passive crown fire + if (currentPatch%FI >= currentPatch%passive_crown_FI) then + passive_canopy_fuel_flg = 1 !enough passive canopy fuels for vertical spread + + ! Calculate rate of spread using FM 10 as in Rothermel 1977 + ! fuel characteristics + fuel_1hr = fuel_1hr_ton * tonnes_acre_to_kg_m2 + fuel_10hr = fuel_10hr_ton * tonnes_acre_to_kg_m2 + fuel_100hr = fuel_100hr_ton * tonnes_acre_to_kg_m2 + fuel_live = fuel_live_ton * tonnes_acre_to_kg_m2 + + total_fuel = fuel_1hr + fuel_10hr + fuel_100hr + fuel_live !total fuel (kg/m2) + + SAV_1hr = sav_1hr_ft * sqft_cubicft_to_sqm_cubicm + SAV_10hr = sav_10hr_ft * sqft_cubicft_to_sqm_cubicm + SAV_100hr = sav_100hr_ft * sqft_cubicft_to_sqm_cubicm + SAV_live = sav_live_ft * sqft_cubicft_to_sqm_cubicm + + fuel_moist1hr = exp(-1.0_r8 * ((SAV_1hr/SF_val_drying_ratio) * currentSite%acc_NI)) + fuel_moist10hr = exp(-1.0_r8 * ((SAV_10hr/SF_val_drying_ratio) * currentSite%acc_NI)) + fuel_moist100hr = exp(-1.0_r8 * ((SAV_100hr/SF_val_drying_ratio) * currentSite%acc_NI)) + fuel_moistlive = exp(-1.0_r8 * ((SAV_live/SF_val_drying_ratio) * currentSite%acc_NI)) + + fuel_depth = fuel_depth_ft *0.3048 !convert to meters + fuel_bd = total_fuel / fuel_depth !fuel bulk density (kg/m3) + + fuel_sav = SAV_1hr *(fuel_1hr/total_fuel) + SAV_10hr*(fuel_10hr/total_fuel) + & + SAV_100hr*(fuel_100hr/total_fuel) + SAV_live*(fuel_live/total_fuel) + + fuel_eff_moist = fuel_moist1hr *(fuel_1hr/total_fuel) + fuel_moist10hr*(fuel_10hr/total_fuel) + & + fuel_moist100hr*(fuel_100hr/total_fuel) + fuel_moistlive*(fuel_live/total_fuel) + + ! remove mineral content from net fuel load + net_fuel = total_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals + + ! ---start spreading--- + !beta = packing ratio (unitless) + beta = fuel_bd / SF_val_part_dens + beta_op = 0.200395_r8 *(fuel_sav**(-0.8189_r8)) + beta_ratio = beta/beta_op + + ! -- heat of pre-ignition -- + q_ig = q_dry + 2594.0_r8 * fuel_eff_moist + + ! ---effective heating number--- + ! Eq A3 in Thonicke et al. 2010. + eps = exp(-4.528_r8 / fuel_sav) + ! Eq A7 in Thonicke et al. 2010 per Eq 49, Rothermel 1972 + b = 0.15988_r8 * (fuel_sav**0.54_r8) + ! Eq A8 in Thonicke et al. 2010 per Eq 48, Rothermel 1972 + c = 7.47_r8 * (exp(-0.8711_r8 * (fuel_sav**0.55_r8))) + ! Eq A9 in Thonicke et al. 2010. (typo in Eq A9, using coefficient Eq 50, Rothermel 1972) + e = 0.715_r8 * (exp(-0.01094_r8 * fuel_sav)) + + midflame_wind = currentSite%wind *0.40_r8 !Scott & Reinhardt 2001 40% open wind speed + + ! Eq A5 in Thonicke et al. 2010 + ! include convert wind from m/min to ft/min for Rothermel ROS eqn + phi_wind = c * ((3.281_r8*midflame_wind)**b)*(beta_ratio**(-e)) !unitless + + ! ---propagating flux = xi (unitless) + ! Eq A2 in Thonicke et al.2010 and Eq 42 Rothermel 1972 + xi = (exp((0.792_r8 + 3.7597_r8 * (fuel_sav**0.5_r8)) * (beta+0.1_r8))) / & + (192_r8+7.9095_r8 * fuel_sav) + + ! ---reaction intensity---- + ! Eq in table A1 Thonicke et al. 2010. + a = 8.9033_r8 * (fuel_sav**(-0.7913_r8)) + a_beta = exp(a*(1.0_r8-beta_ratio)) !dummy variable for reaction_v_opt equation + + ! Eq in table A1 Thonicke et al. 2010. + ! reaction_v_max and reaction_v_opt = reaction velocity in units of per min + ! reaction_v_max = Eq 36 in Rothermel 1972 and Fig 12 + reaction_v_max = 1.0_r8 / (0.0591_r8 + 2.926_r8* (fuel_sav**(-1.5_r8))) + ! reaction_v_opt = Eq 38 in Rothermel 1972 and Fig 11 + reaction_v_opt = reaction_v_max*(beta_ratio**a)*a_beta + + ! mw_weight = relative fuel moisture/fuel moisture of extinction + mw_weight = fuel_eff_moist/fuel_mef + + ! Eq in table A1 Thonicke et al. 2010. (unitless) + moist_damp = max(0.0_r8,(1.0_r8 - (2.59_r8 * mw_weight) + (5.11_r8 * (mw_weight**2.0_r8)) - & + (3.52_r8*(mw_weight**3.0_r8)))) + + ! ir = reaction intenisty in kJ/m2/min + ! sum_fuel as kgBiomass/m2 for ir calculation + ir = reaction_v_opt*(net_fuel)*SF_val_fuel_energy*moist_damp*SF_val_miner_damp + + ! actual ROS (m/min) for FM 10 fuels for open windspeed, Eq 8 Scott & Reinhardt 2001 + ROS_active = 3.34_r8 * ((ir*xi*(1.0_r8+phi_wind)) / (fuel_bd * eps * q_ig)) + + ! critical min rate of spread (m/min) for active crowning + ROS_active_min = (critical_mass_flow_rate / fuel_bd) * 60.0_r8 + + ! check threshold intensity and rate of spread + if (active_crown_fire_switch .and. & + ROS_active >= ROS_active_min) then + currentPatch%active_crown_fire_flg = 1 ! active crown fire ignited + !ROS_final = ROS_surface+CFB(ROS_active - ROS_surface), Eq 21 Scott & Reinhardt 2001 + !with active crown fire CFB (canopy fraction burned) = 100% + canopy_frac_burnt = 1.0_r8 + + else + currentPatch%active_crown_fire_flg = 0 ! only passive crown fire with partial crown burnt + + ! phi_slope is not used yet. consider adding with later + ! development + ! calculate open wind speed critical to sustain active crown + ! fire Eq 20 Scott & Reinhardt + if (ir > 0._r8 .and. currentPatch%canopy_bulk_density > 0._r8) then + CI_temp = ((164.8_r8 * eps * q_ig)/(ir * currentPatch%canopy_bulk_density)) - 1.0_r8 + else + CI_temp = 0._r8 + end if + + wind_active_min = 0.0457_r8 * (CI_temp / 0.001612_r8)**0.7_r8 + + ! use open wind speed "wind_active_min" for ROS surface fire + ! where ROS_SA=ROS_active_min + ROS_SA = (ir * xi * (1.0_r8 + wind_active_min)) / (fuel_bd * eps * q_ig) + + ! canopy fraction burnt, Eq 28 Scott & Reinhardt Appendix A + canopy_frac_burnt = max(0._r8, min(1.0_r8, & + (currentPatch%ROS_front - ROS_active_min) / (ROS_SA - ROS_active_min))) + + endif !check intensity & ROS for active crown fire thresholds + + !ROS_final = ROS_surface+CFB(ROS_active - ROS_surface), Eq 21 Scott & Reinhardt 2001 + ROS_final = currentPatch%ROS_front + & + canopy_frac_burnt * (ROS_active - currentPatch%ROS_front) + + ! recalculate area burned with new ROS_front value from ROS_final + ! ---- re-calculate length of major axis for df using new ROS_front value from ROS final--- + db = currentPatch%ROS_back * currentPatch%FD !(m) + df = ROS_final * currentPatch%FD !(m) update with ROS final + + ! update ROS_front with ROS_final for output variable + ! if changing, expect only an increase in ROS_front + currentPatch%ROS_front = max(ROS_final, currentPatch%ROS_front) + + ! --- calculate updated area burnt using df from ROS final--- + if (currentPatch%lb > 0.0_r8) then + + ! Eq 1 in Thonicke et al. 2010 + ! To Do: Connect here with the Li & Levis GDP fire suppression algorithm. + ! Eq 16 in arora and boer model JGR 2005 + ! AB = AB *3.0_r8 + + !size of fire = Eq 14 Arora and Boer JGR 2005 (area of an ellipse) + size_of_fire = (df + db) * (df + db) * pi_const / (4.0_r8 * currentPatch%lb) + + ! AB = daily area burnt = size fires in m2 * num ignitions per day per km2 * prob ignition starts fire + ! AB = m2 per km2 per day + ! the denominator in the units of currentSite%NF is total gridcell area, but since we assume that ignitions + ! are equally probable across patches, currentSite%NF is equivalently per area of a given patch + ! thus AB has units of m2 burned area per km2 patch area per day + AB = size_of_fire * currentSite%NF * currentSite%FDI + + ! frac_burnt + ! just a unit conversion from AB, to become area burned per area patch per day, + ! or just the fraction of the patch burned on that day + currentPatch%frac_burnt = (min(0.99_r8, AB / km2_to_m2)) + + if(write_SF == itrue)then + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt + endif + + else + currentPatch%frac_burnt = 0.0_r8 + endif ! lb + + !final fireline intensity (kJ/m/sec or kW/m), Eq 22 Scott & Reinhardt 2001 + FI_final = (currentPatch%heat_per_area + & + currentPatch%canopy_fuel_load * canopy_ignite_energy * canopy_frac_burnt) * & + currentPatch%ROS_front / 60._r8 + ! update patch FI to adjust according to potential canopy fuel consumed (passive and active) + ! if changing, expect only an increase in FI + currentPatch%FI = max(FI_final, currentPatch%FI) + + endif !check if passive crown fire? + endif !fire? + + currentPatch => currentPatch%younger; + + enddo !end patch loop + + end subroutine active_crown_fire !***************************************************************** @@ -898,6 +1377,7 @@ subroutine crown_scorching ( currentSite ) end subroutine crown_scorching + !***************************************************************** subroutine crown_damage ( currentSite ) !***************************************************************** @@ -909,42 +1389,40 @@ subroutine crown_damage ( currentSite ) type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - real(r8) :: crown_depth ! Depth of crown in meters - + + real(r8) :: crown_depth ! depth of crown (m) + real(r8) :: height_cbb ! clear branch bole height or crown base height (m) for cohort + currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) + do while(associated(currentPatch)) + !zero Patch level variables + if (currentPatch%fire == 1) then currentCohort=>currentPatch%tallest do while(associated(currentCohort)) currentCohort%fraction_crown_burned = 0.0_r8 - if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees only - ! Flames lower than bottom of canopy. - ! c%hite is height of cohort + if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees + ! height_cbb = clear branch bole height at base of crown (m) + ! inst%crown = crown_depth_frac (PFT) call CrownDepth(currentCohort%hite,currentCohort%pft,crown_depth) - - if (currentPatch%Scorch_ht(currentCohort%pft) < & - (currentCohort%hite-crown_depth)) then - currentCohort%fraction_crown_burned = 0.0_r8 - else - ! Flames part of way up canopy. - ! Equation 17 in Thonicke et al. 2010. - ! flames over bottom of canopy but not over top. - if ((currentCohort%hite > 0.0_r8).and.(currentPatch%Scorch_ht(currentCohort%pft) >= & - (currentCohort%hite-crown_depth))) then - - currentCohort%fraction_crown_burned = (currentPatch%Scorch_ht(currentCohort%pft) - & - (currentCohort%hite - crown_depth))/crown_depth - - else - ! Flames over top of canopy. - currentCohort%fraction_crown_burned = 1.0_r8 - endif - - endif + height_cbb = currentCohort%hite - crown_depth + + ! Equation 17 in Thonicke et al. 2010 + ! flames over bottom of canopy, and potentially over top of + ! canopy + if (currentCohort%hite > 0.0_r8 .and. & + currentPatch%Scorch_ht(currentCohort%pft) >= height_cbb) then + if (currentPatch%active_crown_fire_flg == 0) then + currentCohort%fraction_crown_burned = min(1.0_r8, & + ((currentPatch%Scorch_ht(currentCohort%pft) - height_cbb) / crown_depth)) + else ! active crown fire occurring + currentCohort%fraction_crown_burned = 1.0_r8 + end if + endif !SH frac crown burnt calculation ! Check for strange values. currentCohort%fraction_crown_burned = min(1.0_r8, max(0.0_r8,currentCohort%fraction_crown_burned)) endif !trees only @@ -954,6 +1432,7 @@ subroutine crown_damage ( currentSite ) currentCohort => currentCohort%shorter; enddo !end cohort loop + endif !fire? currentPatch => currentPatch%younger; @@ -983,7 +1462,8 @@ subroutine cambial_damage_kill ( currentSite ) if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; - do while(associated(currentCohort)) + do while(associated(currentCohort)) + currentCohort%cambial_mort = 0.0_r8 if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees only ! Equation 21 in Thonicke et al 2010 bt = EDPftvarcon_inst%bark_scaler(currentCohort%pft)*currentCohort%dbh ! bark thickness. diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 06bcd1858a..edc4c0cf7c 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -642,6 +642,12 @@ subroutine init_patches( nsites, sites, bc_in) currentPatch%scorch_ht(:) = 0._r8 currentPatch%frac_burnt = 0._r8 currentPatch%burnt_frac_litter(:) = 0._r8 + currentPatch%canopy_bulk_density = 0._r8 + currentPatch%canopy_fuel_load = 0._r8 + currentPatch%passive_crown_FI = 0._r8 + currentPatch%ros_torch = 0._r8 + currentPatch%heat_per_area = 0._r8 + currentPatch%lb = 0._r8 currentPatch => currentPatch%older enddo diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 51926d275a..9c494fbdf3 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -61,10 +61,9 @@ module EDParamsMod real(r8),protected, public :: ED_val_patch_fusion_tol real(r8),protected, public :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry integer,protected, public :: stomatal_model !switch for choosing between stomatal conductance models, 1 for Ball-Berry, 2 for Medlyn - - logical,protected, public :: active_crown_fire ! flag, 1=active crown fire 0=no active crown fire - character(len=param_string_length),parameter :: fates_name_active_crown_fire = "fates_fire_active_crown_fire" + logical,protected, public :: active_crown_fire_switch ! flag, 1=active crown fire 0=no active crown fire + character(len=param_string_length),parameter :: fates_name_active_crown_fire = "fates_fire_active_crown_fire" real(r8), protected, public :: cg_strikes ! fraction of cloud to ground lightning strikes (0-1) character(len=param_string_length),parameter :: fates_name_cg_strikes="fates_fire_cg_strikes" @@ -454,7 +453,7 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=fates_name_active_crown_fire, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - + call fates_params%RegisterParameter(name=fates_name_cg_strikes, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -623,7 +622,7 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=fates_name_active_crown_fire, & data=tmpreal) - active_crown_fire = (abs(tmpreal-1.0_r8)