diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 943c5abf37..a12406c907 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -534,11 +534,14 @@ subroutine nan_cohort(cc_p) currentCohort%ddbhdt = nan ! time derivative of dbh ! FIRE + + currentCohort%sh = nan ! scorch height affecting crown (m) currentCohort%fraction_crown_burned = nan ! proportion of crown affected by fire 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 !-------------------------------------------------------------------------------------! @@ -579,6 +582,7 @@ 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%sh = 0._r8 currentcohort%fraction_crown_burned = 0._r8 currentCohort%size_class = 1 currentCohort%seed_prod = 0._r8 @@ -1652,7 +1656,8 @@ subroutine copy_cohort( currentCohort,copyc ) n%dhdt = o%dhdt n%ddbhdt = o%ddbhdt - ! FIRE + ! FIRE + n%sh = o%sh n%fraction_crown_burned = o%fraction_crown_burned n%fire_mort = o%fire_mort n%crownfire_mort = o%crownfire_mort diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 1f6501aa44..f3aee3ca67 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1980,7 +1980,6 @@ subroutine zero_patch(cp_p) currentPatch%fire = 999 ! sr decide_fire.1=fire hot enough to proceed. 0=stop everything- no fires today currentPatch%fd = 0.0_r8 ! fire duration (mins) currentPatch%ros_back = 0.0_r8 ! backward ros (m/min) - currentPatch%sh = 0.0_r8 ! average scorch height for the patch(m) currentPatch%frac_burnt = 0.0_r8 ! fraction burnt daily currentPatch%burnt_frac_litter(:) = 0.0_r8 currentPatch%btran_ft(:) = 0.0_r8 @@ -2296,7 +2295,6 @@ subroutine fuse_2_patches(csite, dp, rp) rp%fi = (dp%fi*dp%area + rp%fi*rp%area) * inv_sum_area rp%fd = (dp%fd*dp%area + rp%fd*rp%area) * inv_sum_area rp%ros_back = (dp%ros_back*dp%area + rp%ros_back*rp%area) * inv_sum_area - rp%sh = (dp%sh*dp%area + rp%sh*rp%area) * inv_sum_area rp%frac_burnt = (dp%frac_burnt*dp%area + rp%frac_burnt*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 diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index ca1077e5bd..17401a0f85 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -22,6 +22,7 @@ module SFMainMod use EDtypesMod , only : AREA use EDtypesMod , only : DL_SF use EDtypesMod , only : FIRE_THRESHOLD + use EDtypesMod , only : crown_fire_threshold use EDTypesMod , only : TW_SF use EDtypesMod , only : LB_SF use EDtypesMod , only : LG_SF @@ -50,9 +51,8 @@ module SFMainMod public :: charecteristics_of_fuel public :: rate_of_spread public :: ground_fuel_consumption - public :: fire_intensity public :: wind_effect - public :: area_burnt + public :: area_burnt_intensity public :: crown_scorching public :: crown_damage public :: cambial_damage_kill @@ -97,8 +97,7 @@ subroutine fire_model( currentSite, bc_in) call charecteristics_of_fuel(currentSite) call rate_of_spread(currentSite) call ground_fuel_consumption(currentSite) - call fire_intensity(currentSite) - call area_burnt(currentSite) + call area_burnt_intensity(currentSite) call crown_scorching(currentSite) call crown_damage(currentSite) call cambial_damage_kill(currentSite) @@ -113,7 +112,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 @@ -166,7 +165,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. @@ -438,7 +437,9 @@ 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 + 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; @@ -481,7 +482,9 @@ subroutine rate_of_spread ( currentSite ) ! Rothermal EQ12= 250 Btu/lb + 1116 Btu/lb * fuel_eff_moist ! conversion of Rothermal (1972) EQ12 in BTU/lb to current kJ/kg ! q_ig in kJ/kg - q_ig = 581.0_r8 +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. @@ -571,8 +574,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 @@ -650,108 +653,76 @@ subroutine ground_fuel_consumption ( currentSite ) end subroutine ground_fuel_consumption + !***************************************************************** - subroutine fire_intensity ( currentSite ) - !***************************************************************** + subroutine area_burnt_intensity ( currentSite ) + !***************************************************************** + !returns the updated currentPatch%FI value for each patch. !currentPatch%FI avg fire intensity of flaming front during day. Backward ROS plays no role here. kJ/m/s or kW/m. !currentSite%FDI probability that an ignition will start a fire + !currentSite%NF number of lighting strikes per day per km2 !currentPatch%ROS_front forward ROS (m/min) !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, & + use EDParamsMod, only : ED_val_nignitions + 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 type(ed_site_type), intent(inout), target :: currentSite - type(ed_patch_type), pointer :: currentPatch real(r8) ROS !m/s real(r8) W !kgBiomass/m2 - - currentPatch => currentSite%oldest_patch; - - do while(associated(currentPatch)) - ROS = currentPatch%ROS_front / 60.0_r8 !m/min to m/sec - W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 to kgbiomass/m2 - - !units of fire intensity = (kJ/kg)*(kgBiomass/m2)*(m/min) - 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 - endif - !'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 - - 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. - !currentPatch%FD = 60.0_r8 * 24.0_r8 !no minutes in a day - else - currentPatch%fire = 0 ! No fire... :-/ - 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. - if( hlm_use_spitfire == ifalse )then - currentPatch%fire = 0 !fudge to turn fire off - endif - - currentPatch => currentPatch%younger; - enddo !end patch loop - - end subroutine fire_intensity - - - !***************************************************************** - subroutine area_burnt ( currentSite ) - !***************************************************************** - - use EDParamsMod, only : ED_val_nignitions - use FatesConstantsMod, only : years_per_day - - type(ed_site_type), intent(inout), target :: currentSite - type(ed_patch_type), pointer :: currentPatch - 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) NF !number of lightning strikes per day per km2 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 + real(r8),parameter :: km2_to_m2 = 1000000.0_r8 !area conversion for square km to square m + real(r8),parameter :: CG_strikes = 0.20_r8 !cloud to ground lightning strikes + !Latham and Williams (2001) ! ---initialize site parameters to zero--- - currentSite%frac_burnt = 0.0_r8 + currentSite%frac_burnt = 0.0_r8 - !NF = number of lighting strikes per day per km2 - NF = ED_val_nignitions * years_per_day + + ! 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) + + !NF = number of lighting strikes per day per km2 scaled by cloud to ground strikes + currentSite%NF = ED_val_nignitions * years_per_day * CG_strikes ! If there are 15 lightning strikes per year, per km2. (approx from NASA product for S.A.) - ! then there are 15 * 1/365 strikes/km2 each day. + ! then there are 15 * 1/365 strikes/km2 each day + currentPatch => currentSite%oldest_patch; do while(associated(currentPatch)) ! ---initialize patch parameters to zero--- + currentPatch%fire = 0 + currentPatch%FD = 0.0_r8 currentPatch%frac_burnt = 0.0_r8 + - if (currentPatch%fire == 1) then + if (currentSite%NF > 0.0_r8) then + + ! Equation 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. + !currentPatch%FD = 60.0_r8 * 24.0_r8 !no minutes in a day + + ! The feedback between vegetation structure and ellipse size if turned off for now, ! to reduce the positive feedback in the syste, ! This will also be investigated by William Hoffmans proposal. @@ -785,38 +756,64 @@ subroutine area_burnt ( currentSite ) !AB = daily area burnt = size fires in m2 * num ignitions per day per km2 * prob ignition starts fire !AB = m2 per km2 per day - AB = size_of_fire * NF * currentSite%FDI + AB = size_of_fire * currentSite%NF * currentSite%FDI - !frac_burnt in units of m2 here. - currentPatch%frac_burnt = min(0.99_r8, AB / km2_to_m2) + !frac_burnt + currentPatch%frac_burnt = (min(0.99_r8, AB / km2_to_m2)) * currentPatch%area/area if(write_SF == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt endif - endif - endif! fire - ! convert frac_burnt to % prior to accumulating at site level - currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt * currentPatch%area/area + endif ! lb + + ROS = currentPatch%ROS_front / 60.0_r8 !m/min to m/sec + W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 to kgbiomass/m2 + + ! EQ 15 Thonicke et al 2010 + !units of fire intensity = (kJ/kg)*(kgBiomass/m2)*(m/min)*unitless_fraction + currentPatch%FI = SF_val_fuel_energy * W * ROS * currentPatch%frac_burnt !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 + endif + + !'decide_fire' subroutine + if (currentPatch%FI > fire_threshold) then !track fires greater than kW/m2 energy threshold + currentPatch%fire = 1 ! Fire... :D + + else + currentPatch%fire = 0 ! No fire... :-/ + currentPatch%FD = 0.0_r8 + currentPatch%frac_burnt = 0.0_r8 + endif + + endif! NF ignitions check + + + ! accumulate frac_burnt % at site level + currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt currentPatch => currentPatch%younger enddo !end patch loop - end subroutine area_burnt + end subroutine area_burnt_intensity + + !***************************************************************** subroutine crown_scorching ( currentSite ) !***************************************************************** - !currentPatch%SH !average scorch height for the patch(m) - !currentPatch%FI average fire intensity of flaming front during day. kW/m. + + !currentPatch%FI average fire intensity of flaming front during day. kW/m. + !currentCohort%SH scorch height for the cohort(m) 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 real(r8) :: tree_ag_biomass ! total amount of above-ground tree biomass in patch. kgC/m2 real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] @@ -827,7 +824,6 @@ subroutine crown_scorching ( currentSite ) do while(associated(currentPatch)) tree_ag_biomass = 0.0_r8 - f_ag_bmass = 0.0_r8 if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; do while(associated(currentCohort)) @@ -840,36 +836,18 @@ subroutine crown_scorching ( currentSite ) tree_ag_biomass = tree_ag_biomass + & currentCohort%n * (leaf_c + & EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)*(sapw_c + struct_c)) - endif !trees only - - currentCohort=>currentCohort%shorter; - - enddo !end cohort loop - - !This loop weights the scorch height for the contribution of each cohort to the overall biomass. - - ! does this do anything? I think it might be redundant? RF. - currentPatch%SH = 0.0_r8 - currentCohort => currentPatch%tallest; - do while(associated(currentCohort)) - if (EDPftvarcon_inst%woody(currentCohort%pft) == 1 & - .and. (tree_ag_biomass > 0.0_r8)) then !trees only + - 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) - - f_ag_bmass = currentCohort%n * (leaf_c + & - EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)*(sapw_c + struct_c)) & - / tree_ag_biomass + currentCohort%SH = 0.0_r8 + if (tree_ag_biomass > 0.0_r8) then - !equation 16 in Thonicke et al. 2010 - if(write_SF == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'currentPatch%SH',currentPatch%SH,f_ag_bmass - endif - !2/3 Byram (1959) - currentPatch%SH = currentPatch%SH + f_ag_bmass * & - EDPftvarcon_inst%fire_alpha_SH(currentCohort%pft) * (currentPatch%FI**0.667_r8) + !Equation 16 in Thonicke et al. 2010 !Van Wagner 1973 EQ8 based on 2/3 Byram (1959) + currentCohort%SH = EDPftvarcon_inst%fire_alpha_SH(currentCohort%pft) * (currentPatch%FI**0.667_r8) + + if(write_SF == itrue)then + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'currentCohort%SH',currentCohort%SH + endif + endif ! tree biomass endif !trees only currentCohort=>currentCohort%shorter; @@ -887,43 +865,169 @@ 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 + + + 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 + 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) :: passive_crown_FI ! fire intensity for ignition from passive canopy fuel (kW/m), EQ 8 + real(r8) :: ignite_passive_crown ! ratio for ignition from passive canopy fuel,EQ 14 Bessie & Johnson 1995 + 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) :: canopy_bulk_density ! density of canopy fuel on patch + real(r8) :: height_base_canopy ! lowest height of fuels to carry fire in crown + integer :: ih ! counter + integer :: passive_canopy_fuel_flg ! flag if canopy fuel true for vertical spread + + 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 :: crown_ignite_energy = 3060_r8 !crown ignition energy (kJ/kg) Van Wagner 1977 + + currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) + do while(associated(currentPatch)) + + !zero Patch level variables + passive_crown_FI = 0.0_r8 + ignite_passive_crown = 0.0_r8 + biom_matrix = 0.0_r8 + canopy_bulk_density = 0.0_r8 + max_height = 0.0_r8 + height_base_canopy = 0.0_r8 + if (currentPatch%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 (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then !trees only + + crown_depth = currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft) + 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 * & + (EDPftvarcon_inst%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%hite) + biom_matrix(ih) = biom_matrix(ih) + crown_fuel_per_m + end do + + 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 + 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 EQ3 Van Wagner 1977 + ! in place of crown_ignite_energy parameter + + ! EQ 3 Van Wagner 1977 + ! h = crown_ignite_energy (kJ/kg), m = foliar moisture content based on dry fuel (%) + ! crown_ignite_energy = 460 + 26 * m + + ! Crown fuel ignition potential, EQ 8 Bessie and Johnson 1995, EQ 4 Van Wagner 1977 + ! FI = (Czh)**3/2 where z=canopy base height,h=heat of crown ignite energy, FI=fire intensity + ! 0.01 = C from Van Wagner 1977 EQ4 for canopy base height 6m, 100% FMC, and FI 2500kW/m + passive_crown_FI = (0.01_r8 * height_base_canopy * crown_ignite_energy)**1.5_r8 + + passive_canopy_fuel_flg = 0 !does patch have canopy fuels for vertical spread? + + ! Initiation of passive crown fire, EQ 14a Bessie and Johnson 1995 + ! Are the canopy fuels in the stand large enough to support vertical spread of fire? + ignite_passive_crown = currentPatch%FI/passive_crown_FI + + if (ignite_passive_crown >= 1.0_r8) then + passive_canopy_fuel_flg = 1 !enough passive canopy fuels for vertical spread + endif + + !evaluate active crown fire conditions + + currentCohort=>currentPatch%tallest + do while(associated(currentCohort)) currentCohort%fraction_crown_burned = 0.0_r8 + if (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then !trees only - ! Flames lower than bottom of canopy. - ! c%hite is height of cohort - if (currentPatch%SH < (currentCohort%hite-currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft))) 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%SH >= & - (currentCohort%hite-currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft)))) then - - currentCohort%fraction_crown_burned = (currentPatch%SH-currentCohort%hite*(1.0_r8- & - EDPftvarcon_inst%crown(currentCohort%pft)))/(currentCohort%hite* & - EDPftvarcon_inst%crown(currentCohort%pft)) - - else - ! Flames over top of canopy. - currentCohort%fraction_crown_burned = 1.0_r8 - endif + + ! height_cbb = clear branch bole height at base of crown (m) + ! inst%crown = crown_depth_frac (PFT) + crown_depth = currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft) + height_cbb = currentCohort%hite - crown_depth - endif + + ! 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. currentCohort%SH >= height_cbb) then + !add active flag !.and. currentCohort%active_crown_fire_flg == 0) then + + currentCohort%fraction_crown_burned = min(1.0_r8, ((currentCohort%SH - height_cbb)/crown_depth)) + !add active flag !.and. currentCohort%active_crown_fire_flg == 1) then + !currentCohort%fraction_crown_burned = 1.0_r8 !1 for active fire = true + 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 @@ -962,7 +1066,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 !testing this removal 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. @@ -1019,7 +1124,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._r8,min(1.0_r8,currentCohort%crownfire_mort+currentCohort%cambial_mort- & + currentCohort%fire_mort = max(0.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/EDInitMod.F90 b/main/EDInitMod.F90 index 9613d79957..7e89c71ade 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -174,6 +174,7 @@ subroutine zero_site( site_in ) ! FIRE site_in%acc_ni = 0.0_r8 ! daily nesterov index accumulating over time. time unlimited theoretically. + site_in%NF = 0.0_r8 ! daily lightning strikes per km2 site_in%frac_burnt = 0.0_r8 ! burn area read in from external file do el=1,num_elements @@ -280,6 +281,7 @@ subroutine set_site_properties( nsites, sites ) sites(s)%dstatus = dstat sites(s)%acc_NI = acc_NI + sites(s)%NF = 0.0_r8 sites(s)%frac_burnt = 0.0_r8 diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 0e4a9f7720..807747a252 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -48,9 +48,10 @@ module EDPftvarcon real(r8), allocatable :: crown(:) ! fraction of the height of the plant ! that is occupied by crown. For fire model. real(r8), allocatable :: bark_scaler(:) ! scaler from dbh to bark thickness. For fire model. - real(r8), allocatable :: crown_kill(:) ! scaler on fire death. For fire model. + real(r8), allocatable :: crown_kill(:) ! crown resistance to fire. (1 = none) For fire model. + real(r8), allocatable :: active_crown_fire(:) ! Is plant susceptible to active crown fire? real(r8), allocatable :: initd(:) ! initial seedling density - real(r8), allocatable :: seed_suppl(:) ! seeds that come from outside the gridbox. + real(r8), allocatable :: seed_suppl(:) ! seeds that come from outside the gridbox. real(r8), allocatable :: BB_slope(:) ! ball berry slope parameter real(r8), allocatable :: seed_alloc_mature(:) ! fraction of carbon balance allocated to @@ -381,6 +382,10 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_fire_active_crown_fire' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_recruit_initd' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -813,6 +818,10 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%crown_kill) + name = 'fates_fire_active_crown_fire' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%active_crown_fire) + name = 'fates_recruit_initd' call fates_params%RetreiveParameterAllocate(name=name, & data=this%initd) @@ -1722,6 +1731,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'crown = ',EDPftvarcon_inst%crown write(fates_log(),fmt0) 'bark_scaler = ',EDPftvarcon_inst%bark_scaler write(fates_log(),fmt0) 'crown_kill = ',EDPftvarcon_inst%crown_kill + write(fates_log(),fmt0) 'active_crown_fire = ',EDPftvarcon_inst%active_crown_fire write(fates_log(),fmt0) 'initd = ',EDPftvarcon_inst%initd write(fates_log(),fmt0) 'seed_suppl = ',EDPftvarcon_inst%seed_suppl write(fates_log(),fmt0) 'BB_slope = ',EDPftvarcon_inst%BB_slope diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 6ed7598d63..29b201db21 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -146,6 +146,8 @@ module EDTypesMod integer, parameter, public :: lg_sf = 6 ! array index of live grass pool for spitfire real(r8), parameter, public :: fire_threshold = 50.0_r8 ! threshold for fires that spread or go out. KWm-2 (Pyne 1986) + real(r8), parameter, public :: crown_fire_threshold = 200.0_r8 ! threshold for passive crown fire ignition. KWm-2 (Bessie & Johnson 1995) + ! PATCH FUSION real(r8), parameter, public :: force_patchfuse_min_biomass = 0.005_r8 ! min biomass (kg / m2 patch area) below which to force-fuse patches @@ -342,7 +344,9 @@ module EDTypesMod real(r8) :: dbdeaddt ! time derivative of dead biomass : KgC/year ! FIRE + real(r8) :: sh ! scorch height: m 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) 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 @@ -519,7 +523,6 @@ module EDTypesMod 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 real(r8) :: fd ! fire duration: mins - real(r8) :: sh ! average scorch height: m ! FIRE EFFECTS real(r8) :: frac_burnt ! fraction burnt: frac gridcell/day @@ -687,7 +690,8 @@ module EDTypesMod real(r8) :: wind ! daily wind in m/min for Spitfire units real(r8) :: acc_ni ! daily nesterov index accumulating over time. real(r8) :: fdi ! daily probability an ignition event will start a fire - real(r8) :: frac_burnt ! fraction of soil burnt in this day. + real(r8) :: NF ! daily ignitions in km2 + real(r8) :: frac_burnt ! fraction of area burnt in this day. ! PLANT HYDRAULICS type(ed_site_hydr_type), pointer :: si_hydr @@ -1009,6 +1013,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%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 write(fates_log(),*) 'co%cambial_mort = ', ccohort%cambial_mort diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index b1870beee0..57d30aa4c3 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -141,7 +141,7 @@ module FatesHistoryInterfaceMod integer :: ih_TFC_ROS_pa integer :: ih_fire_intensity_pa integer :: ih_fire_area_pa - integer :: ih_scorch_height_pa + ! integer :: ih_scorch_height_pa integer :: ih_fire_fuel_bulkd_pa integer :: ih_fire_fuel_eff_moist_pa integer :: ih_fire_fuel_sav_pa @@ -1631,7 +1631,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_effect_wspeed_pa => this%hvars(ih_effect_wspeed_pa)%r81d, & hio_fire_intensity_pa => this%hvars(ih_fire_intensity_pa)%r81d, & hio_fire_area_pa => this%hvars(ih_fire_area_pa)%r81d, & - hio_scorch_height_pa => this%hvars(ih_scorch_height_pa)%r81d, & + ! hio_scorch_height_pa => this%hvars(ih_scorch_height_pa)%r81d, & hio_fire_fuel_bulkd_pa => this%hvars(ih_fire_fuel_bulkd_pa)%r81d, & hio_fire_fuel_eff_moist_pa => this%hvars(ih_fire_fuel_eff_moist_pa)%r81d, & hio_fire_fuel_sav_pa => this%hvars(ih_fire_fuel_sav_pa)%r81d, & @@ -2438,7 +2438,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_tfc_ros_pa(io_pa) = cpatch%TFC_ROS hio_fire_intensity_pa(io_pa) = cpatch%FI hio_fire_area_pa(io_pa) = cpatch%frac_burnt - hio_scorch_height_pa(io_pa) = cpatch%SH + !hio_scorch_height_pa(io_pa) = cpatch%SH hio_fire_fuel_bulkd_pa(io_pa) = cpatch%fuel_bulkd hio_fire_fuel_eff_moist_pa(io_pa) = cpatch%fuel_eff_moist hio_fire_fuel_sav_pa(io_pa) = cpatch%fuel_sav @@ -3906,10 +3906,10 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_fire_area_pa ) - call this%set_history_var(vname='SCORCH_HEIGHT', units='m', & - long='spitfire flame height:m', use_default='active', & - avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & - ivar=ivar, initialize=initialize_variables, index = ih_scorch_height_pa ) + !call this%set_history_var(vname='SCORCH_HEIGHT', units='m', & + ! long='spitfire flame height:m', use_default='active', & + ! avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ! ivar=ivar, initialize=initialize_variables, index = ih_scorch_height_pa ) call this%set_history_var(vname='fire_fuel_mef', units='m', & long='spitfire fuel moisture', use_default='active', & diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 41a49c5e42..87e1659863 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -122,7 +122,7 @@ variables: fates_displar:long_name = "Ratio of displacement height to canopy top height" ; double fates_fire_alpha_SH(fates_pft) ; fates_fire_alpha_SH:units = "NA" ; - fates_fire_alpha_SH:long_name = "spitfire parameter, alpha scorch height, Equation 16 Thonicke et al 2010" ; + fates_fire_alpha_SH:long_name = "spitfire parameter, alpha scorch height, EQ 16 Thonicke et al 2010" ; double fates_fire_bark_scaler(fates_pft) ; fates_fire_bark_scaler:units = "fraction" ; fates_fire_bark_scaler:long_name = "the thickness of a cohorts bark as a fraction of its dbh" ; @@ -131,7 +131,10 @@ variables: fates_fire_crown_depth_frac:long_name = "the depth of a cohorts crown as a fraction of its height" ; double fates_fire_crown_kill(fates_pft) ; fates_fire_crown_kill:units = "NA" ; - fates_fire_crown_kill:long_name = "fire parameter, see equation 22 in Thonicke et al 2010" ; + fates_fire_crown_kill:long_name = "crown resistance to fire (1 = none), EQ 22 in Thonicke et al 2010" ; + double fates_fire_active_crown_fire(fates_pft) ; + fates_fire_active_crown_fire:units = "fraction" ; + fates_fire_active_crown_fire:long_name = "resistance to active crown fire (1=none)"; double fates_fr_fcel(fates_pft) ; fates_fr_fcel:units = "fraction" ; fates_fr_fcel:long_name = "Fine root litter cellulose fraction" ; @@ -718,6 +721,8 @@ data: fates_fire_crown_kill = 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775 ; + fates_fire_active_crown_fire = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + fates_fr_fcel = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 ; fates_fr_flab = 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25,