diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 6f022ccfbd..e84658ef86 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -2500,6 +2500,7 @@ 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%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/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 8a38366217..d8ee395f62 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -196,6 +196,8 @@ module FatesPatchMod real(r8) :: fuel_eff_moist ! effective avearage fuel moisture content of the ground fuel ! (incl. live grasses. omits 1000hr fuels) real(r8) :: litter_moisture(nfsc) ! moisture of litter [m3/m3] + real(r8) :: canopy_bulk_density + ! fire spread real(r8) :: ros_front ! rate of forward spread of fire [m/min] @@ -204,6 +206,7 @@ module FatesPatchMod real(r8) :: tau_l ! duration of lethal heating [min] real(r8) :: fi ! average fire intensity of flaming front [kJ/m/s] or [kW/m] integer :: fire ! is there a fire? [1=yes; 0=no] + integer :: active_crown_fire_flg real(r8) :: fd ! fire duration [min] ! fire effects @@ -382,9 +385,11 @@ subroutine NanValues(this) this%tau_l = nan this%fi = nan this%fire = fates_unset_int + this%active_crown_fire_flg = fates_unset_int this%fd = nan this%scorch_ht(:) = nan this%frac_burnt = nan + this%canopy_bulk_density = nan this%tfc_ros = nan this%burnt_frac_litter(:) = nan @@ -459,7 +464,8 @@ subroutine ZeroValues(this) this%fi = 0.0_r8 this%fd = 0.0_r8 this%scorch_ht(:) = 0.0_r8 - this%frac_burnt = 0.0_r8 + this%frac_burnt = 0.0_r8 + this%canopy_bulk_density = 0.0_r8 this%tfc_ros = 0.0_r8 this%burnt_frac_litter(:) = 0.0_r8 @@ -732,6 +738,7 @@ subroutine Dump(this) write(fates_log(),*) 'pa%c_stomata = ',this%c_stomata write(fates_log(),*) 'pa%c_lblayer = ',this%c_lblayer write(fates_log(),*) 'pa%disturbance_rates = ',this%disturbance_rates(:) + write(fates_log(),*) 'pa%active_crown_fire_flg = ', this%active_crown_fire_flg write(fates_log(),*) 'pa%anthro_disturbance_label = ',this%anthro_disturbance_label write(fates_log(),*) '----------------------------------------' diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index ef245b04f9..b056fcc6b8 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -36,7 +36,6 @@ module SFMainMod use FatesLitterMod , only : TR_SF use FatesLitterMod , only : litter_type - use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : fnrt_organ @@ -53,11 +52,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 +88,21 @@ subroutine fire_model( currentSite, bc_in) type (fates_patch_type), pointer :: currentPatch + real(r8) :: canopy_fuel_load ! available canopy fuel load in patch (kg biomass) + real(r8) :: passive_crown_FI ! fire intensity for ignition of passive canopy fuel (kW/m) + real(r8) :: ROS_torch ! ROS for crown torch initation (m/min) + real(r8) :: lb !length to breadth ratio of fire ellipse (unitless) + real(r8) :: heat_per_area ! heat release per unit area (kJ/m2) for surface fuel + + !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 +113,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 rate_of_spread(currentSite) + call characteristics_of_fuel(currentSite) + call characteristics_of_crown(currentSite, canopy_fuel_load, passive_crown_FI) + call rate_of_spread(currentSite, passive_crown_FI, ROS_torch, heat_per_area) call ground_fuel_consumption(currentSite) - call area_burnt_intensity(currentSite, bc_in) + call area_burnt_intensity(currentSite, bc_in, lb) + call active_crown_fire (currentSite,passive_crown_FI,canopy_fuel_load,ROS_torch,heat_per_area,lb) call crown_scorching(currentSite) call crown_damage(currentSite) call cambial_damage_kill(currentSite) @@ -163,21 +176,21 @@ 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 @@ -186,73 +199,75 @@ 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 currentPatch => currentSite%oldest_patch; - do while(associated(currentPatch)) + + do while (associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then + if (currentPatch%nocomp_pft_label .ne. nocomp_bareground) then - litt_c => currentPatch%litter(element_pos(carbon12_element)) - - ! How much live grass is there? - currentPatch%livegrass = 0.0_r8 - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - ! for grasses sum all aboveground tissues - if( prt_params%woody(currentCohort%pft) == ifalse)then - - currentPatch%livegrass = currentPatch%livegrass + & - ( currentCohort%prt%GetState(leaf_organ, carbon12_element) + & - currentCohort%prt%GetState(sapw_organ, carbon12_element) + & - currentCohort%prt%GetState(struct_organ, carbon12_element) ) * & - currentCohort%n/currentPatch%area + litt_c => currentPatch%litter(element_pos(carbon12_element)) + + ! How much live grass is there? + currentPatch%livegrass = 0.0_r8 + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + ! for grasses sum all aboveground tissues + if( prt_params%woody(currentCohort%pft) == ifalse)then + + currentPatch%livegrass = currentPatch%livegrass + & + ( currentCohort%prt%GetState(leaf_organ, carbon12_element) + & + currentCohort%prt%GetState(sapw_organ, carbon12_element) + & + currentCohort%prt%GetState(struct_organ, carbon12_element) ) * & + currentCohort%n/currentPatch%area endif currentCohort => currentCohort%shorter - enddo - - ! There are SIX fuel classes - ! 1:4) four CWD_AG pools (twig, s branch, l branch, trunk), 5) dead leaves and 6) live grass - ! NCWD =4 NFSC = 6 - ! tw_sf = 1, lb_sf = 3, tr_sf = 4, dl_sf = 5, lg_sf = 6, - + enddo - if(write_sf == itrue)then + ! There are SIX fuel classes + ! 1:4) four CWD_AG pools (twig, s branch, l branch, trunk), 5) dead leaves and 6) live grass + ! NCWD =4 NFSC = 6 + ! tw_sf = 1, lb_sf = 3, tr_sf = 4, dl_sf = 5, lg_sf = 6, + + + if(write_sf == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) ' leaf_litter1 ',sum(litt_c%leaf_fines(:)) if ( hlm_masterproc == itrue ) write(fates_log(),*) ' leaf_litter2 ',sum(litt_c%ag_cwd(:)) if ( hlm_masterproc == itrue ) write(fates_log(),*) ' leaf_litter3 ',currentPatch%livegrass - endif + endif - currentPatch%sum_fuel = sum(litt_c%leaf_fines(:)) + & - sum(litt_c%ag_cwd(:)) + & - currentPatch%livegrass - if(write_SF == itrue)then + currentPatch%sum_fuel = sum(litt_c%leaf_fines(:)) + & + sum(litt_c%ag_cwd(:)) + & + currentPatch%livegrass + if (write_SF == itrue) then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sum fuel', currentPatch%sum_fuel,currentPatch%area - endif - ! =============================================== - ! Average moisture, bulk density, surface area-volume and moisture extinction of fuel - ! ================================================ - - if (currentPatch%sum_fuel > 0.0) then + endif + + ! =============================================== + ! Average moisture, bulk density, surface area-volume and moisture extinction of fuel + ! ================================================ + + if (currentPatch%sum_fuel > 0.0) then ! Fraction of fuel in litter classes currentPatch%fuel_frac(dl_sf) = sum(litt_c%leaf_fines(:))/ currentPatch%sum_fuel currentPatch%fuel_frac(tw_sf:tr_sf) = litt_c%ag_cwd(:) / currentPatch%sum_fuel - if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff2a ', & - lg_sf,currentPatch%livegrass,currentPatch%sum_fuel + if (write_sf == itrue) then + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff2a ', & + lg_sf,currentPatch%livegrass,currentPatch%sum_fuel endif 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 @@ -264,28 +279,28 @@ subroutine charecteristics_of_fuel ( currentSite ) ! dead leaves and twigs included in 1hr pool per Thonicke (2010) ! Calculate fuel moisture for trunks to hold value for fuel consumption alpha_FMC(tw_sf:dl_sf) = SF_val_SAV(tw_sf:dl_sf)/SF_val_drying_ratio - + fuel_moisture(tw_sf:dl_sf) = exp(-1.0_r8 * alpha_FMC(tw_sf:dl_sf) * currentSite%acc_NI) - + if(write_SF == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fm ',fuel_moisture - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%acc_NI - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fm ',fuel_moisture + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%acc_NI + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC endif - + ! live grass moisture is a function of SAV and changes via Nesterov Index ! along the same relationship as the 1 hour fuels (live grass has same SAV as dead grass, ! but retains more moisture with this calculation.) fuel_moisture(lg_sf) = exp(-1.0_r8 * ((SF_val_SAV(tw_sf)/SF_val_drying_ratio) * currentSite%acc_NI)) - + ! Average properties over the first three litter pools (twigs, s branches, l branches) currentPatch%fuel_bulkd = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * SF_val_FBD(tw_sf:lb_sf)) currentPatch%fuel_sav = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * SF_val_SAV(tw_sf:lb_sf)) currentPatch%fuel_mef = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * MEF(tw_sf:lb_sf)) currentPatch%fuel_eff_moist = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * fuel_moisture(tw_sf:lb_sf)) - if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff4 ',currentPatch%fuel_eff_moist + if (write_sf == itrue) then + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff4 ',currentPatch%fuel_eff_moist endif ! Add on properties of dead leaves and live grass pools (5 & 6) currentPatch%fuel_bulkd = currentPatch%fuel_bulkd + sum(currentPatch%fuel_frac(dl_sf:lg_sf) * SF_val_FBD(dl_sf:lg_sf)) @@ -299,28 +314,25 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_sav = currentPatch%fuel_sav * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) currentPatch%fuel_mef = currentPatch%fuel_mef * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) currentPatch%fuel_eff_moist = currentPatch%fuel_eff_moist * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) - + ! Pass litter moisture into the fuel burning routine (all fuels: twigs,s branch,l branch,trunk,dead leaves,live grass) ! (wo/me term in Thonicke et al. 2010) currentPatch%litter_moisture(tw_sf:lb_sf) = fuel_moisture(tw_sf:lb_sf)/MEF(tw_sf:lb_sf) currentPatch%litter_moisture(tr_sf) = fuel_moisture(tr_sf)/MEF(tr_sf) currentPatch%litter_moisture(dl_sf) = fuel_moisture(dl_sf)/MEF(dl_sf) currentPatch%litter_moisture(lg_sf) = fuel_moisture(lg_sf)/MEF(lg_sf) - - else - - if(write_SF == itrue)then - - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'no litter fuel at all',currentPatch%patchno, & - currentPatch%sum_fuel,sum(litt_c%ag_cwd(:)),sum(litt_c%leaf_fines(:)) + else + if (write_SF == itrue) then + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'no litter fuel at all',currentPatch%patchno, & + currentPatch%sum_fuel,sum(litt_c%ag_cwd(:)),sum(litt_c%leaf_fines(:)) endif currentPatch%fuel_sav = sum(SF_val_SAV(1:nfsc))/(nfsc) ! make average sav to avoid crashing code. - if ( hlm_masterproc == itrue .and. write_SF == itrue)then - write(fates_log(),*) 'problem with spitfire fuel averaging' + if ( hlm_masterproc == itrue .and. write_SF == itrue) then + write(fates_log(),*) 'problem with spitfire fuel averaging' end if - + ! FIX(SPM,032414) refactor...should not have 0 fuel unless everything is burnt ! off. currentPatch%fuel_eff_moist = 0.0000000001_r8 @@ -329,27 +341,184 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_mef = 0.0000000001_r8 currentPatch%sum_fuel = 0.0000000001_r8 - endif - ! check values. - ! FIX(SPM,032414) refactor... - if(write_SF == itrue.and.currentPatch%fuel_sav <= 0.0_r8.or.currentPatch%fuel_bulkd <= & - 0.0_r8.or.currentPatch%fuel_mef <= 0.0_r8.or.currentPatch%fuel_eff_moist <= 0.0_r8)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'problem with spitfire fuel averaging' - endif - endif !nocomp_pft_label check - currentPatch => currentPatch%younger + endif + ! check values. + ! FIX(SPM,032414) refactor... + if (write_SF == itrue.and.currentPatch%fuel_sav <= 0.0_r8.or.currentPatch%fuel_bulkd <= & + 0.0_r8.or.currentPatch%fuel_mef <= 0.0_r8.or.currentPatch%fuel_eff_moist <= 0.0_r8) then + 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 + + end if enddo !end patch loop - end subroutine charecteristics_of_fuel + end subroutine characteristics_of_fuel + + !**************************************************************** + subroutine characteristics_of_crown ( currentSite, canopy_fuel_load, passive_crown_FI) + !****************************************************************. + + !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(fates_patch_type) , pointer :: currentPatch + type(fates_cohort_type), pointer :: currentCohort + + ! ARGUMENTS + real(r8), intent(out) :: canopy_fuel_load ! available canopy fuel load in patch (kg biomass) + real(r8), intent(out) :: passive_crown_FI ! min fire intensity to ignite canopy fuel (kW/m) + + ! 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 + + !! check to see if active_crown_fire is enabled? + + do while(associated(currentPatch)) + !zero Patch level variables + height_base_canopy = 0.0_r8 + canopy_fuel_load = 0.0_r8 + passive_crown_FI = 0.0_r8 + currentPatch%canopy_bulk_density = 0.0_r8 + +! if (currentPatch%active_crown_fire == 1) then + + 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%height,currentCohort%pft,crown_depth) + height_cbb = currentCohort%height - crown_depth + + !find patch max height for stand canopy fuel + if (currentCohort%height > max_height) then + max_height = currentCohort%height + endif + + leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) + sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element) + struct_c = currentCohort%prt%GetState(struct_organ, carbon12_element) + + 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 = int(height_cbb), int(currentCohort%height) + 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 + canopy_fuel_load = 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_density (kg/m3) for Patch + currentPatch%canopy_bulk_density = sum(biom_matrix) / (max_height - height_base_canopy) + + ! 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) + passive_crown_FI = (0.01_r8 * height_base_canopy * crown_ignite_energy)**1.5_r8 + +! endif !active crown fire? + + 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 @@ -445,13 +614,13 @@ subroutine wind_effect ( currentSite, bc_in) end subroutine wind_effect - !***************************************************************** - subroutine rate_of_spread ( currentSite ) + !******************************************************************* + subroutine rate_of_spread ( currentSite, ROS_torch, passive_crown_FI, heat_per_area) !*****************************************************************. !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 @@ -460,6 +629,13 @@ subroutine rate_of_spread ( currentSite ) type(fates_patch_type), pointer :: currentPatch + ! ARGUMENTS + real(r8), intent(out) :: ROS_torch ! ROS for crown torch initation (m/min) + real(r8), intent(out) :: heat_per_area ! heat release per unit area (kJ/m2) for surface fuel + real(r8), intent(in) :: passive_crown_FI ! min fire intensity to ignite canopy fuel (kW/m or kJ/m/s) + + ! 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) @@ -470,9 +646,11 @@ 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),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; @@ -480,8 +658,11 @@ subroutine rate_of_spread ( currentSite ) if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - ! 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--- @@ -494,7 +675,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)) @@ -507,20 +688,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 @@ -532,35 +713,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)))) @@ -571,15 +752,28 @@ 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 currentPatch%ROS_front = 0.0_r8 - else ! Equation 9. Thonicke et al. 2010. + ROS_torch = 0.0_r8 + 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 / currentPatch%fuel_sav + 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 + ROS_torch = (1.0 / 54.683 * wind_reduce)* & + ((((60.0*passive_crown_FI*currentPatch%fuel_bulkd*eps*q_ig)/heat_per_area*ir*xi)-1.0) & + / (c*beta_ratio)**(-1*e))**1/b 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) @@ -618,7 +812,7 @@ subroutine ground_fuel_consumption ( currentSite ) 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 @@ -684,7 +878,7 @@ end subroutine ground_fuel_consumption !***************************************************************** - subroutine area_burnt_intensity ( currentSite, bc_in ) + subroutine area_burnt_intensity ( currentSite, bc_in, lb) !***************************************************************** !returns the updated currentPatch%FI value for each patch. @@ -695,42 +889,46 @@ 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 FatesInterfaceTypesMod, only : hlm_spitfire_mode ! TODO slevis: redundant? 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(fates_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 + real(r8), intent(out) :: lb !length to breadth ratio of fire ellipse (unitless) + + ! 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) @@ -780,14 +978,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 @@ -803,12 +1001,12 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) if ((currentPatch%effect_wspeed*m_per_min__to__km_per_hour) < 1._r8) then !16.67m/min = 1km/hr 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) + 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) 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 Wotton et al., 2009 for a typo in CFFBPS Ont.Inf.Rep. ST-X-3, 1992) + 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) lb = (1.1_r8*((m_per_min__to__km_per_hour * currentPatch%effect_wspeed)**0.464_r8)) endif endif @@ -822,13 +1020,13 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) ! --- calculate area burnt--- if(lb > 0.0_r8) then - - ! Equation 1 in Thonicke et al. 2010 + + ! 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 = Eq 14 Arora and Boer JGR 2005 (area of an ellipse) size_of_fire = ((pi_const/(4.0_r8*lb))*((df+db)**2.0_r8)) ! AB = daily area burnt = size fires in m2 * num ignitions per day per km2 * prob ignition starts fire @@ -851,12 +1049,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 @@ -884,6 +1082,282 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) end subroutine area_burnt_intensity + !***************************************************************** + subroutine active_crown_fire ( currentSite, canopy_fuel_load, ROS_torch, & + lb, heat_per_area, passive_crown_FI) + !***************************************************************** + + !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 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(fates_patch_type) , pointer :: currentPatch + type(fates_cohort_type), pointer :: currentCohort + + ! ARGUMENTS + real(r8), intent(in) :: ROS_torch ! ROS for crown torch initation (m/min) + real(r8), intent(in) :: canopy_fuel_load ! available canopy fuel load in patch (kg biomass) + real(r8), intent(in) :: lb !length to breadth ratio of fire ellipse (unitless) + real(r8), intent(in) :: heat_per_area ! heat release per unit area (kJ/m2) for surface fuel + real(r8), intent(in) :: passive_crown_FI ! fire intensity for ignition of passive canopy fuel (kW/m) + + ! 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 + real(r8) fuel_10hr + real(r8) fuel_100hr + real(r8) fuel_live + 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_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 + + !! check to see if active_crown_fire is enabled + + 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 >= 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 (currentPatch%FI >= passive_crown_FI .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 + ROS_final = currentPatch%ROS_front + canopy_frac_burnt*(ROS_active-currentPatch%ROS_front) + + 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 + CI_temp = ((164.8_r8 * eps * q_ig)/(ir * currentPatch%canopy_bulk_density)) - 1.0_r8 + + 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 = (min(1.0_r8, ((currentPatch%ROS_front - ROS_active_min) & + /(ROS_SA - ROS_active_min)))) + + !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) + + endif !check intensity & ROS for active crown fire thresholds + + ! 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 + currentPatch%ROS_front = ROS_final + + ! --- calculate updated area burnt using df from ROS final--- + if(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 = ((pi_const/(4.0_r8*lb))*((df+db)**2.0_r8)) + + ! 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 = ((heat_per_area + (canopy_fuel_load*canopy_ignite_energy*canopy_frac_burnt))& + *currentPatch%ROS_front)/60.0 + ! update patch FI to adjust according to potential canopy fuel consumed (passive and active) + currentPatch%FI = FI_final + + endif !check if passive crown fire? + endif !fire? + + currentPatch => currentPatch%younger; + + enddo !end patch loop + + end subroutine active_crown_fire !***************************************************************** @@ -950,6 +1424,7 @@ subroutine crown_scorching ( currentSite ) end subroutine crown_scorching + !***************************************************************** subroutine crown_damage ( currentSite ) !***************************************************************** @@ -961,12 +1436,14 @@ subroutine crown_damage ( currentSite ) type(fates_patch_type) , pointer :: currentPatch type(fates_cohort_type), pointer :: currentCohort - real(r8) :: crown_depth ! Depth of crown in meters - - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) + 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)) + !zero Patch level variables if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then if (currentPatch%fire == 1) then @@ -974,31 +1451,25 @@ subroutine crown_damage ( currentSite ) do while(associated(currentCohort)) currentCohort%fraction_crown_burned = 0.0_r8 - if ( prt_params%woody(currentCohort%pft) == itrue) then !trees only - ! Flames lower than bottom of canopy. - ! c%height 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%height,currentCohort%pft,crown_depth) - - if (currentPatch%Scorch_ht(currentCohort%pft) < & - (currentCohort%height-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%height > 0.0_r8).and.(currentPatch%Scorch_ht(currentCohort%pft) >= & - (currentCohort%height-crown_depth))) then - - currentCohort%fraction_crown_burned = (currentPatch%Scorch_ht(currentCohort%pft) - & - (currentCohort%height - crown_depth))/crown_depth - - else - ! Flames over top of canopy. - currentCohort%fraction_crown_burned = 1.0_r8 - endif - - endif + height_cbb = currentCohort%height - crown_depth + + ! Equation 17 in Thonicke et al. 2010 + ! flames over bottom of canopy, and potentially over top of + ! canopy + if (currentCohort%height > 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 @@ -1008,6 +1479,7 @@ subroutine crown_damage ( currentSite ) currentCohort => currentCohort%shorter; enddo !end cohort loop + endif !fire? endif !nocomp_pft_label check @@ -1040,8 +1512,9 @@ subroutine cambial_damage_kill ( currentSite ) if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; - do while(associated(currentCohort)) - if ( prt_params%woody(currentCohort%pft) == itrue) then !trees only + 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. ! Equation 20 in Thonicke et al. 2010. diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 438d387213..413b2fadb7 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -605,9 +605,6 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_history_height_bin_edges, dimension_shape=dimension_shape_1d, & dimension_names=dim_names_height) - 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) @@ -815,10 +812,6 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=name_dev_arbitrary, & data=dev_arbitrary) - call fates_params%RetrieveParameter(name=fates_name_active_crown_fire, & - data=tmpreal) - active_crown_fire = (abs(tmpreal-1.0_r8)