diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 9efacbf416..3d8e6ddf9b 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -534,11 +534,12 @@ subroutine nan_cohort(cc_p) currentCohort%ddbhdt = nan ! time derivative of dbh ! FIRE - currentCohort%fraction_crown_burned = nan ! proportion of crown affected by fire - currentCohort%passive_crown_fire_flg = nan ! flag to indicate passive crown fire ignition - currentCohort%cambial_mort = nan ! probability that trees dies due to cambial char P&R (1986) - currentCohort%crownfire_mort = nan ! probability of tree post-fire mortality due to crown scorch - currentCohort%fire_mort = nan ! post-fire mortality from cambial and crown damage assuming two are independent + currentCohort%fraction_crown_burned = nan ! proportion of crown affected by fire + currentCohort%passive_crown_fire_flg = 9999 ! flag to indicate passive crown fire ignition + currentCohort%active_crown_fire_flg = 9999 ! flag to indicate active crown fire ignition + currentCohort%cambial_mort = nan ! probability that trees dies due to cambial char P&R (1986) + currentCohort%crownfire_mort = nan ! probability of tree post-fire mortality due to crown scorch + currentCohort%fire_mort = nan ! post-fire mortality from cambial and crown damage assuming two are independent end subroutine nan_cohort @@ -580,7 +581,8 @@ subroutine zero_cohort(cc_p) currentcohort%year_net_uptake(:) = 999._r8 ! this needs to be 999, or trimming of new cohorts will break. currentcohort%ts_net_uptake(:) = 0._r8 - currentCohort%size_class = 1 + currentcohort%fraction_crown_burned = 0._r8 + currentCohort%size_class = 1 currentCohort%seed_prod = 0._r8 currentCohort%size_class_lasttimestep = 0 currentcohort%npp_acc_hold = 0._r8 diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index a47e93f903..a77a84d001 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -31,7 +31,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 @@ -114,7 +113,7 @@ subroutine fire_danger_index ( currentSite, bc_in) !***************************************************************** ! currentSite%acc_NI is the accumulated Nesterov fire danger index - use SFParamsMod, only : SF_val_fdi_a, SF_val_fdi_b + use SFParamsMod, only : SF_val_fdi_a, SF_val_fdi_b use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm use FatesConstantsMod , only : sec_per_day @@ -167,7 +166,7 @@ subroutine charecteristics_of_fuel ( currentSite ) type(ed_patch_type), pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - type(litter_type), pointer :: litt_c + 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. @@ -439,8 +438,8 @@ subroutine rate_of_spread ( currentSite ) 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 - logical,parameter :: debug_windspeed = .false. !for debugging - real(r8),parameter :: q_dry = 581.0_r8 !heat of pre-ignition of dry fuels (kJ/kg) + 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; @@ -573,8 +572,8 @@ subroutine ground_fuel_consumption ( currentSite ) !returns the the hypothetic fuel consumed by the fire use SFParamsMod, only : SF_val_miner_total, SF_val_min_moisture, & - SF_val_mid_moisture, SF_val_low_moisture_Coeff, SF_val_low_moisture_Slope, & - SF_val_mid_moisture_Coeff, SF_val_mid_moisture_Slope + SF_val_mid_moisture, SF_val_low_moisture_Coeff, SF_val_low_moisture_Slope, & + SF_val_mid_moisture_Coeff, SF_val_mid_moisture_Slope type(ed_site_type) , intent(in), target :: currentSite type(ed_patch_type), pointer :: currentPatch @@ -663,8 +662,8 @@ subroutine fire_intensity ( currentSite ) !currentPatch%TFC_ROS total fuel consumed by flaming front (kgC/m2) use FatesInterfaceMod, only : hlm_use_spitfire - use SFParamsMod, only : SF_val_fdi_alpha,SF_val_fuel_energy, & - SF_val_max_durat, SF_val_durat_slope + use SFParamsMod, only : SF_val_fdi_alpha,SF_val_fuel_energy, & + SF_val_max_durat, SF_val_durat_slope type(ed_site_type), intent(inout), target :: currentSite @@ -688,11 +687,11 @@ subroutine fire_intensity ( currentSite ) !'decide_fire' subroutine shortened and put in here... if (currentPatch%FI >= fire_threshold) then ! 50kW/m is the threshold for a self-sustaining fire currentPatch%fire = 1 ! Fire... :D - + ! Equation 7 from Venevsky et al GCB 2002 (modification of equation 8 in Thonicke et al. 2010) ! FDI 0.1 = low, 0.3 moderate, 0.75 high, and 1 = extreme ignition potential for alpha 0.000337 currentSite%FDI = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI) - + ! Equation 14 in Thonicke et al. 2010 ! fire duration in minutes @@ -706,7 +705,7 @@ subroutine fire_intensity ( currentSite ) !currentPatch%FD = 60.0_r8 * 24.0_r8 !no minutes in a day else currentPatch%fire = 0 ! No fire... :-/ - currentPatch%FD = 0.0_r8 + currentPatch%FD = 0.0_r8 endif ! FIX(SPM,032414) needs a refactor ! FIX(RF,032414) : should happen outside of SF loop - doing all spitfire code is inefficient otherwise. @@ -737,7 +736,6 @@ subroutine area_burnt ( currentSite ) real(r8) AB !daily area burnt in m2 per km2 real(r8) size_of_fire !in m2 real(r8),parameter :: km2_to_m2 = 1000000.0_r8 !area conversion for square km to square m - integer g, p ! ---initialize site parameters to zero--- currentSite%frac_burnt = 0.0_r8 @@ -776,7 +774,7 @@ subroutine area_burnt ( currentSite ) ! --- calculate area burnt--- if(lb > 0.0_r8) then - + ! Equation 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 @@ -815,7 +813,7 @@ subroutine crown_scorching ( currentSite ) type(ed_site_type), intent(in), target :: currentSite - type(ed_patch_type), pointer :: currentPatch + type(ed_patch_type), pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort real(r8) :: f_ag_bmass ! fraction of tree cohort's above-ground biomass as a proportion of total patch ag tree biomass @@ -889,7 +887,7 @@ subroutine crown_damage ( currentSite ) !returns the updated currentCohort%fraction_crown_burned for each tree cohort within each patch. !currentCohort%fraction_crown_burned is the proportion of crown affected by fire - + type(ed_site_type), intent(in), target :: currentSite @@ -897,6 +895,12 @@ subroutine crown_damage ( currentSite ) type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort + real(r8) :: leaf_c ! leaf carbon [kg] + real(r8) :: crown_fuel_bulkd ! leaf biomass per unit area divided by canopy depth [kg/m3] Van Wagner used this method according to Scott (2006) p. 12 + real(r8), parameter :: low_heat_of_combustion = 12700.0_r8 ! [kJ/kg] + real(r8), parameter :: critical_mass_flow_rate = 0.05_r8 ! [kg/m2/s] value for conifer forests; if available for other vegetation, move to the params file? + real(r8) active_crown_FI ! critical fire intensity for active crown fire ignition (kW/m) + real(r8) ignite_active_crown ! ratio for ignition of active crown fire,EQ 14b Bessie & Johnson 1995 real(r8) :: crown_depth ! depth of crown (m) real(r8) :: height_cbb ! clear branch bole height or crown base height (m) real(r8) :: passive_crown_FI ! critical fire intensity for passive crown fire ignition (kW/m) @@ -912,11 +916,13 @@ subroutine crown_damage ( currentSite ) do while(associated(currentCohort)) currentCohort%fraction_crown_burned = 0.0_r8 - if (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then !trees only ! height_cbb = clear branch bole height at base of crown (m) ! inst%crown = crown_depth_frac (PFT) + active_crown_FI = 0.0_r8 !critical fire intensity for active crown fire + ignite_active_crown = 0.0_r8 !ratio for ignition of active crown fire,EQ 14b Bessie & Johnson 1995 + currentCohort%active_crown_fire_flg = 0 !flag for active crown fire ignition crown_depth = currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft) height_cbb = currentCohort%hite - crown_depth passive_crown_FI = 0.0_r8 @@ -925,7 +931,7 @@ subroutine crown_damage ( currentSite ) ! Evaluate for passive crown fire ignition - if (EDPftvarcon_inst%crown_fire(currentCohort%pft) > 0.0_r8) then + if (EDPftvarcon_inst%crown_fire(currentCohort%pft) > 0.001_r8) then ! Note: crown_ignition_energy to be calculated based on foliar moisture content from FATES-Hydro ! EQ 3 Van Wagner 1977 @@ -944,22 +950,56 @@ subroutine crown_damage ( currentSite ) if (ignite_crown >= 1.0_r8) then currentCohort%passive_crown_fire_flg = 1 ! passive crown fire ignited - ! "...in passive crown fires and high intensity surface fires trees can survive. Jack, red - ! and white pine survive...scars used in dating fires" + write(fates_log(),*) 'SF currentCohort%passive_crown_fire_flg = ', currentCohort%passive_crown_fire_flg ! slevis diag + ! "...in passive crown fires and high intensity surface + ! fires trees can survive. Jack, red + ! and white pine survive...scars used in dating fires" ! Johnson, E.A. 1992 Fire and Veg Dynamics: North American boreal forest. Cambridge Press - currentCohort%fraction_crown_burned = EDPftvarcon_inst%crown_fire(currentCohort%pft) - ! else ! evaluate crown damage based on scorch height + currentCohort%fraction_crown_burned = EDPftvarcon_inst%crown_fire(currentCohort%pft) + + ! Critical intensity for active crowning (kW/m) + ! EQ 12 Bessie and Johnson 1995 + ! Fuels / 0.45 to get biomass but note that the 0.45 + ! cancels out and could be removed. Also dividing + ! critical_mass_flow_rate by 3.34, an empirical + ! constant in Bessie & Johnson 1995 + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + crown_fuel_bulkd = currentCohort%n * leaf_c / & + (0.45_r8 * currentCohort%c_area * crown_depth) + active_crown_FI = critical_mass_flow_rate * & + low_heat_of_combustion * currentPatch%sum_fuel / & + (0.45_r8 * 3.34_r8 * crown_fuel_bulkd) + + ! Initiate active crown fire? + ! EQ 14b Bessie & Johnson 1995 + ignite_active_crown = currentPatch%FI / active_crown_FI + + if (ignite_active_crown >= 1.0_r8) then + currentCohort%active_crown_fire_flg = 1 ! active crown fire ignited + write(fates_log(),*) 'SF currentCohort%active_crown_fire_flg = ', currentCohort%active_crown_fire_flg ! slevis diag + currentCohort%fraction_crown_burned = 1.0_r8 + ! TODO Additional effects of active crown fire. + ! Currently in code design phase; see + ! https://github.com/NGEET/fates/issues/573 + else + write(fates_log(),*) 'SF currentPatch%FI < active_crown_FI = ', currentPatch%FI, active_crown_FI ! slevis diag + endif ! ignite active crown fire + else ! crown damage based on scorch height done below + write(fates_log(),*) 'SF currentPatch%FI < passive_crown_FI = ', currentPatch%FI, passive_crown_FI ! slevis diag endif ! ignite passive crown fire - ! else no crown fire today - endif ! crown fire intensity - ! else ! not crown fire plant - endif ! evaluate passive crown fire + else ! crown fire not possible + write(fates_log(),*) 'SF currentPatch%FI < crown_fire_threshold = ', currentPatch%FI, crown_fire_threshold ! slevis diag + endif ! fire intensity vs. crown fire threshold + else ! not a crown-fire pft + write(fates_log(),*) 'SF Not a crown-fire pft' ! slevis diag + endif ! crown-fire pft ! For surface fires, are flames in the canopy? ! height_cbb is clear branch bole height or height of bottom of canopy ! Equation 17 in Thonicke et al. 2010 if (currentCohort%hite > 0.0_r8 .and. currentPatch%SH > height_cbb & - .and. currentCohort%passive_crown_fire_flg == 0) then + .and. currentCohort%passive_crown_fire_flg == 0 & + .and. currentCohort%active_crown_fire_flg == 0) then currentCohort%fraction_crown_burned = min(1.0_r8, ((currentPatch%SH - height_cbb)/crown_depth)) endif !SH frac crown burnt calculation @@ -1002,7 +1042,7 @@ subroutine cambial_damage_kill ( currentSite ) if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; do while(associated(currentCohort)) - !currentCohort%cambial_mort = 0.0_r8 !testing this removal + currentCohort%cambial_mort = 0.0_r8 if (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then !trees only ! Equation 21 in Thonicke et al 2010 bt = EDPftvarcon_inst%bark_scaler(currentCohort%pft)*currentCohort%dbh ! bark thickness. @@ -1059,7 +1099,7 @@ subroutine post_fire_mortality ( currentSite ) ! Equation 22 in Thonicke et al. 2010. currentCohort%crownfire_mort = EDPftvarcon_inst%crown_kill(currentCohort%pft)*currentCohort%fraction_crown_burned**3.0_r8 ! Equation 18 in Thonicke et al. 2010. - currentCohort%fire_mort = max(0.0_r8,min(1.0_r8,currentCohort%crownfire_mort+currentCohort%cambial_mort- & + currentCohort%fire_mort = max(0._r8,min(1.0_r8,currentCohort%crownfire_mort+currentCohort%cambial_mort- & (currentCohort%crownfire_mort*currentCohort%cambial_mort))) !joint prob. else currentCohort%fire_mort = 0.0_r8 !Set to zero. Grass mode of death is removal of leaves. diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 131931b5c2..d35db2a910 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -344,7 +344,8 @@ module EDTypesMod ! FIRE real(r8) :: fraction_crown_burned ! proportion of crown affected by fire:- - real(r8) :: passive_crown_fire_flg ! flag for passive crown fire ignition (1=ignition) + integer :: active_crown_fire_flg ! flag for active crown fire ignition + integer :: passive_crown_fire_flg ! flag for passive crown fire ignition (1=ignition) real(r8) :: cambial_mort ! probability that trees dies due to cambial char ! (conditional on the tree being subjected to the fire) real(r8) :: crownfire_mort ! probability of tree post-fire mortality @@ -1011,6 +1012,7 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%ddbhdt = ', ccohort%ddbhdt write(fates_log(),*) 'co%dbdeaddt = ', ccohort%dbdeaddt write(fates_log(),*) 'co%fraction_crown_burned = ', ccohort%fraction_crown_burned + write(fates_log(),*) 'co%active_crown_fire_flg = ', ccohort%active_crown_fire_flg write(fates_log(),*) 'co%passive_crown_fire_flg = ', ccohort%passive_crown_fire_flg write(fates_log(),*) 'co%fire_mort = ', ccohort%fire_mort write(fates_log(),*) 'co%crownfire_mort = ', ccohort%crownfire_mort diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index b1870beee0..44767ed085 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -4540,12 +4540,12 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m5_si_scpf ) call this%set_history_var(vname='CROWNFIREMORT_SCPF', units = 'N/ha/yr', & - long='crown fire mortality by pft/size',use_default='inactive', & + long='crown fire mortality by pft/size',use_default='active', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_crownfiremort_si_scpf ) call this%set_history_var(vname='CAMBIALFIREMORT_SCPF', units = 'N/ha/yr', & - long='cambial fire mortality by pft/size',use_default='inactive', & + long='cambial fire mortality by pft/size',use_default='active', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_cambialfiremort_si_scpf )