diff --git a/biogeochem/FatesCohortMod.F90 b/biogeochem/FatesCohortMod.F90 index 84ad0f2247..c541743bdc 100644 --- a/biogeochem/FatesCohortMod.F90 +++ b/biogeochem/FatesCohortMod.F90 @@ -279,6 +279,8 @@ module FatesCohortMod real(r8) :: rx_cambial_mort ! cambial kill mortality due to prescribed fire real(r8) :: rx_crown_mort ! crown fire mortality due to prescribed fire real(r8) :: rx_fire_mort ! post-fire mortality due to prescribed fire + real(r8) :: lfmc ! live fuel moisture content [%] + real(r8) :: canopy_fuel_1h ! leaf + 1 hour woody live fuel [kg biomass] !--------------------------------------------------------------------------- @@ -464,6 +466,8 @@ subroutine NanValues(this) this%rx_cambial_mort = nan this%rx_crown_mort = nan this%rx_fire_mort = nan + this%lfmc = nan + this%canopy_fuel_1h = nan end subroutine NanValues @@ -556,6 +560,8 @@ subroutine ZeroValues(this) this%rx_cambial_mort = 0._r8 this%rx_crown_mort = 0._r8 this%rx_fire_mort = 0._r8 + this%lfmc = 0._r8 + this%canopy_fuel_1h = 0._r8 end subroutine ZeroValues @@ -807,6 +813,8 @@ subroutine Copy(this, copyCohort) copyCohort%rx_cambial_mort = this%rx_cambial_mort copyCohort%rx_crown_mort = this%rx_crown_mort copyCohort%rx_fire_mort = this%rx_fire_mort + copyCohort%lfmc = this%lfmc + copyCohort%canopy_fuel_1h = this%canopy_fuel_1h ! HYDRAULICS if (hlm_use_planthydro .eq. itrue) then @@ -1151,6 +1159,8 @@ subroutine Dump(this) write(fates_log(),*) 'cohort%rx_crown_mort = ', this%rx_crown_mort write(fates_log(),*) 'cohort%rx_cambial_mort = ', this%rx_cambial_mort write(fates_log(),*) 'cohort%rx_fire_mort = ', this%rx_fire_mort + write(fates_log(),*) 'cohort%lfmc = ', this%lfmc + write(fates_log(),*) 'cohort%canopy_fuel_1h = ', this%canopy_fuel_1h write(fates_log(),*) 'cohort%size_class = ', this%size_class write(fates_log(),*) 'cohort%size_by_pft_class = ', this%size_by_pft_class diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 4efea0133b..79b152ab06 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -247,6 +247,11 @@ module FatesPatchMod ! fire effects real(r8) :: scorch_ht(maxpft) ! scorch height [m] real(r8) :: tfc_ros ! total intensity-relevant fuel consumed - no trunks [kgC/m2 of burned ground/day] + + ! crown fire + integer :: passive_crown_fire ! is there a passive crown fire [1=yes; 0=no] + integer :: active_crown_fire ! is there an active crown fire [1=yes; 0=no] + !--------------------------------------------------------------------------- ! PLANT HYDRAULICS (not currently used in hydraulics RGK 03-2018) @@ -537,6 +542,8 @@ subroutine NanValues(this) this%fire = fates_unset_int this%nonrx_fire = fates_unset_int this%rx_fire = fates_unset_int + this%passive_crown_fire = fates_unset_int + this%active_crown_fire = fates_unset_int this%nonrx_fi = nan this%nonrx_frac_burnt = nan this%rx_fi = nan diff --git a/fire/CMakeLists.txt b/fire/CMakeLists.txt index b018217dd8..dc27e3c0d1 100644 --- a/fire/CMakeLists.txt +++ b/fire/CMakeLists.txt @@ -6,6 +6,7 @@ list(APPEND fates_sources FatesFuelMod.F90 FatesFuelClassesMod.F90 SFEquationsMod.F90 + CrownFireEquationsMod.F90 ) sourcelist_to_parent(fates_sources) diff --git a/fire/CrownFireEquationsMod.F90 b/fire/CrownFireEquationsMod.F90 new file mode 100644 index 0000000000..72f82be061 --- /dev/null +++ b/fire/CrownFireEquationsMod.F90 @@ -0,0 +1,406 @@ +module CrownFireEquationsMod + + ! ============================================================================ + ! Helper methods for FATES crown fire model + ! Most are from Scott & Reinhardt 2001 & + ! Giuseppe 2024 + ! ============================================================================ + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : nearzero + + + implicit none + private + + public :: PassiveCrownFireIntensity + public :: HeatReleasePerArea + public :: CrowningIndex + public :: CrownFireIntensity + public :: LiveFuelMoistureContent + public :: MaxHeight + public :: CrownFireBehaveFM10 + public :: BiomassBin + public :: CrownFireCFB + + contains + + real(r8) function PassiveCrownFireIntensity(canopy_base_height, canopy_water_content) + ! DESCRIPTION: + ! Calculate the energy threshold for igniting crown fuels [kW/m or kJ/m/s] + ! EQ. 11 in Scott & Reinhardt 2001 + + ! ARGUMENTS: + real(r8), intent(in) :: canopy_base_height ! canopy base height at which biomass density > minimum density 0.011 kg/m3 + real(r8), intent(in) :: canopy_water_content ! canopy water content [%] + + + ! Locals: + real(r8) :: crown_ignition_energy ! surface fire intensity required to ignite crown fuels [kJ/kg] + + ! 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 canopy_water_content parameter + + ! Eq 3 Van Wagner 1977, Eq 11 Scott & Reinhardt 2001 + crown_ignition_energy = 460.0_r8 + 25.9_r8 * canopy_water_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) + PassiveCrownFireIntensity = (0.01_r8 * canopy_base_height * crown_ignition_energy)**1.5_r8 + + end function PassiveCrownFireIntensity + + !--------------------------------------------------------------------------------------- + real(r8) function HeatReleasePerArea(SAV, i_r) + ! + ! DESCRIPTION: + ! Calculate heat release per unit area of surface fire (HPA, kJ/m2) + ! EQ. 2 in Scott & Reinhardt 2001 + + ! ARGUMENTS: + real(r8), intent(in) :: SAV ! fuel surface area to volume ratio [/cm] + real(r8), intent(in) :: i_r ! reaction intensity [kJ/m2/min] + + + ! Locals: + real(r8) :: time_r ! fire residence time [min] + if (SAV < nearzero) then + time_r = 0.0_r8 + else + time_r = 12.595_r8 / SAV ! EQ 3 in Scott & Reinhardt 2001 + end if + + HeatReleasePerArea = i_r * time_r + + end function HeatReleasePerArea + + + !--------------------------------------------------------------------------------------- + real(r8) function CrowningIndex(eps, q_ig, i_r, xi, beta_ratio, SAV, FBD, canopy_bulk_density) + ! + ! DESCRIPTION: + ! Calculate open wind speed [km/hr] at which a fully active crown fire is sustained + ! EQ. 19 in Scott & Reinhardt 2001 + ! note that U'active calculated in EQ 19 is the mid-flame wind speed, which is open wind speed * 0.4 + ! so CI = U' active / 0.4, 0.4 is the wind reduction factor + ! 1/54.683/0.4 = 0.0457 + ! XLG: currently we are ignoring slope effect + ! + ! ARGUMENTS: + real(r8), intent(in) :: eps ! effective heating number [unitless] + real(r8), intent(in) :: q_ig ! heat of preignition [kJ/kg] + real(r8), intent(in) :: i_r ! reaction intensity [kJ/m2/min] + real(r8), intent(in) :: xi ! propagating flux [unitless] + real(r8), intent(in) :: beta_ratio ! ratio of packing ratio to optimum packing ratio [0-1] + real(r8), intent(in) :: SAV ! average fuel bed SA:V [cm-1] + real(r8), intent(in) :: FBD ! average fuel bed bulk density [kg m-3] + real(r8), intent(in) :: canopy_bulk_density ! canopy fuel bulk density [kg biomass / m3] + + ! Locals: + real(r8) :: heat_sink ! energy required to ignite per unit volume fuel bed [kJ m-3] + real(r8) :: CI_temp ! temporary variables + real(r8) :: b,c,e ! constants for calculating wind factor + + heat_sink = FBD*eps*q_ig + + if(i_r <= nearzero .or. xi <= nearzero .or. canopy_bulk_density <= nearzero) then + CI_temp = 0.0_r8 + else + CI_temp = (3.0_r8/canopy_bulk_density*heat_sink) / (3.34_r8*i_r*xi) - 1.0_r8 + end if + b = 0.15988_r8*(SAV**0.54_r8) + c = 7.47_r8*(exp(-0.8711_r8*(SAV**0.55_r8))) + e = 0.715_r8*(exp(-0.01094_r8*SAV)) + CrowningIndex = (CI_temp/(c*(beta_ratio)**(-1.0_r8*e)))**(1.0_r8/b) * 0.0457_r8 + + end function CrowningIndex + + + !--------------------------------------------------------------------------------------- + + real(r8) function CrownFireIntensity(HPA, canopy_fuel,patch_area, CFB, ROS_final) + ! + ! Description + ! Calculate fire intentisy for crown fire using + ! EQ. 22 in Scott & Reinhardt 2001 + ! + ! ARGUMENTS: + real(r8), intent(in) :: HPA ! heat release per unit area [kJ/m2] + real(r8), intent(in) :: canopy_fuel ! canopy fuel load [kg biomass] + real(r8), intent(in) :: patch_area ! current patch area, to calculate fuel load density in kg/m2 + real(r8), intent(in) :: CFB ! crown fraction burnt [fraction] + real(r8), intent(in) :: ROS_final ! final rate of spread after a crown fire happens [m/min] + ! Locals: + real(r8), parameter :: H_canopy = 18000.0_r8 ! heat yield for canopy fuels [kJ/kg biomass] + + CrownFireIntensity = (HPA + (canopy_fuel / patch_area * H_canopy * CFB)) * & + ROS_final / 60.0_r8 + + end function CrownFireIntensity + + !--------------------------------------------------------------------------------------- + + real(r8) function LiveFuelMoistureContent(lai, smp, min_lfmc, coeff_lfmc, & + smp_alpha, lai_beta, gamma_int) + ! + ! DESCRIPTION + ! Calculates live fuel moisture content + ! using EQ. 4 in McNorton and Giuseppe 2024 'A global fuel characteristic model + ! and dataset for wildfire prediction' + ! LFMC = max_lfmc - min_lfmc * exp(-swc_alpha*swc + lai_beta*lai + gamma_int*swc*lai) + ! swc is mass-based soil water content in the original Eq., we now switch to soil matric potential (smp). + ! This switch require refitting this equation to min_lfmc + coeff_lfmc*exp(smp_alpha*smp + lai_beta*lai + + ! gamma_int*smp*lai) + ! we can make all parameters to be PFT specific + ! + ! ARGUMENTS: + real(r8), intent(in) :: lai ! cohort level leaf area index [m2 m-2] + real(r8), intent(in) :: smp ! soil matric potential [MPa] + real(r8), intent(in) :: min_lfmc ! min. live fuel moisure content [%] + real(r8), intent(in) :: coeff_lfmc ! this will adds to the min to increase LMFC based on soil water and plant phenology [unitelss] + real(r8), intent(in) :: smp_alpha ! model coef. associated with swc [unitless] + real(r8), intent(in) :: lai_beta ! model coef. associated with lai [unitless] + real(r8), intent(in) :: gamma_int ! model coef. for interaction effect of swc and lai on LFMC [unitless] + + ! Locals: + real(r8) :: effect_temp ! the temporary variable for calculating effect of SWC and LAI + + effect_temp = smp_alpha*smp + lai_beta*lai + gamma_int*smp*lai + LiveFuelMoistureContent = min_lfmc + coeff_lfmc*exp(effect_temp) + + + end function LiveFuelMoistureContent + + !--------------------------------------------------------------------------------------- + subroutine MaxHeight(height, max_height) + ! + ! DESCRIPTION + ! Search for max height across cohorts on patch + ! + ! ARGUMENTS: + real(r8), intent(in) :: height ! cohort height [m] + real(r8), intent(inout) :: max_height ! max height of all cohorts [m] + + if(height > max_height)then + max_height = height + end if + + end subroutine MaxHeight + + + !--------------------------------------------------------------------------------------- + + + subroutine CrownFireBehaveFM10(fireWeatherClass, drying_ratio, wind, canopy_bulk_density, ROS_active, CI) + ! + ! DESCRIPTION + ! Calculate theoretical rate of spread for a active crown fire using + ! fuel model 10 and Rothermel's ROS model + ! + use SFEquationsMod, only : OptimumPackingRatio, ReactionIntensity + use SFEquationsMod, only : HeatofPreignition, EffectiveHeatingNumber + use SFEquationsMod, only : WindFactor, PropagatingFlux + use SFEquationsMod, only : ForwardRateOfSpread + use FatesFuelMod, only : fuel_type + use SFFireWeatherMod, only : fire_weather + use FatesFuelClassesMod, only : num_fuel_classes, fuel_classes + use SFParamsMod, only : SF_val_part_dens, SF_val_miner_total + + ! ARGUMENTS: + class(fire_weather), intent(in) :: fireWeatherClass ! fireWeatherClass + real(r8), intent(in) :: drying_ratio ! SPITFIRE fuel parameters controlling fuel dying rate [unitless] + real(r8), intent(in) :: wind ! Site wind speed [m/s] + real(r8), intent(in) :: canopy_bulk_density ! Canopy fuel bulk density [kg biomass / m3] ! + real(r8), intent(out) :: ROS_active ! Theoretical rate of spread a fully active crown fire using fuel model 10 [m/min] + real(r8), intent(out) :: CI ! Open wind speed to sustain an active crown fire using fuel model 10 [m/min] + + ! Local variables: + type(fuel_type) :: fuel_fm10 ! fuel object initialized using fuel model 10 data + real(r8) :: fuel_1h ! 1 hour fuel load using FM 10 [kg] + real(r8) :: fuel_10h ! 10 hour fuel load using FM 10 [kg] + real(r8) :: fuel_100h ! 100 hour fuel load using FM 10 [kg] + real(r8) :: fuel_live ! live fuel load using FM 10 [kg] + real(r8) :: fuel_depth ! fuel bed depth using FM 10 [m] + real(r8) :: fuel_bd ! fuel bulk density using FM 10 [kg biomass/m3] + real(r8) :: fuel_sav1h ! SAV for 1 hour fuel for FM 10 [/cm] + real(r8) :: fuel_sav10h ! SAV for 10 hour fuel for FM 10 [/cm] + real(r8) :: fuel_sav100h ! SAV for 100 hour fuel for FM 10 [/cm] + real(r8) :: fuel_savlive ! SAV for live fuel for FM 10 [/cm] + real(r8) :: fuel_sav(num_fuel_classes) ! fuel surface area to volume ratio by fuel size [/cm] + real(r8) :: midflame_wind ! 40% of open wind speed + real(r8) :: beta_fm10 ! packing ratio derived for fuel model 10 [unitless] + real(r8) :: beta_op_fm10 ! optimum packing ratio for FM 10 [unitless] + real(r8) :: beta_ratio_fm10 ! relative packing ratio for FM 10 [unitless] + real(r8) :: ir_dead ! reaction intensity of dead fuel [kJ m-2 min-1] + real(r8) :: ir_live ! reaction intensity of live fuel [kJ m-2 min-1] + real(r8) :: i_r_fm10 ! reaction intensity for FM 10 [kJ/m2/min] + real(r8) :: xi_fm10 ! propagating flux ratio for FM 10 [unitless] + real(r8) :: phi_wind_fm10 ! wind factor for FM 10 [unitless] + real(r8) :: q_ig_fm10 ! heat of pre-ignition for FM 10 [kJ/kg] + real(r8) :: eps_fm10 ! effective heating number for FM 10 [unitless] + + ! Parameters for fuel model 10 to describe fuel characteristics; and some constants + ! fuel loading, MEF, and depth from Anderson 1982 Aids to determining fuel models for fire behavior + ! SAV values from BEHAVE model Burgan & Rothermel (1984) + ! XLG: Can consider change all FM 10 parameters to SF parameters so users can use costumized fuel model for their study + + real(r8),parameter :: fuel_1h_ton = 3.01_r8 ! FM 10 1-hr fuel loading (US tons/acre) + real(r8),parameter :: fuel_10h_ton = 2.0_r8 ! FM 10 10-hr fuel loading (US tons/acre) + real(r8),parameter :: fuel_100h_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), XLG: should we use this from FM 10?? + real(r8),parameter :: fuel_depth_ft= 1.0_r8 ! FM 10 fuel depth (ft) + real(r8),parameter :: sav_1h_ft = 2000.0_r8 ! BEHAVE model 1-hr SAV (ft2/ft3) + real(r8),parameter :: sav_10h_ft = 109.0_r8 ! BEHAVE model 10-hr SAV (ft2/ft3) + real(r8),parameter :: sav_100h_ft = 30.0_r8 ! BEHAVE model 100-hr SAV (ft2/ft3) + real(r8),parameter :: sav_live_ft = 1650.0_r8 ! BEHAVE model live SAV (ft2/ft3) + real(r8),parameter :: tonnes_acre_to_kg_m2 = 0.2241701_r8 ! convert tons/acre to kg/m2 + real(r8),parameter :: sqft_cubicft_to_sqm_cubicm = 0.03280844_r8 ! convert ft2/ft3 to m2/m3 + real(r8),parameter :: ft_to_meter = 0.3048_r8 ! convert ft to meter + real(r8),parameter :: km_per_hr_to_m_per_min = 16.6667_r8 ! convert km/hour to m/min for wind speed + real(r8),parameter :: mef_a = 0.524_r8 + real(r8),parameter :: mef_b = 0.066_r8 + + ! set up fuels. Since currently FATES fuel classes don't include + ! live woody fuels, we assign the live fuels from FM10 (only live woody no live grass) + ! to live grass. Also, 1h to dead leaf, 10h to small branch, and 100h to large branch + + fuel_1h = fuel_1h_ton * tonnes_acre_to_kg_m2 + fuel_10h = fuel_10h_ton * tonnes_acre_to_kg_m2 + fuel_100h = fuel_100h_ton * tonnes_acre_to_kg_m2 + fuel_live = fuel_live_ton * tonnes_acre_to_kg_m2 + fuel_sav1h = sav_1h_ft * sqft_cubicft_to_sqm_cubicm + fuel_sav10h = sav_10h_ft * sqft_cubicft_to_sqm_cubicm + fuel_sav100h = sav_100h_ft * sqft_cubicft_to_sqm_cubicm + fuel_savlive = sav_live_ft * sqft_cubicft_to_sqm_cubicm + + ! since live fuel moisture is calculated using SA:V of twig, we assign fuel_savlive to twig too + ! SA:V cannot be zero when used for calculation of moisture, so we set trunk SA:V to a small value + fuel_sav = (/fuel_savlive, fuel_sav10h, fuel_sav100h, 0.002_r8, fuel_sav1h, fuel_savlive/) + + ! update fuel chracteristics + call fuel_fm10%Init() + call fuel_fm10%UpdateLoading(fuel_1h, 0.0_r8, fuel_10h, fuel_100h, 0.0_r8, fuel_live) + call fuel_fm10%SumLoading() + call fuel_fm10%CalculateFractionalLoading() + call fuel_fm10%UpdateFuelMoisture(fuel_sav, drying_ratio, fireWeatherClass) + call fuel_fm10%AverageSAV_NoTrunks(fuel_sav) + + ! use total fuel and fuel bed depth to calculate fuel bulk density + fuel_depth = fuel_depth_ft * ft_to_meter !convert to meters + fuel_bd = fuel_fm10%non_trunk_loading/fuel_depth !fuel bulk density (kg biomass/m3) + + ! remove mineral content, this fuel is in kg biomass m-2 + fuel_fm10%non_trunk_loading = fuel_fm10%non_trunk_loading*(1.0_r8 - SF_val_miner_total) + + beta_fm10 = fuel_bd / SF_val_part_dens + beta_op_fm10 = OptimumPackingRatio(fuel_fm10%SAV_notrunks) + if(beta_op_fm10 < nearzero) then + beta_ratio_fm10 = 0.0_r8 + else + beta_ratio_fm10 = beta_fm10 / beta_op_fm10 + end if + + i_r_fm10 = ReactionIntensity(fuel_fm10%non_trunk_loading, fuel_fm10%SAV_notrunks, & + beta_ratio_fm10, fuel_fm10%average_moisture_notrunks, fuel_fm10%MEF_notrunks) + + ! XLG: I used calculated MEF instead of the constant MEF from fm10 + q_ig_fm10 = HeatofPreignition(fuel_fm10%average_moisture_notrunks) + + eps_fm10 = EffectiveHeatingNumber(fuel_fm10%SAV_notrunks) + + ! Scott & Reinhardt 2001 use 40% of open wind speed as effective wind speed + midflame_wind = wind * 0.40_r8 + phi_wind_fm10 = WindFactor(midflame_wind, beta_ratio_fm10, fuel_fm10%SAV_notrunks) + + xi_fm10 = PropagatingFlux(beta_fm10, fuel_fm10%SAV_notrunks) + + ! Calculate ROS_active, used for determining whether there is active or passtive crown fire + ROS_active = ForwardRateOfSpread(fuel_bd, eps_fm10, q_ig_fm10, i_r_fm10, & + xi_fm10, phi_wind_fm10) + ! apply the multiplier of 3.34 for the correlation between durface ROS and active crown fire ROS + ! EQ. 8 + ROS_active = ROS_active * 3.34_r8 + + ! Calculate crowning index, which is used for calculating ROS_SA + CI = CrowningIndex(eps_fm10, q_ig_fm10, i_r_fm10, xi_fm10, beta_ratio_fm10, & + fuel_fm10%SAV_notrunks, fuel_bd, canopy_bulk_density) + CI = CI * km_per_hr_to_m_per_min ! convert to m/min + + end subroutine CrownFireBehaveFM10 + + !--------------------------------------------------------------------------------------- + + + subroutine BiomassBin(cbh, height, crown_depth, canopy_fuel_1h, biom_matrix) + ! + ! DESCRIPTION: + ! accumulate biomass at 1 m interval across cohorts and return + ! the biomass array + ! + ! ARGUMENTS: + real(r8), intent(in) :: cbh ! canopy base height [m] + real(r8), intent(in) :: height ! cohort height [m] + real(r8), intent(in) :: crown_depth ! crown length of a cohort [m] + real(r8), intent(in) :: canopy_fuel_1h ! leaf + 1-hour woody biomass of a cohort [kg biomass] + real(r8), intent(inout) :: biom_matrix(0:) ! array that holds biomass by 1m interval across all cohorts on a patch [kg biomass] + + ! LOCALS: + integer :: h_idx ! looping index + real(r8) :: crown_fuel_per_m ! biomass by 1m interval [kg biomass] + + ! calculate biomass by 1m interval + crown_fuel_per_m = canopy_fuel_1h / crown_depth + + do h_idx = int(cbh), int(height) + biom_matrix(h_idx) = biom_matrix(h_idx) + crown_fuel_per_m + end do + + + end subroutine BiomassBin + + !--------------------------------------------------------------------------------------- + + subroutine CrownFireCFB(ROS_active, ROS_critical, ROS_front, & + ROS_init, ROS_SA, active_crownfire, passive_crownfire, crown_frac_burnt) + ! + ! DESCRIPTION: + ! Determine if it is passive or active crown fire and + ! return crown fraction burned EQ. 28 in Scott & Reinhardt 2001 + ! + ! ARGUMENTS: + real(r8), intent(in) :: ROS_active ! active crown fire ROS using FM10 [m/min] + real(r8), intent(in) :: ROS_critical ! critical ROS for sustaining active crown fire [m/min] + real(r8), intent(in) :: ROS_front ! surface ROS [m/min] + real(r8), intent(in) :: ROS_init ! ROS to initiate a crown fire [m/min] + real(r8), intent(in) :: ROS_SA ! ROS at crowning index wind speed when ROS_active = ROS_SA [m/min] + integer, intent(out) :: active_crownfire ! active crown fire = 1 + integer, intent(out) :: passive_crownfire ! passive crown fire = 1 + real(r8), intent(out) :: crown_frac_burnt ! crown fraction burnt [0-1] + ! + if(ROS_active >= ROS_critical)then + active_crownfire = 1 + passive_crownfire = 0 + ! for active crown fire we set CFB to 1 + crown_frac_burnt = 1.0_r8 + else if(ROS_active < ROS_critical .and. & + ROS_front > ROS_init .and. ROS_front < ROS_SA)then ! calculation of ROSs are disconnected so check to make sure these are true + active_crownfire = 0 + passive_crownfire = 1 + ! calculate crown fraction burnt EQ. 28 in Scott & Reinhardt 2001 + crown_frac_burnt = min(1.0_r8, (ROS_front - ROS_init) / (ROS_SA - ROS_init)) + else + active_crownfire = 0 + passive_crownfire = 0 + crown_frac_burnt = 0.0_r8 + end if + + + end subroutine CrownFireCFB + + +end module CrownFireEquationsMod \ No newline at end of file diff --git a/fire/FatesFuelMod.F90 b/fire/FatesFuelMod.F90 index b05702fc1a..44f3281b96 100644 --- a/fire/FatesFuelMod.F90 +++ b/fire/FatesFuelMod.F90 @@ -5,6 +5,7 @@ module FatesFuelMod use FatesConstantsMod, only : nearzero use SFNesterovMod, only : nesterov_index use SFFireWeatherMod, only : fire_weather + use FatesLitterMod, only : ncwd use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun use shr_log_mod, only : errMsg => shr_log_errMsg @@ -23,6 +24,11 @@ module FatesFuelMod real(r8) :: bulk_density_notrunks ! weighted average of bulk density across non-trunk fuel classes [kg/m3] real(r8) :: SAV_notrunks ! weighted average of surface area to volume ratio across non-trunk fuel classes [/cm] real(r8) :: MEF_notrunks ! weighted average of moisture of extinction across non-trunk fuel classes [m3/m3] + real(r8) :: canopy_fuel_load ! patch level total canopy fuel load [kg biomass] + real(r8) :: canopy_base_height ! patch level canopy base height at which biomass density > minimum density required for canopy fire spread [m] + real(r8) :: canopy_bulk_density ! patch level canopy fuel bulk density [kg biomass m-3] + real(r8) :: canopy_water_content ! patch level canopy water content [%] + contains @@ -34,8 +40,11 @@ module FatesFuelMod procedure :: UpdateFuelMoisture procedure :: AverageBulkDensity_NoTrunks procedure :: AverageSAV_NoTrunks + procedure :: CalculateCanopyFuelLoad + procedure :: CalculateCanopyBulkDensity procedure :: CalculateFuelBurnt procedure :: CalculateResidenceTime + procedure :: CanopyWaterContent end type fuel_type @@ -58,6 +67,10 @@ subroutine Init(this) this%bulk_density_notrunks = 0.0_r8 this%SAV_notrunks = 0.0_r8 this%MEF_notrunks = 0.0_r8 + this%canopy_fuel_load = 0.0_r8 + this%canopy_base_height = 0.0_r8 + this%canopy_bulk_density = 0.0_r8 + this%canopy_water_content = 0.0_r8 end subroutine Init @@ -99,7 +112,15 @@ subroutine Fuse(this, self_area, donor_area, donor_fuel) this%bulk_density_notrunks = this%bulk_density_notrunks*self_weight + & donor_fuel%bulk_density_notrunks*donor_weight this%SAV_notrunks = this%SAV_notrunks*self_weight + donor_fuel%SAV_notrunks*donor_weight - this%MEF_notrunks = this%MEF_notrunks*self_weight + donor_fuel%MEF_notrunks*donor_weight + this%MEF_notrunks = this%MEF_notrunks*self_weight + donor_fuel%MEF_notrunks*donor_weight + this%canopy_fuel_load = this%canopy_fuel_load*self_weight + & + donor_fuel%canopy_fuel_load*donor_weight + this%canopy_base_height = this%canopy_base_height*self_weight + & + donor_fuel%canopy_base_height*donor_weight + this%canopy_bulk_density = this%canopy_bulk_density*self_weight + & + donor_fuel%canopy_bulk_density*donor_weight + this%canopy_water_content = this%canopy_water_content*self_weight + & + donor_fuel%canopy_water_content*donor_weight end subroutine Fuse @@ -377,7 +398,117 @@ subroutine AverageSAV_NoTrunks(this, sav_fuel) end subroutine AverageSAV_NoTrunks !--------------------------------------------------------------------------------------- - + + subroutine CalculateCanopyFuelLoad(this, leaf_c, woody_c, & + cwd_frac_adj, canopy_fuel_1h) + ! DESCRIPTION: + ! Calculate canopy fuel load by summing up leaf biomass and 1 hour woody biomass + ! XLG: I did not change how Sam L. calculate canopy fuel load here except + ! switch to the adjusted CWD instead of using the default parameter + ! Seems only including leaf biomass is recommend by Gruz et al. 2003: + ! https://www.sierraforestlegacy.org/Resources/Conservation/FireForestEcology/FireScienceResearch/FireEcology/FireEcology-Cruz03.pdf + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: leaf_c ! leaf biomass of each cohort [kg biomass] + real(r8), intent(in) :: woody_c ! 1-hour woody fuels [kg biomass] + real(r8), intent(in) :: cwd_frac_adj(ncwd) ! adjusted coarse woody debris fraction [fraction] + real(r8), intent(out) :: canopy_fuel_1h ! leaf + 1-hour woody biomass of current cohort [kg biomass] + + ! calculate canopy fuel load for FATES crown fire + canopy_fuel_1h = woody_c * cwd_frac_adj(1) + leaf_c + this%canopy_fuel_load = this%canopy_fuel_load + canopy_fuel_1h + + + end subroutine CalculateCanopyFuelLoad + + + !--------------------------------------------------------------------------------------- + subroutine CalculateCanopyBulkDensity(this, biom_matrix, max_height) + ! DESCRIPTION: + ! Calculate canopy fuel bulk density [kg biomass / m3] by searching for + ! the canopy height (canopy base height CBH) at which canopy fuel bed density + ! is > minimum canopy fuel density [0.011 kg biomass/m3] for propogating canopy fire; + ! CBD is then calculated as canopy fuel load devided by + ! crown depth ( CD= max. canopy height - CBH). + ! XLG: one modification I made compare to Sam L. version is that I recalculate + ! canopy fuel load after knowing CD by excluding fuels below CBH. As if including + ! fuels below CBH then we are artificially "condense" canopy fuel bed by deviding total + ! fuel load by an average, shallower CD. + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: biom_matrix(0:) ! 1m biomass bin (kg/m3) in the vertical space to sort the canopy fule + ! across cohorts, used to search for CBH + real(r8), intent(in) :: max_height ! the max. cohort height at current patch [m] + + ! Locals: + real(r8) :: canopy_top_height ! the highest point at which biomass density is > minimum canopy fuel density [m] + integer :: ih ! biomass bin counter + + ! Params: + real(r8), parameter :: min_density_canopy_fuel = 0.011_r8 ! minimum canopy fuel density for propogating canopy fire + ! vertically through canopy. + ! Scott and Reinhardt 2001 RMRS-RP-29 + + ! loop from 1m to 70m to find CBH + do ih=0,int(max_height) + if (biom_matrix(ih) > min_density_canopy_fuel) then + this%canopy_base_height = dble(ih) + 1.0_r8 ! the dimension index of biom_matrix is a ronded-down integer of cohort height + ! add 1 to be conservative when searching for CBH + + exit + else ! when this is an open stand that there is no such a height, use max_height + this%canopy_base_height = max_height + end if + end do + + ! now loop from top to bottom to find the highest point at which the minimum bulk density is met + ! XLG: I modified the way how canopy bulk density is calculated by reducing the maximum canopy height + ! to the highest point where the minimum bulk density is met. + + do ih = int(max_height), 0, -1 + if (biom_matrix(ih) > min_density_canopy_fuel) then + canopy_top_height = dble(ih) + 1.0_r8 + exit + else + canopy_top_height = max_height + end if + end do + + ! XLG: We now only calculate canopy bulk density for fuels between + ! canopy base and top height + if ((canopy_top_height - this%canopy_base_height) > nearzero) then + this%canopy_bulk_density = sum(biom_matrix(int(this%canopy_base_height-1.0_r8):int(canopy_top_height-1.0_r8))) / & + (canopy_top_height - this%canopy_base_height) + else + this%canopy_bulk_density = 0.0_r8 + end if + + + end subroutine CalculateCanopyBulkDensity + + !--------------------------------------------------------------------------------------- + + subroutine CanopyWaterContent(this, co_cwc, co_fuel) + ! DESCRIPTION: + ! Calculates live canopy water content + ! as sum of fuel load weighted cohort level live fuel moisture content + ! + ! ARGUMENTS + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: co_cwc ! coohort live fuel moisture content [%] + real(r8), intent(in) :: co_fuel ! cohort canopy fuel, only 1 hour woody + leaf + + this%canopy_water_content = this%canopy_water_content + & + co_cwc * co_fuel/this%canopy_fuel_load + + + end subroutine CanopyWaterContent + + !--------------------------------------------------------------------------------------- + + subroutine CalculateFuelBurnt(this, fuel_consumed) ! DESCRIPTION: ! Calculates the fraction and total amount of fuel burnt diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 7e6de3323b..479ee87407 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -29,6 +29,9 @@ module SFMainMod use FatesLitterMod, only : litter_type use FatesFuelClassesMod, only : num_fuel_classes use PRTGenericMod, only : carbon12_element + use PRTGenericMod, only : leaf_organ + use PRTGenericMod, only : sapw_organ + use PRTGenericMod, only : struct_organ use FatesInterfaceTypesMod, only : numpft use FatesAllometryMod, only : CrownDepth use FatesFuelClassesMod, only : fuel_classes @@ -58,9 +61,11 @@ subroutine DailyFireModel(currentSite, bc_in) if (hlm_spitfire_mode > hlm_sf_nofire_def) then call UpdateFireWeather(currentSite, bc_in) call UpdateFuelCharacteristics(currentSite) + call UpdateCanopyFuelCharacteristics(currentSite, bc_in) call CalculateIgnitionsandFDI(currentSite, bc_in) call CalculateSurfaceRateOfSpread(currentSite) call CalculateSurfaceFireIntensity(currentSite) + call PassiveActiveCrownFireCheck(currentSite) call CalculateAreaBurnt(currentSite) call CalculateRxFireAreaBurnt(currentSite) call CalculatePostFireMortality(currentSite) @@ -148,7 +153,7 @@ subroutine UpdateFuelCharacteristics(currentSite) ! Updates fuel characteristics on each patch of the site ! - 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 ! ARGUMENTS: type(ed_site_type), intent(in), target :: currentSite ! site object @@ -190,6 +195,182 @@ subroutine UpdateFuelCharacteristics(currentSite) end subroutine UpdateFuelCharacteristics + !--------------------------------------------------------------------------------------- + + subroutine UpdateCanopyFuelCharacteristics(currentSite, bc_in) + ! + ! DESCRIPTION: + ! Calculate canopy fuel load (sum of all 1 hour live fuels, in kg biomass), + ! canopy base height (minimum canopy height at which canopy fuel bed density > 0.011 kg/m3 for + ! vertical fire spread through canopy), and canopy bulk desity (kg biomass/m3) + ! + use FatesLitterMod, only : ncwd + use SFParamsMod, only : SF_val_CWD_frac + use FatesLitterMod, only : adjust_SF_CWD_frac + use CrownFireEquationsMod, only : MaxHeight + use CrownFireEquationsMod, only : BiomassBin + use CrownFireEquationsMod, only : LiveFuelMoistureContent + use EDTypesMod, only : numWaterMem + use FatesInterfaceTypesMod, only : hlm_use_planthydro + use FatesConstantsMod, only : denh2o => dens_fresh_liquid_water + use FatesHydraulicsMemMod, only : n_hypool_leaf + use FatesHydraulicsMemMod, only : ed_cohort_hydr_type + + ! ARGUMENTS: + type(ed_site_type), intent(inout), target :: currentSite ! site object + type(bc_in_type), intent(in) :: bc_in + + + type(fates_patch_type), pointer :: currentPatch ! FATES patch + type(fates_cohort_type), pointer :: currentCohort ! FATES cohort + type(ed_cohort_hydr_type), pointer :: ccohort_hydr ! current cohort hyro-derived type + + ! Locals: + real(r8) :: crown_depth ! depth of crown [m] + real(r8) :: cbh_co ! canopy base height of cohort [m] + real(r8) :: max_height ! max cohort height on patch (m) + real(r8) :: woody_c ! above-ground tree struct and sapwood 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) :: canopy_fuel_1h ! leaf + 1-hour woody fuel [kg biomass] + real(r8) :: SF_val_CWD_frac_adj(ncwd) ! adjusted fractional allocation of woody biomass to coarse wood debris pool + real(r8) :: mean_10day_smp(numpft) ! averaged 10 day soil matric potential for each PFT + real(r8) :: water_mass ! canopy water mass per individual plant (kgH2O/plant) + integer :: ipft ! pft index + real(r8), allocatable :: biom_matrix(:) ! matrix to track biomass from bottom to top + + ! Parameters + real(r8), parameter :: carbon_2_biomass = 0.45_r8 + ! LFMC parameters for testing + real(r8), parameter :: min_lfmc = 80.0_r8 + real(r8), parameter :: coeff_lfmc = 30.0_r8 + real(r8), parameter :: smp_alpha = 4E-4_r8 + real(r8), parameter :: lai_beta = 0.15_r8 + real(r8), parameter :: gamma_int = 0.0_r8 + + ! update site level soil matric potential for each PFT + mean_10day_smp(:) = 0.0_r8 + + do ipft=1,numpft + if(int(prt_params%woody(ipft)) == itrue) then + mean_10day_smp(ipft) = sum(currentSite%smp_memory(1:numWaterMem,ipft)) / & + real(numWaterMem,r8) + + end if + end do + + currentPatch => currentSite%oldest_patch + + do while(associated(currentPatch)) + if (currentPatch%nocomp_pft_label /= nocomp_bareground) then + !zero Patch level variables + max_height = 0.0_r8 + currentPatch%fuel%canopy_fuel_load = 0.0_r8 + currentPatch%fuel%canopy_bulk_density = 0.0_r8 + currentPatch%fuel%canopy_base_height = 0.0_r8 + currentPatch%fuel%canopy_water_content = 0.0_r8 + + ! find the max cohort height to set the upper bounds of biom_matrix + currentCohort=>currentPatch%tallest + do while(associated(currentCohort)) + + if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees + + call MaxHeight(currentCohort%height, max_height) + end if ! trees only + currentCohort => currentCohort%shorter; + end do ! end cohort loop + + !allocate and initialize biom_matrix + allocate(biom_matrix(0:int(max_height))) + biom_matrix(:) = 0.0_r8 + + + !loop across cohorts to calculate canopy fuel load by 1m height bin + currentCohort => currentPatch%tallest + + do while(associated(currentCohort)) + !zero cohort level variables + woody_c = 0.0_r8 + leaf_c = 0.0_r8 + sapw_c = 0.0_r8 + struct_c = 0.0_r8 + crown_depth = 0.0_r8 + cbh_co = 0.0_r8 + SF_val_CWD_frac_adj(ncwd) = 0.0_r8 + + + ! Calculate crown 1hr fuel biomass (leaf, twig sapwood, twig structural biomass) + if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees + ! calculate canopy base height for cohort, this is a different concept than the canopy fuel bed base height + call CrownDepth(currentCohort%height,currentCohort%pft,crown_depth) + cbh_co = currentCohort%height - crown_depth + + 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) + + woody_c = currentCohort%n * (prt_params%allom_agb_frac(currentCohort%pft)* & + (sapw_c + struct_c)) / carbon_2_biomass + leaf_c = currentCohort%n * leaf_c / carbon_2_biomass + + call adjust_SF_CWD_frac(currentCohort%dbh, ncwd, SF_val_CWD_frac, SF_val_CWD_frac_adj) + + ! update canopy fuel load + call currentPatch%fuel%CalculateCanopyFuelLoad(leaf_c, woody_c, & + SF_val_CWD_frac_adj, canopy_fuel_1h) + currentCohort%canopy_fuel_1h = canopy_fuel_1h + ! 1m biomass bin + call BiomassBin(cbh_co, currentCohort%height, crown_depth, canopy_fuel_1h, biom_matrix) + + ! calculate live canopy water content for current cohort + ! either using statistic model- or hydro-derived cwc + if (hlm_use_planthydro == itrue) then + ccohort_hydr => currentCohort%co_hydr + water_mass = ccohort_hydr%th_ag(n_hypool_leaf) * ccohort_hydr%v_ag(n_hypool_leaf) * & + denh2o * currentCohort%n + currentCohort%lfmc = water_mass / leaf_c * 100.0_r8 ! leaf water content in percent + else + currentCohort%lfmc = LiveFuelMoistureContent(currentCohort%treelai, & + mean_10day_smp(currentCohort%pft), & + min_lfmc, coeff_lfmc, smp_alpha, lai_beta, gamma_int) + end if + write(fates_log(),*) 'current cohort cwc is ', currentCohort%lfmc + + end if ! trees only + currentCohort => currentCohort%shorter; + end do ! end cohort loop + write(fates_log(),*) 'current patch canopy fuel is ', currentPatch%fuel%canopy_fuel_load + + ! loop across cohorts to calculate patch level canopy water content + currentCohort => currentPatch%tallest + + do while(associated(currentCohort)) + if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees + call currentPatch%fuel%CanopyWaterContent(currentCohort%lfmc, & + currentCohort%canopy_fuel_1h) + end if + currentCohort => currentCohort%shorter; + end do + + biom_matrix(:) = biom_matrix(:) / currentPatch%area ! kg biomass / m3 + ! update canopy fuel bulk density and canopy base height + call currentPatch%fuel%CalculateCanopyBulkDensity(biom_matrix, max_height) + write(fates_log(),*) 'current patch CBD is ', currentPatch%fuel%canopy_bulk_density + write(fates_log(),*) 'current patch canopy base height is ', currentPatch%fuel%canopy_base_height + + deallocate(biom_matrix) + + end if ! nocomp_bareground + currentPatch => currentPatch%younger; + end do ! end patch loop + + end subroutine UpdateCanopyFuelCharacteristics + + + + !--------------------------------------------------------------------------------------- subroutine CalculateIgnitionsandFDI(currentSite, bc_in) @@ -443,6 +624,185 @@ subroutine CalculateSurfaceFireIntensity(currentSite) end subroutine CalculateSurfaceFireIntensity + !--------------------------------------------------------------------------------------- + + subroutine PassiveActiveCrownFireCheck(currentSite) + ! + ! DESCRIPTION: + ! 1) Calculate theoretical active crown fire ROS (R_active) using fuel model 10 and Rothermel's ROS model + ! XLG: the calculation of R_active can also incluede canopy bulk density (CBD) and canopy water content (see EQ. 10), + ! but those are currently ignored in Scott & Reinhardt 2001 due to contrasting effects of CBD on R_active in lit., + ! and the difficulty of obtaining a base-line FME (see discussion in Scott & Reinhardt 2001 pg12-13) + ! 2) Calculate the critical mass flow rate for sustaining a fully active crown fire given current + ! canopy bulk density, which is R'_active in Scott & Reinhardt 2001; + ! 3) Calculate the critical open wind speed for sustaining a fully active crown fire, AKA crowning index (CI) + ! using FM10 fuel characteristics + ! and the theoretical surface fire ROS (R'_SA, will be used in calculating final ROS for crown fire) using CI + ! XLG: I'm not sure what SAV value they used in Scott & Reinhardt 2001 to calculate CI, + ! Sam Levis used SAV values from the BEHAVE model, but when use that to calculate some constants and comapre + ! to EQ. 20 it doesn't seem right to me. This is something we have to discuss + ! 4) Calculate the critical surface fire ROS (R'_initiation) for initiating a crown fire + ! 5) Determine passive or active crown fire by comparing current surface ROS_front to R'_initiation, and + ! R_active to R'_active. Passive: ROS_front > R'_initiation but R_active < R'_active; + ! Active: ROS_front > R'_initiation and R_active > R'_active. Note: FI > FI_init is equivalent to ROS_front > R'_initiation + ! 6) Calculate new ROS and FI and update ROS_front and FI + ! + use SFParamsMod, only : SF_val_miner_total, SF_val_part_dens, SF_val_drying_ratio + use EDParamsMod, only : crown_fire_switch + use EDTypesMod, only : CalculateTreeGrassAreaSite + use SFEquationsMod, only : OptimumPackingRatio, ReactionIntensity + use SFEquationsMod, only : HeatofPreignition, EffectiveHeatingNumber + use SFEquationsMod, only : WindFactor, PropagatingFlux + use SFEquationsMod, only : ForwardRateOfSpread + use CrownFireEquationsMod, only : CrownFireBehaveFM10, CrownFireCFB + use CrownFireEquationsMod, only : PassiveCrownFireIntensity, HeatReleasePerArea + use CrownFireEquationsMod, only : CrowningIndex, CrownFireIntensity + + ! ARGUMENTS: + type(ed_site_type), intent(inout), target :: currentSite + + ! LOCALS: + type(fates_patch_type), pointer :: currentPatch ! patch object + + + real(r8) :: FI_init ! critical surface fire intensity for initiating a crown fire [kW/m or kJ/m/s] + real(r8) :: HPA ! heat release per unit area [kW/m2] + real(r8) :: ROS_init ! critical surface ROS for initiating a crown fire [m/min] + real(r8) :: ROS_active_min ! critical ROS for sustaining a fully active crown fire [m/min] + real(r8) :: ROS_active ! theoretical crown fire rate of spread using FM 10 fuels [m/min] + real(r8) :: ROS_final ! final ROS when a crown fire happens [m/min] + real(r8) :: FI_final ! final fireline intensity with crown consumption [kW/m or kJ/m/s] + real(r8) :: crown_frac_burnt ! patch level crown fraction burnt to calculate final crown fire intensty [fraction] + integer :: active_crownfire ! 1 = active crown fire + integer :: passive_crownfire ! 1 = passive crown fire + + + ! local variables for calculating ROS_SA based on current patch fuel condition + real(r8) :: ROS_SA ! surface ROS calculated at open wind speed of CI [m/min] + real(r8) :: beta ! packing ratio of current patch [unitless] + real(r8) :: beta_op ! optimum packing ratio of current patch [unitless] + real(r8) :: beta_ratio ! relative packing ratio of current patch [unitless] + real(r8) :: i_r ! reaction intensity of current patch [kJ/m2/min] + real(r8) :: xi ! propagating flux ratio of current patch [unitless] + real(r8) :: eps ! effective heating number of current patch [unitless] + real(r8) :: tree_fraction ! site-level tree fraction [0-1] + real(r8) :: grass_fraction ! site-level grass fraction [0-1] + real(r8) :: bare_fraction ! site-level bare ground fraction [0-1] + real(r8) :: CI ! crowning index: open wind speed to sustain active crown fire [km/hr] + real(r8) :: CI_effective ! effective wind speed using CI [m/min] + real(r8) :: phi_wind_SA ! wind factor for calculating ROS_SA [unitless] + real(r8) :: q_ig ! heat of pre-ignition of current patch [kJ/kg] + + ! Local parameters + real(r8), parameter :: wind_atten_tree = 0.4_r8 ! wind attenuation factor for tree fraction + real(r8), parameter :: wind_atten_grass = 0.6_r8 ! wind attenuation factor for grass fraction + + + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) + ! initialize patch level variables + currentPatch%passive_crown_fire = 0 + currentPatch%active_crown_fire = 0 + + if (currentPatch%nocomp_pft_label /= nocomp_bareground .and. & + currentPatch%fire == itrue .and. crown_fire_switch) then + ! calculate passive crown fire intensity, the minimum surface FI to initiate a crown fire + FI_init = PassiveCrownFireIntensity(currentPatch%fuel%canopy_base_height, & + currentPatch%fuel%canopy_water_content) + + write(fates_log(),*) 'FI_init is ', FI_init + write(fates_log(),*) 'FI is ', currentPatch%FI + + ! check if there is a crown fire + if (currentPatch%FI > FI_init ) then + ! calculate ROS_active and CI using fuel characteristics from fuel model 10 + call CrownFireBehaveFM10(currentSite%fireWeather, SF_val_drying_ratio, & + currentSite%wind, currentPatch%fuel%canopy_bulk_density, ROS_active, CI) + + ! Calculate ROS_acitive_min, EQ 14 in Scott & Reinhardt 2001 + ROS_active_min = 3.0_r8 / currentPatch%fuel%canopy_bulk_density + + ! Calculate ROS_SA using current patch fuel conditions + beta = currentPatch%fuel%bulk_density_notrunks/SF_val_part_dens + beta_op = OptimumPackingRatio(currentPatch%fuel%SAV_notrunks) + + if (beta_op < nearzero) then + beta_ratio = 0.0_r8 + else + beta_ratio = beta/beta_op + end if + + i_r = ReactionIntensity(currentPatch%fuel%non_trunk_loading/0.45_r8, & + currentPatch%fuel%SAV_notrunks, beta_ratio, & + currentPatch%fuel%average_moisture_notrunks, currentPatch%fuel%MEF_notrunks) + q_ig = HeatofPreignition(currentPatch%fuel%average_moisture_notrunks) + eps = EffectiveHeatingNumber(currentPatch%fuel%SAV_notrunks) + + + ! calculate effective wind speed at CI + ! XLG: this is the disconnection from ROS_active calculation, where they used midflame wind speed + ! as effective wind speed. But it might also make sense as ROS_active is theoretical crown fire + ! rate of spread, and ROS_SA is just the surface fire rate of spread at crowning index wind speed, + ! so it makes senses to consider wind attenuation for grass cover as it's still surface fire?? + + call CalculateTreeGrassAreaSite(currentSite,tree_fraction, grass_fraction, bare_fraction) + CI_effective = CI * (tree_fraction*wind_atten_tree + & + (grass_fraction + bare_fraction)*wind_atten_grass) + ! phi_wind + phi_wind_SA = WindFactor(CI_effective, beta_ratio, & + currentPatch%fuel%SAV_notrunks) + xi = PropagatingFlux(beta, currentPatch%fuel%SAV_notrunks) + ! calculate surface fire spread rate at CI + ROS_SA = ForwardRateOfSpread(currentPatch%fuel%bulk_density_notrunks, & + eps, q_ig, i_r, xi, phi_wind_SA) + + ! Calculate ROS_init, EQ. 12 in Scott & Reinhardt 2001 + ! first calculate heat release per unit area [kW/m2] + HPA = HeatReleasePerArea(currentPatch%fuel%SAV_notrunks, i_r) + ROS_init = (60.0_r8 * FI_init) / HPA + + ! Now check if there is passive or active crown fire and calculate crown fraction burnt (CFB) + ! XLG: there are alternative ways to calculate CFB, see pg 39-41 in Scott & Reinhardt 2001 + call CrownFireCFB(ROS_active, ROS_active_min, currentPatch%ROS_front, & + ROS_init, ROS_SA, active_crownfire, passive_crownfire, crown_frac_burnt) + + currentPatch%active_crown_fire = active_crownfire + currentPatch%passive_crown_fire = passive_crownfire + + ! calculate ROS_final EQ. 21 in Scott & Reinhardt 2001 + ROS_final = currentPatch%ROS_front + crown_frac_burnt * & + ( ROS_active - currentPatch%ROS_front) + + ! calculate final fire intensity by accounting for burned canopy fuels + ! EQ. 22 in Scott & Reinhardt 2001 + FI_final = CrownFireIntensity(HPA, currentPatch%fuel%canopy_fuel_load, & + currentPatch%area, crown_frac_burnt, ROS_final) + + write(fates_log(),*) 'FI_final is ', FI_final + write(fates_log(),*) 'passive crown fire is ', currentPatch%passive_crown_fire + write(fates_log(),*) 'active crown fire is ', currentPatch%active_crown_fire + write(fates_log(),*) 'FI final is ', FI_final + write(fates_log(),*) 'ROS final is ', ROS_final + write(fates_log(),*) 'ROS_active is ', ROS_active + write(fates_log(),*) 'ROS_active_min is ', ROS_active_min + write(fates_log(),*) 'HPA is ', HPA + + + ! only update FI and ROS_front when CFB > 0 + if (crown_frac_burnt > 0.0_r8) then + currentPatch%ROS_front = ROS_final + currentPatch%FI = FI_final + end if + + end if ! end check if there is a crown fire + end if ! if there is a fire + + + currentPatch => currentPatch%younger + end do ! end patch loop + + end subroutine PassiveActiveCrownFireCheck + !--------------------------------------------------------------------------------------- subroutine CalculateAreaBurnt(currentSite) @@ -611,7 +971,14 @@ subroutine CalculatePostFireMortality(currentSite) ! calculate crown fraction burned [0-1] call CrownDepth(currentCohort%height, currentCohort%pft, crown_depth) currentCohort%fraction_crown_burned = CrownFractionBurnt(currentPatch%Scorch_ht(currentCohort%pft), & - currentCohort%height, crown_depth) + currentCohort%height, crown_depth) + + ! If there is an active crown fire, we set CFB to 1 + if(currentPatch%active_crown_fire == 1)then + currentCohort%fraction_crown_burned = 1.0_r8 + end if + + ! shrink canopy to account for burnt section ! currentCohort%canopy_trim = min(currentCohort%canopy_trim, 1.0_r8 - currentCohort%fraction_crown_burned) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 1518f75be5..334d683c92 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -63,10 +63,9 @@ module EDParamsMod real(r8),protected, public :: ED_val_cohort_age_fusion_tol ! minimum fraction in differece in cohort age between cohorts real(r8),protected, public :: ED_val_patch_fusion_tol ! minimum fraction in difference in profiles between patches real(r8),protected, public :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry - - logical,protected, public :: active_crown_fire ! flag, 1=active crown fire 0=no active crown fire - - character(len=param_string_length),parameter :: fates_name_active_crown_fire = "fates_fire_active_crown_fire" + + logical,protected, public :: crown_fire_switch ! flag, 1=active crown fire 0=no active crown fire + character(len=param_string_length),parameter :: fates_name_crown_fire_switch = "fates_crown_fire_switch" real(r8), protected, public :: cg_strikes ! fraction of cloud to ground lightning strikes (0-1) character(len=param_string_length),parameter :: fates_name_cg_strikes="fates_fire_cg_strikes" @@ -545,7 +544,7 @@ 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, & + call fates_params%RegisterParameter(name=fates_name_crown_fire_switch, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) call fates_params%RegisterParameter(name=fates_name_cg_strikes, dimension_shape=dimension_shape_scalar, & @@ -724,9 +723,9 @@ 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, & + call fates_params%RetrieveParameter(name=fates_name_crown_fire_switch, & data=tmpreal) - active_crown_fire = (abs(tmpreal-1.0_r8) fates_r8 + use FatesUnitTestParamReaderMod, only : fates_unit_test_param_reader + use FatesArgumentUtils, only : command_line_arg + use FatesCohortMod, only : fates_cohort_type + use FatesPatchMod, only : fates_patch_type + use FatesFactoryMod, only : InitializeGlobals, GetSyntheticPatch + use SyntheticPatchTypes, only : synthetic_patch_array_type + use FatesTestFireMod, only : SetUpFuel + use SFFireWeatherMod, only : fire_weather + use SFNesterovMod, only : nesterov_index + use SyntheticFuelModels, only : fuel_models_array_class + use FatesFuelClassesMod, only : num_fuel_classes, fuel_classes + use FatesFuelMod, only : fuel_type + use FatesTestCrownFireMod, only : ROSWrapper, EffectiveWindWrapper, WriteCanopyFuelData + use CrownFireEquationsMod, only : MaxHeight , BiomassBin + use CrownFireEquationsMod, only : PassiveCrownFireIntensity, CrownFireBehaveFM10 + use CrownFireEquationsMod, only : HeatReleasePerArea, CrownFireCFB + use CrownFireEquationsMod, only : CrownFireIntensity + use SFEquationsMod, only : ReactionIntensity, FireIntensity + use FatesAllometryMod, only : CrownDepth + use PRTGenericMod, only : leaf_organ + use PRTGenericMod, only : sapw_organ + use PRTGenericMod, only : struct_organ + use PRTGenericMod, only : carbon12_element + use PRTParametersMod, only : prt_params + use FatesConstantsMod, only : itrue + use SFParamsMod, only : SF_val_CWD_frac + use FatesLitterMod, only : adjust_SF_CWD_frac + use FatesLitterMod, only : ncwd + use SFParamsMod, only : SF_val_SAV + use SFParamsMod, only : SF_val_FBD, SF_val_fire_threshold + use SFParamsMod, only : SF_val_part_dens, SF_val_miner_total + + + implicit none + + ! LOCALS: + type(fates_unit_test_param_reader) :: param_reader ! param reader instance + type(synthetic_patch_array_type) :: patch_data ! array of synthetic patches + type(fuel_models_array_class) :: fuel_models_array ! array of fuel models + class(fire_weather), pointer :: fireWeather ! fire weather object + type(fuel_type), allocatable :: fuel(:) ! fuel objects + character(len=:), allocatable :: param_file ! input parameter file + character(len=100), allocatable :: fuel_names(:) ! names of fuel models + character(len=2), allocatable :: carriers(:) ! carriers of fuel models + character(len=100), allocatable :: patch_names(:) ! names of patch types + type(fates_patch_type), pointer :: patch ! patch + type(fates_cohort_type), pointer :: cohort ! cohort + + real(r8) :: effect_wind ! effective wind speed [m/min] + real(r8) :: max_height ! max cohort height [m] + real(r8) :: leaf_c ! leaf carbon [kg C] + real(r8) :: sapw_c ! sap wood carbon [kg C] + real(r8) :: struct_c ! strcutural carbon [kg C] + real(r8) :: woody_c ! aboveground sap wood + structural carbon [kg C] + real(r8) :: canopy_fuel_1h ! leaf + 1-hour woody fuel [kg C] + real(r8), allocatable :: biom_matrix(:) ! array to hold biomass at 1m interval [kg biomass] + real(r8) :: crown_depth ! crown length of a cohort [m] + real(r8) :: cbh_co ! cohort base height, different from path level base height [m] + real(r8) :: cwd_frac_adj(ncwd) ! adjusted fractional allocation of woody biomass to coarse wood debris pool + real(r8) :: ROS ! surface fire forward rate of spread [m/min] + real(r8) :: fuel_consumed(num_fuel_classes) ! fuel consumed by fuel class [kgC/m2] + real(r8) :: tfc_ros ! overall fuel consumed by spreading fire ignoring trunks [kgC/m2] + real(r8) :: ROS_active ! active crown fire ROS using FM10 [m/min] + real(r8) :: CI ! crowning index [m/min] + real(r8) :: CI_cp ! a copy of CI to calculate ROS at CI wind speed + real(r8) :: i_r ! reaction intensity [kJ/m2/min] + real(r8) :: HPA ! heat release per area [kW/m2] + real(r8) :: ROS_init ! ROS for initiating crown fire [m/min] + real(r8) :: crown_frac_burnt ! crown fraction burnt by crown fire [0-1] + real(r8), allocatable :: Wind(:) ! wind speed [m/min] + real(r8), allocatable :: NI(:) ! Nestrove index + real(r8), allocatable :: CWC(:) ! canopy water content [%] + real(r8), allocatable :: FI(:,:,:,:) ! surface fire intensity [kW/m] + real(r8), allocatable :: FI_init(:,:) ! initiation fire intensity by stand type [kW/m] + real(r8), allocatable :: canopy_fuel_load(:) ! patch level canopy fuel load by stand type [kg biomass] + real(r8), allocatable :: CBD(:) ! patch level canopy bulk density [kg biomass / m3] + real(r8), allocatable :: CBH(:) ! patch level canopy base height [m] + real(r8), allocatable :: ROS_front(:,:,:,:) ! surface fire forward rate of spread [m/min] + real(r8), allocatable :: ROS_actfm10(:,:,:,:,:) ! active crown fire ROS [m/min] + real(r8), allocatable :: ROS_actCI(:,:,:,:,:) ! active crown fire ROS [m/min] + real(r8), allocatable :: ROS_critical(:) ! critical ROS for active crown fire to occur, by stand type [m/min] + real(r8), allocatable :: ROS_final(:,:,:,:,:) ! final ROS by stand type and surface fuel model [m/min] + real(r8), allocatable :: FI_final(:,:,:,:,:) ! final fire intensity by stand type and surface fuel model [kW/m] + real(r8), allocatable :: CFB(:,:,:,:,:) ! crown fraction burned by stand type and surface fuel model [fraction] + integer :: p, f, w, n, c, i ! looping indices + integer :: num_fuel_models ! number of fuel models to test + integer :: num_patch_types ! number of patch types to test + integer :: num_wind ! size of wind speed + integer :: active_crownfire ! 1 = active crown fire + integer :: passive_crownfire ! 1 = passive crown fire + + ! CONSTANTS: + integer, parameter :: num_levsoil = 10 ! number of soil layers + real(r8), parameter :: step_size = 1800.0_r8 ! step-size [s] + real(r8), parameter :: biomass_2_carbon = 0.45_r8 ! biomass to carbon multiplier + real(r8), parameter, dimension(4) :: NI_vals = (/300.0_r8, 1000.0_r8, 3000.0_r8, 5000.0_r8/) ! fire weather index to use + real(r8), parameter, dimension(5) :: CWC_vals = (/120.0_r8, 100.0_r8, 90.0_r8, 79.0_r8, 69.0_r8/) ! canopy water content [%] + real(r8), parameter :: wind_max = 800.0_r8 ! max. wind speed to use [m/min] + real(r8), parameter :: wind_min = 100.0_r8 ! min. wind speed to use [m/min] + real(r8), parameter :: wind_inc = 25.0_r8 ! wind speed increment to scale [m/min] + real(r8), parameter :: drying_ratio = 3000.0_r8 + real(r8), parameter :: wind_atten_tree = 0.4_r8 ! wind attenuation factor for tree fraction + real(r8), parameter :: wind_atten_grass = 0.6_r8 ! wind attenuation factor for grass fraction + character(len=*), parameter :: out_file = 'canopyfuel_out.nc' ! output file + + + ! fuel models and patch types to test + integer, parameter, dimension(1) :: fuel_models = (/10/) + integer, parameter, dimension(4) :: patch_ids = (/9, 10, 11, 12/) + + ! number of fuel models, patch types, and wind speed to test + num_fuel_models = size(fuel_models) + num_patch_types = size(patch_ids) + num_wind = int((wind_max - wind_min) / wind_inc + 1) + + ! read in parameter file name from command line + param_file = command_line_arg(1) + + ! allocate arrays + allocate(fuel_names(num_fuel_models)) + allocate(carriers(num_fuel_models)) + allocate(patch_names(num_patch_types)) + allocate(Wind(num_wind)) + allocate(NI(size(NI_vals))) + allocate(CWC(size(CWC_vals))) + allocate(FI(num_wind, size(NI_vals), num_patch_types, num_fuel_models)) + allocate(FI_init(num_patch_types, size(CWC_vals))) + allocate(canopy_fuel_load(num_patch_types)) + allocate(CBD(num_patch_types)) + allocate(CBH(num_patch_types)) + allocate(ROS_front(num_wind, size(NI), num_patch_types, num_fuel_models)) + allocate(ROS_actfm10(num_wind, size(NI), num_patch_types, size(CWC),num_fuel_models)) + allocate(ROS_actCI(num_wind, size(NI), num_patch_types, size(CWC),num_fuel_models)) + allocate(ROS_critical(num_patch_types)) + allocate(ROS_final(num_wind, size(NI), num_patch_types, size(CWC), num_fuel_models)) + allocate(FI_final(num_wind, size(NI), num_patch_types, size(CWC), num_fuel_models)) + allocate(CFB(num_wind, size(NI), num_patch_types, size(CWC), num_fuel_models)) + + ! read in parameter file + call param_reader%Init(param_file) + call param_reader%RetrieveParameters() + + ! set up fire weather class + allocate(nesterov_index :: fireWeather) + call fireWeather%Init() + + ! set up fuel objects and calculate loading + allocate(fuel(num_fuel_models)) + call fuel_models_array%GetFuelModels() + + do f = 1, num_fuel_models + ! uses data from fuel_models to initialize fuel + call SetUpFuel(fuel(f), fuel_models_array, fuel_models(f), fuel_names(f), carriers(f)) + + ! sum up fuel and calculate loading + call fuel(f)%SumLoading() + call fuel(f)%CalculateFractionalLoading() + + ! calculate geometric properties + call fuel(f)%AverageBulkDensity_NoTrunks(SF_val_FBD) + call fuel(f)%AverageSAV_NoTrunks(SF_val_SAV) + end do + + ! initialize some global data we need + call InitializeGlobals(step_size) + + ! get all the patch data + call patch_data%GetSyntheticPatchData() + + ! loop through wind speed, NI and patch type to do all + ! calculation; inside, we also have to loop through fuel + ! models and cwc for certain variables but not all + + do w = 1, num_wind + Wind(w) = wind_min + wind_inc*(w-1) + do n = 1, size(NI_vals) + NI(n) = NI_vals(n) + ! update fire weather index + fireWeather%fire_weather_index = NI(n) + do p = 1, num_patch_types + i = patch_data%PatchDataPosition(patch_id=patch_ids(p)) + call GetSyntheticPatch(patch_data%patches(i), num_levsoil, patch) + + ! update patch tree and grass area + call patch%UpdateTreeGrassArea() + + ! effective wind speed + call EffectiveWindWrapper(patch%total_tree_area, patch%total_grass_area, & + patch%area, Wind(w), effect_wind) + + + ! Calculate ROS for each fuel model + do f = 1, num_fuel_models + ! first update fuel moisture content + call fuel(f)%UpdateFuelMoisture(SF_val_SAV, drying_ratio, fireWeather) + + call ROSWrapper(fuel(f)%bulk_density_notrunks, fuel(f)%SAV_notrunks, & + fuel(f)%non_trunk_loading, fuel(f)%average_moisture_notrunks, & + fuel(f)%MEF_notrunks, NI(n), effect_wind, ROS, i_r) + ROS_front(w, n, p, f) = ROS + + end do + + cohort => patch%tallest + ! initialize max_height + max_height = 0.0_r8 + ! search for max height + do while (associated(cohort)) + if (prt_params%woody(cohort%pft) == itrue) then + call MaxHeight(cohort%height, max_height) + end if + cohort => cohort%shorter + end do + + ! allocate and initialize biom_martix + allocate(biom_matrix(0:int(max_height))) + biom_matrix = 0.0_r8 + + ! derive canopy fuel load + cohort => patch%tallest + do while (associated(cohort)) + if(prt_params%woody(cohort%pft) == itrue) then + call CrownDepth(cohort%height, cohort%pft, crown_depth) + cbh_co = cohort%height - crown_depth + leaf_c = cohort%prt%GetState(leaf_organ, carbon12_element) + sapw_c = cohort%prt%GetState(sapw_organ, carbon12_element) + struct_c = cohort%prt%GetState(struct_organ, carbon12_element) + woody_c = cohort%n * prt_params%allom_agb_frac(cohort%pft)* & + (struct_c + sapw_c) / biomass_2_carbon + leaf_c = cohort%n * leaf_c / biomass_2_carbon + call adjust_SF_CWD_frac(cohort%dbh, ncwd, SF_val_CWD_frac, cwd_frac_adj) + + ! update canopy fuel load, this has nothing to do with fuel model + ! but we save it for each fuel model for later use + do f = 1, num_fuel_models + call fuel(f)%CalculateCanopyFuelLoad(leaf_c, woody_c, & + cwd_frac_adj, canopy_fuel_1h) + end do + + cohort%canopy_fuel_1h = canopy_fuel_1h + + ! 1m biomass bin + call BiomassBin(cbh_co, cohort%height, crown_depth, canopy_fuel_1h, biom_matrix) + end if + cohort => cohort%shorter + end do + + ! calculate canopy bulk density + biom_matrix(:) = biom_matrix(:) / patch%area ! kg biomass / m3 + + do f = 1, num_fuel_models + call fuel(f)%CalculateCanopyBulkDensity(biom_matrix, max_height) + end do + + ! save canopy fuel characteristics + canopy_fuel_load(p) = fuel(1)%canopy_fuel_load + CBD(p) = fuel(1)%canopy_bulk_density + CBH(p) = fuel(1)%canopy_base_height + + deallocate(biom_matrix) + + ! calculate critical ROS + ROS_critical(p) = 3.0_r8 / CBD(p) + + ! loop through CWC and fuel models to determine + ! crown fire occurence and calculate relevant variales + + do c = 1, size(CWC_vals) + CWC(c) = CWC_vals(c) + ! calculate initiation surface fire intensity + FI_init(p,c) = PassiveCrownFireIntensity(CBH(p), CWC(c)) + + do f = 1, num_fuel_models + call fuel(f)%CalculateFuelBurnt(fuel_consumed) + tfc_ros = sum(fuel_consumed) - fuel_consumed(fuel_classes%trunks()) + FI(w,n,p,f) = FireIntensity(tfc_ros/0.45_r8, ROS_front(w,n,p,f)/60.0_r8) + + if(FI(w,n,p,f) > SF_val_fire_threshold .and. & + FI(w,n,p,f) > FI_init(p,c))then + call CrownFireBehaveFM10(drying_ratio, NI(n), & + SF_val_miner_total, SF_val_part_dens, Wind(w), & + CBD(p), ROS_active, CI) + ROS_actfm10(w,n,p,c,f) = ROS_active + CI_cp = CI + call CrownFireBehaveFM10(drying_ratio, NI(n), & + SF_val_miner_total, SF_val_part_dens, CI_cp, & + CBD(p), ROS_active, CI) + ROS_actCI(w,n,p,c,f) = ROS_active + + ! calculate effective wind speed using CI + call EffectiveWindWrapper(patch%total_tree_area, patch%total_grass_area, & + patch%area, CI, effect_wind) + + ! calculate ROS_SA, which is the returned ROS + call ROSWrapper(fuel(f)%bulk_density_notrunks, fuel(f)%SAV_notrunks, & + fuel(f)%non_trunk_loading, fuel(f)%average_moisture_notrunks, & + fuel(f)%MEF_notrunks, NI(n), effect_wind, ROS, i_r) + + HPA = HeatReleasePerArea(fuel(f)%SAV_notrunks, i_r) + ROS_init = (60.0_r8 * FI_init(p,c)) / HPA + + ! calculate crown fraction burnt + call CrownFireCFB(ROS_active, ROS_critical(p), ROS_front(w,n,p,f), & + ROS_init, ROS, active_crownfire, passive_crownfire, crown_frac_burnt) + CFB(w,n,p,c,f) = crown_frac_burnt + else + ROS_actfm10(w,n,p,c,f) = 0.0_r8 + crown_frac_burnt = 0.0_r8 + CFB(w,n,p,c,f) = 0.0_r8 + end if + + ! update ROS and FI + if(crown_frac_burnt > 0.0_r8)then + ROS_final(w,n,p,c,f) = ROS_front(w,n,p,f) + crown_frac_burnt* & + (ROS_actfm10(w,n,p,c,f) - ROS_front(w,n,p,f)) + FI_final(w,n,p,c,f) = CrownFireIntensity(HPA, fuel(f)%canopy_fuel_load, & + patch%area, crown_frac_burnt, ROS_final(w,n,p,c,f)) + else + ROS_final(w,n,p,c,f) = ROS_front(w,n,p,f) + FI_final(w,n,p,c,f) = FI(w,n,p,f) + end if + + end do + + end do + end do + end do + end do + + ! write out data + call WriteCanopyFuelData(out_file, num_fuel_models, num_patch_types, num_wind, Wind, & + NI, CWC, CBD, CBH, canopy_fuel_load, ROS_front, FI, FI_init, ROS_actfm10, ROS_critical, & + ROS_actCI, CFB, ROS_final, FI_final, fuel_models, patch_ids) + + +end program FatesTestCanopyFuel diff --git a/testing/functional_testing/fire/canopyfuel/canopyfuel_test.py b/testing/functional_testing/fire/canopyfuel/canopyfuel_test.py new file mode 100644 index 0000000000..7fe58f2bbc --- /dev/null +++ b/testing/functional_testing/fire/canopyfuel/canopyfuel_test.py @@ -0,0 +1,378 @@ + +""" +Concrete class for running the canopyfuel functional test for FATES. +""" +import os +import numpy as np +import xarray as xr +import pandas as pd +import matplotlib.pyplot as plt +from functional_class import FunctionalTest +from utils_plotting import get_color_palette + +class CanopyFuelTest(FunctionalTest): + """Canopy fuel test class""" + + name = "canopyfuel" + + def __init__(self, test_dict): + super().__init__( + CanopyFuelTest.name, + test_dict["test_dir"], + test_dict["test_exe"], + test_dict["out_file"], + test_dict["use_param_file"], + test_dict["other_args"], + ) + self.plot = True + + def plot_output(self, run_dir: str, save_figs: bool, plot_dir: str): + """Plot output associated with canopy fuel tests + + Args: + run_dir (str): run directory + out_file (str): output file + save_figs (bool): whether or not to save the figures + plot_dir (str): plot directory + """ + + # read in canopy fuel data + cfuel_dat = xr.open_dataset(os.path.join(run_dir, self.out_file)) + + ## 1. some barchart plots for variables with one dimension + + self.plot_barchart( + cfuel_dat, + "CBD", + "Canopy bulk density", + "kg m$^{-3}$", + save_figs, + plot_dir, + by_fuel_model=False, + stacked= False, + ) + self.plot_barchart( + cfuel_dat, + "CBH", + "Canopy base height", + "m", + save_figs, + plot_dir, + by_fuel_model=False, + stacked= False, + ) + self.plot_barchart( + cfuel_dat, + "ROS_min", + "Critical rate of spread", + "m min$^{-1}$", + save_figs, + plot_dir, + by_fuel_model=False, + stacked= False, + ) + + ## 2. plots for variables with at least 2 dimensions + self.plot_facet( + cfuel_dat, + x="canopy_water", y="FI_init", + hue="patch_type",row=None,col=None,style=None, + kind="line", agg=None, + x_label="Canopy water content(%)", + y_label="Fire intensity to ignite crown fuel(kW m-1)", + save_figs=save_figs, + plot_dir=plot_dir, + filename="FI_init_plot.png", + ) + self.plot_facet( + cfuel_dat, + x="wind", y="ROS_front", + hue="fire_weather",row="fuel_model",col="patch_type",style=None, + kind="line", agg=None, + x_label="Wind speed(m min-1)", + y_label="Surface fire rate of spread(m min-1)", + save_figs=save_figs, + plot_dir=plot_dir, + filename="ROS_front_plot.png", + ) + self.plot_facet( + cfuel_dat, + x="wind", y="ROS_active", + hue="fire_weather",row="canopy_water",col="patch_type",style="fuel_model", + kind="line", agg=None, + x_label="Wind speed(m min-1)", + y_label="Active crown fire rate of spread(m min-1)", + save_figs=save_figs, + plot_dir=plot_dir, + filename="ROS_act_plot.png", + ) + self.plot_facet( + cfuel_dat, + x="wind",y="ROS_final", + hue="fire_weather",row="canopy_water",col="patch_type",style="fuel_model", + kind="line", agg=None, + x_label="Wind speed(m min-1)", + y_label="Final rate of spread(m min-1)", + save_figs=save_figs, + plot_dir=plot_dir, + filename="ROS_final_plot.png", + ) + self.plot_facet( + cfuel_dat, + x="wind", y="FI_final", + hue="fire_weather",row="canopy_water",col="patch_type",style="fuel_model", + kind="line", agg=None, + x_label="Wind speed(m min-1)", + y_label="Final fire intensity(kW m-1)", + save_figs=save_figs, + plot_dir=plot_dir, + filename="FI_final_plot.png", + ) + self.plot_facet( + cfuel_dat, + x="wind",y="CFB", + hue="fire_weather",row="canopy_water",col="patch_type",style="fuel_model", + kind="line", agg=None, + x_label="Wind speed(m min-1)", + y_label="Crown fraction burnt", + save_figs=save_figs, + plot_dir=plot_dir, + filename="CFB_plot.png", + + ) + +# ------- Define plot functions ------- + + @staticmethod + def plot_barchart( + cfuel_dat: xr.Dataset, + var: str, + varname: str, + units: str, + save_figs: bool, + plot_dir: str, + by_fuel_model: bool = False, + stacked: bool = False, + ): + """ Plot canopy fuel test outputs as bar plot + Args: + cfuel_dat (xr.Dataset): canopy fuel test data output + var (str): variable to plot + varname (str): variable name for y-axis + units (str): unit expression + save_figs (str): dir to save figure + plot_dir (str): which dir to save figs + by_fuel_model (bool, optional): whether or not the bar plot is grouped by fuel model, default to True + + """ + assert var in cfuel_dat, f"{var!r} not found in dataset variables." + da = cfuel_dat[var] + + # figure out the dimensions we care about + has_patch = "patch_type" in da.dims + has_fm = "fuel_model" in da.dims + + if not has_patch: + raise ValueError(f"{var} is missing 'patch_type' dimension; got dims {da.dims}") + + patch_types = [str(p) for p in cfuel_dat["patch_type"].values] + + + if by_fuel_model and has_fm: + da2 = da.transpose("patch_type", "fuel_model") + arr = da2.to_numpy() + fm_labels = [str(f) for f in cfuel_dat["fuel_model"].values] + + fig, ax = plt.subplots(figsize=(8, 4.6)) + x = np.arange(len(patch_types)) + n_fm = arr.shape[1] + width = 0.8 / n_fm + + if stacked: + bottom = np.zeros(len(patch_types)) + for j in range(n_fm): + ax.bar( + x, + arr[:, j], + width=0.8, + bottom=bottom, + label=fm_labels[j], + ) + bottom += arr[:, j] + else: + for j in range(n_fm): + offset = (j - (n_fm - 1) / 2) * width + ax.bar( + x + offset, + arr[:, j], + width=width, + label=fm_labels[j], + ) + ax.set_xticks(x) + ax.set_xticklabels(patch_types, rotation=90) + ax.legend(title="fuel_model", loc="center left", bbox_to_anchor=(1, 0.5)) + else: + if has_fm: + da_plot = da.mean(dim="fuel_model") + else: + da_plot = da + + da_plot = da_plot.transpose("patch_type", ...) + y = da_plot.to_numpy() + fig, ax = plt.subplots(figsize=(7.5, 4.2)) + ax.bar(patch_types, y) + ax.set_xticklabels(patch_types, rotation=90) + + ax.set_ylabel(f"{varname} ({units})", fontsize=11) + ax.set_xlabel("Patch Type") + ax.set_title(varname) + plt.tight_layout() + + if save_figs: + fig_name = os.path.join(plot_dir, f"{varname}_plot.png") + plt.savefig(fig_name) + + @staticmethod + def plot_facet( + ds, *, + x: str, + y: str, + hue: str | None = None, + row: str | None = None, + col: str | None = None, + style: str | None = None, + kind: str = "line", # "line" or "scatter" + agg: str | None = "mean",# "mean", "median", or None (no aggregation) + orders: dict | None = None, # e.g., {"CWC":[40,50,60,70,80], "NI":[1000,2000,...]} + x_label: str | None = None, + y_label: str | None = None, + title: str | None = None, + save_figs: bool = False, + plot_dir: str | None = None, + filename: str = "facet_plot.png", + ): + """ + General faceted plot using matplotlib. + + Panels: row × col + Inside each panel: one series per `hue`; optionally split each hue by `style`. + + - If `agg` is provided, aggregates y vs x within each group (good for many points). + - If `kind="scatter"`, plots points; if "line", draws lines (with optional markers). + """ + + # Tidy DataFrame + if hasattr(ds, "to_dataframe"): + df = ds.to_dataframe().reset_index() + else: + df = ds.copy() + + # Collect the fields actually used + fields = [x, y] + [v for v in [hue, row, col, style] if v] + missing = [k for k in fields if k not in df.columns] + if missing: + raise KeyError(f"Missing columns: {missing}") + + df = df.dropna(subset=[x, y]) # basic hygiene + + # Helpers for level orders + orders = orders or {} + def _levels(key): + return orders.get(key, sorted(pd.unique(df[key]))) if key else [None] + + row_levels = _levels(row) + col_levels = _levels(col) + hue_levels = _levels(hue) + + n_rows, n_cols = len(row_levels), len(col_levels) + fig, axes = plt.subplots( + n_rows, n_cols, + figsize=(3.8 * n_cols, 2.8 * n_rows), + sharex=True, sharey=True + ) + if n_rows == 1 and n_cols == 1: + axes = np.array([[axes]]) + elif n_rows == 1: + axes = np.array([axes]) + elif n_cols == 1: + axes = axes[:, None] + + # Aggregator + if agg is None: + def _agg(g): + return g.sort_values(x)[[x, y]] + else: + agg_name = {"mean": "mean", "median": "median"}[agg] + def _agg(g): + out = g.groupby(x, as_index=False)[y].agg(agg_name) + return out.sort_values(x) + + # Plot panels + for r, rv in enumerate(row_levels): + for c, cv in enumerate(col_levels): + ax = axes[r, c] + sub = df.copy() + if row: sub = sub[sub[row] == rv] + if col: sub = sub[sub[col] == cv] + + ttl_bits = [] + if row: ttl_bits.append(f"{row}={rv}") + if col: ttl_bits.append(f"{col}={cv}") + ax.set_title(" • ".join(ttl_bits) if ttl_bits else (title or "")) + + if sub.empty: + ax.text(0.5, 0.5, "no data", ha="center", va="center", transform=ax.transAxes) + continue + + # Loop hues + for hv in hue_levels: + sub_h = sub if not hue else sub[sub[hue] == hv] + if sub_h.empty: + continue + + # Optional style split + if style: + for sv, sub_s in sub_h.groupby(style): + cur = _agg(sub_s) + if cur.empty: + continue + if kind == "scatter": + ax.scatter(cur[x].to_numpy(), cur[y].to_numpy(), + label=f"{hue}={hv}, {style}={sv}" if hue else f"{style}={sv}", s=18) + else: + ax.plot(cur[x].to_numpy(), cur[y].to_numpy(), + marker="o", linestyle="-", linewidth=1.25, markersize=3.5, + label=f"{hue}={hv}, {style}={sv}" if hue else f"{style}={sv}") + else: + cur = _agg(sub_h) + if cur.empty: + continue + if kind == "scatter": + ax.scatter(cur[x].to_numpy(), cur[y].to_numpy(), + label=f"{hue}={hv}" if hue else None, s=18) + else: + ax.plot(cur[x].to_numpy(), cur[y].to_numpy(), + marker="o", linestyle="-", linewidth=1.25, markersize=3.5, + label=f"{hue}={hv}" if hue else None) + + if r == n_rows - 1: + ax.set_xlabel(x_label or x) + if c == 0: + ax.set_ylabel(y_label or y) + + # Consolidated legend + handles, labels = axes[0, 0].get_legend_handles_labels() + if handles: + fig.legend(handles, labels, loc="center right", bbox_to_anchor=(1.02, 0.5), + title=(" × ".join([k for k in [hue, style] if k])) or None) + + if title and (n_rows * n_cols > 1): + fig.suptitle(title, y=1.02) + + fig.tight_layout(rect=(0, 0, 0.9, 1)) + if save_figs and plot_dir: + os.makedirs(plot_dir, exist_ok=True) + fig.savefig(os.path.join(plot_dir, filename), dpi=150, bbox_inches="tight") + plt.close(fig) + + diff --git a/testing/functional_testing/fire/crownfire/CMakeLists.txt b/testing/functional_testing/fire/crownfire/CMakeLists.txt new file mode 100644 index 0000000000..f3be016385 --- /dev/null +++ b/testing/functional_testing/fire/crownfire/CMakeLists.txt @@ -0,0 +1,25 @@ +set(crownfire_test_sources + FatesTestCrownFire.F90) + +set(NETCDF_C_DIR ${NETCDF_C_PATH}) +set(NETCDF_FORTRAN_DIR ${NETCDF_F_PATH}) + +FIND_PATH(NETCDFC_FOUND libnetcdf.a ${NETCDF_C_DIR}/lib) +FIND_PATH(NETCDFF_FOUND libnetcdff.a ${NETCDF_FORTRAN_DIR}/lib) + + +include_directories(${NETCDF_C_DIR}/include + ${NETCDF_FORTRAN_DIR}/include) + +link_directories(${NETCDF_C_DIR}/lib + ${NETCDF_FORTRAN_DIR}/lib + ${PFUNIT_TOP_DIR}/lib) + +add_executable(FATES_crownfire_exe ${crownfire_test_sources}) + +target_link_libraries(FATES_crownfire_exe + netcdf + netcdff + fates + csm_share + funit) \ No newline at end of file diff --git a/testing/functional_testing/fire/crownfire/FatesTestCrownFire.F90 b/testing/functional_testing/fire/crownfire/FatesTestCrownFire.F90 new file mode 100644 index 0000000000..59eb3331e8 --- /dev/null +++ b/testing/functional_testing/fire/crownfire/FatesTestCrownFire.F90 @@ -0,0 +1,442 @@ +program FatesTestCrownFire + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesArgumentUtils, only : command_line_arg + use FatesUnitTestParamReaderMod, only : fates_unit_test_param_reader + + implicit none + + ! LOCALS: + type(fates_unit_test_param_reader) :: param_reader ! param reader instance + character(len=:), allocatable :: param_file ! input parameter file + real(r8), allocatable :: CWC(:) ! canopy water content [%] + real(r8), allocatable :: CBD(:) ! canopy bulk density in biomass [kg/m3] + real(r8), allocatable :: CBH(:) ! canopy base height [m] + real(r8), allocatable :: smp(:) ! soil matric potential [MPa] + real(r8), allocatable :: smp_alpha(:) ! coefficient associate with smp for LFMC [unitless] + real(r8), allocatable :: wind(:) ! wind speed [km/hr] + real(r8), allocatable :: drying_ratio(:) ! drying ratio [unitless] + real(r8), allocatable :: passive_crown_fi(:,:) ! min surface fire intensity to ignite crown fuel [kW/m] + real(r8), allocatable :: CI_FM10(:,:) ! open wind speed at which a fully active crown fire is maintained using fule model 10 [km/hr] + real(r8), allocatable :: LFMC(:,:) ! live fuel moisture content [%] + real(r8), allocatable :: ROS_active_FM10(:,:) ! active crown fire spread rate using fuel model 10 [m/min] + + ! CONSTANTS: + character(len=*), parameter :: out_file = 'crownfire_out.nc' ! output file + + interface + + subroutine TestPassCrownFI(CBH, CWC, passive_crown_fi) + + use FatesConstantsMod, only : r8 => fates_r8 + use CrownFireEquationsMod, only : PassiveCrownFireIntensity + implicit none + real(r8), allocatable, intent(out) :: CBH(:) + real(r8), allocatable, intent(out) :: CWC(:) + real(r8), allocatable, intent(out) :: passive_crown_fi(:,:) + + end subroutine TestPassCrownFI + + subroutine TestLiveFuelMoisture(smp, smp_alpha, LFMC) + + use FatesConstantsMod, only : r8 => fates_r8 + use CrownFireEquationsMod, only : LiveFuelMoistureContent + implicit none + real(r8), allocatable, intent(out) :: smp(:) + real(r8), allocatable, intent(out) :: smp_alpha(:) + real(r8), allocatable, intent(out) :: LFMC(:,:) + + end subroutine TestLiveFuelMoisture + + subroutine TestCrownFireFM10(CBD, wind, drying_ratio, ROS_active_FM10, CI_FM10) + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : nearzero + use CrownFireEquationsMod, only : CrownFireBehaveFM10 + implicit none + real(r8), allocatable, intent(out) :: CBD(:) + real(r8), allocatable, intent(out) :: wind(:) + real(r8), allocatable, intent(out) :: drying_ratio(:) + real(r8), allocatable, intent(out) :: ROS_active_FM10(:,:) + real(r8), allocatable, intent(out) :: CI_FM10(:,:) + + end subroutine TestCrownFireFM10 + + subroutine WriteCrownFireData(out_file, CBH, CWC, passive_crown_fi, & + smp, smp_alpha, LFMC, CBD, wind, drying_ratio, ROS_active_FM10, CI_FM10) + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesUnitTestIOMod, only : OpenNCFile, CloseNCFile, RegisterNCDims + use FatesUnitTestIOMod, only : RegisterVar, EndNCDef, WriteVar + use FatesUnitTestIOMod, only : type_double + implicit none + character(len=*), intent(in) :: out_file + real(r8), intent(in) :: CBH(:) + real(r8), intent(in) :: CWC(:) + real(r8), intent(in) :: passive_crown_fi(:,:) + real(r8), intent(in) :: smp(:) + real(r8), intent(in) :: smp_alpha(:) + real(r8), intent(in) :: LFMC(:,:) + real(r8), intent(in) :: CBD(:) + real(r8), intent(in) :: wind(:) + real(r8), intent(in) :: drying_ratio(:) + real(r8), intent(in) :: ROS_active_FM10(:,:) + real(r8), intent(in) :: CI_FM10(:,:) + + end subroutine WriteCrownFireData + + end interface + + ! read in parameter file name + param_file = command_line_arg(1) + + ! read in parameter file + call param_reader%Init(param_file) + call param_reader%RetrieveParameters() + + ! calculate minimum fire intensity for igniting crown fuels + call TestPassCrownFI(CBH, CWC, passive_crown_fi) + + ! calculate live fuel moisture content + call TestLiveFuelMoisture(smp, smp_alpha, LFMC) + + ! calculate active crown fire spread rate and crowning index using fuel model 10 + call TestCrownFireFM10(CBD, wind, drying_ratio, ROS_active_FM10, CI_FM10) + + + ! write out data + call WriteCrownFireData(out_file, CBH, CWC, passive_crown_fi, & + smp, smp_alpha, LFMC, CBD, wind, drying_ratio, ROS_active_FM10, CI_FM10) + + ! deallocate arrays + if (allocated(CBH)) deallocate(CBH) + if(allocated(CWC)) deallocate(CWC) + if(allocated(passive_crown_fi)) deallocate(passive_crown_fi) + if(allocated(smp)) deallocate(smp) + if(allocated(smp_alpha)) deallocate(smp_alpha) + if(allocated(LFMC)) deallocate(LFMC) + if(allocated(CBD)) deallocate(CBD) + if(allocated(wind)) deallocate(wind) + if(allocated(drying_ratio)) deallocate(drying_ratio) + if(allocated(ROS_active_FM10)) deallocate(ROS_active_FM10) + if(allocated(CI_FM10)) deallocate(CI_FM10) + +end program FatesTestCrownFire + +!========================================================================================= + +subroutine TestPassCrownFI(CBH, CWC, passive_crown_fi) + ! + ! DESCRIPTION: + ! Calculate min surface fire intensity to ignite crown fuel over a range of + ! canopy base height and canopy water content + ! + use FatesConstantsMod, only : r8 => fates_r8 + use CrownFireEquationsMod, only : PassiveCrownFireIntensity + + implicit none + + ! ARGUMENTS: + real(r8), allocatable, intent(out) :: CBH(:) ! canopy base height [m] + real(r8), allocatable, intent(out) :: CWC(:) ! canopy water content [%] + real(r8), allocatable, intent(out) :: passive_crown_fi(:,:) ! min energy to ignite crown fuels [kW/m] + + ! CONSTANTS: + real(r8), parameter :: CBH_min = 0.0_r8 ! min canopy base height [m] + real(r8), parameter :: CBH_max = 10.0_r8 ! max canopy base height [m] + real(r8), parameter :: CBH_inc = 1.0_r8 ! CBH increment to scale [m] + real(r8), parameter, dimension(5) :: canopy_water_content = (/10.0_r8, 30.0_r8, 75.0_r8, 85.0_r8, 120.0_r8/) ! CWC to use + + !LOCALS: + integer :: num_CBH ! size of CBH arrays + integer :: i, j ! looping indices + + ! allocate arrays + num_CBH = int((CBH_max - CBH_min) / CBH_inc + 1) + allocate(passive_crown_fi(num_CBH, size(canopy_water_content))) + allocate(CBH(num_CBH)) + allocate(CWC(size(canopy_water_content))) + + do i = 1, num_CBH + + CBH(i) = CBH_min + CBH_inc*(i-1) + + do j = 1, size(canopy_water_content) + CWC(j) = canopy_water_content(j) + passive_crown_fi(i, j) = PassiveCrownFireIntensity(CBH(i), CWC(j)) + end do + end do + +end subroutine TestPassCrownFI + +!========================================================================================= + +subroutine TestLiveFuelMoisture(smp, smp_alpha, LFMC) + ! + ! DESCRIPTION: + ! Calculate live fuel moisture content over a range of soil matirc potential and associated coefficient + ! + + use FatesConstantsMod, only : r8 => fates_r8 + use CrownFireEquationsMod, only : LiveFuelMoistureContent + + implicit none + + ! ARGUMENTS: + real(r8), allocatable, intent(out) :: smp(:) ! soil matric potential [MPa] + real(r8), allocatable, intent(out) :: smp_alpha(:) ! coefficient associate with smp for predicting LFMC [unitless] + real(r8), allocatable, intent(out) :: LFMC(:,:) ! live fule moisture content [%] + + ! CONSTANTS: + real(r8), parameter :: smp_max = 0.0_r8 ! max soil matric potential [MPa] + real(r8), parameter :: smp_min = -10.0_r8 ! min soil matric potential [MPa] + real(r8), parameter :: smp_inc = -1.0_r8 ! smp increment to scale [MPa] + real(r8), parameter, dimension(5) :: smp_coef = (/0.1_r8, 0.09_r8, 0.2_r8, 0.5_r8, 0.05_r8/) ! model coeff associated with smp to use + real(r8), parameter :: lai = 3.0_r8 ! leaf area index for LFMC [m2/m2] + real(r8), parameter :: min_lfmc = 70.0_r8 ! min LFMC [%] + real(r8), parameter :: coef_lfmc = 50.0_r8 ! value add to min LFMC as a response to change in smp and lai [unitless] + real(r8), parameter :: lai_beta = 0.15_r8 ! coefficient associate with lai for LFMC [unitless] + real(r8), parameter :: gamma_int = 0.0_r8 ! coefficient associate with the LAI and SMP interaction term + + ! LOCALS: + integer :: num_smp ! size of soil matric potential + integer :: i, j ! looping indicies + + ! allocate arrays + num_smp = int((smp_max - smp_min) / abs(smp_inc) + 1) + allocate(smp(num_smp)) + allocate(smp_alpha(size(smp_coef))) + allocate(LFMC(num_smp, size(smp_coef))) + + do i = 1, num_smp + + smp(i) = smp_max + smp_inc*(i-1) + + do j = 1, size(smp_coef) + smp_alpha(j) = smp_coef(j) + LFMC(i, j) = LiveFuelMoistureContent(lai, smp(i), min_lfmc, coef_lfmc, & + smp_alpha(j), lai_beta, gamma_int) + end do + end do + +end subroutine TestLiveFuelMoisture + +!========================================================================================= + +subroutine TestCrownFireFM10(CBD, wind, drying_ratio, ROS_active_FM10, CI_FM10) + ! + ! DESCRIPTION: + ! Calculate fully active crown fire spread rate using fule model 10 over a range of CBD and wind speed + ! and open wind speed to initiate crown fire over a range of CBD assuming fule model 10 + ! + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : nearzero + use SFParamsMod, only : SF_val_miner_total + use SFParamsMod, only : SF_val_part_dens + use CrownFireEquationsMod, only : CrownFireBehaveFM10 + + + implicit none + + ! ARGUMENTS: + real(r8), allocatable, intent(out) :: CBD(:) ! canopy bulk density in biomass [kg/m3] + real(r8), allocatable, intent(out) :: wind(:) ! wind speed [km/hr] + real(r8), allocatable, intent(out) :: drying_ratio(:) ! drying ratio that controls how fast fuel dries out [unitless] + real(r8), allocatable, intent(out) :: ROS_active_FM10(:,:) ! fully active crown fire spread rate using fule model 10 + real(r8), allocatable, intent(out) :: CI_FM10(:,:) ! crowning index using fule model 10 [km/hr] + + ! CONSTANTS: + real(r8), parameter :: CBD_min = 0.0_r8 ! min canopy bulk density [kg/m3] + real(r8), parameter :: CBD_max = 0.3_r8 ! max canopy bulk density [kg/m3] + real(r8), parameter :: CBD_inc = 0.01_r8 ! CBD increment to scale [unitless] + real(r8), parameter :: wind_min = 0.0_r8 ! min open wind speed [km/hr] + real(r8), parameter :: wind_max = 50.0_r8 ! max open wind speed [km/hr] + real(r8), parameter :: wind_inc = 1.0_r8 ! wind increment to scale [km/hr] + real(r8), parameter, dimension(5) :: dratio_vals = (/1000.0_r8, 3000.0_r8, 8000.0_r8, 15000.0_r8, 45000.0_r8/) ! drying ratio values to use + real(r8), parameter :: cbd_ref = 0.2_r8 ! reference canopy bulk density [kg/m3] + real(r8), parameter :: wind_ref = 350.0_r8 ! reference wind speed [m/min] + real(r8), parameter :: fire_weather_index = 5000.0_r8 ! Nesterove fire weather index [unitless] + real(r8), parameter :: kmhr_to_mmin = 16.6667_r8 ! convert km/hour to m/min for wind speed + + ! LOCALS: + integer :: num_CBD ! size of canopy bulk density array + integer :: num_wind ! size of wind speed array + integer :: i, j ! looping indicies + real(r8) :: ROS_active ! ROS_active output from CrownFireBehaveFM10 [m/min] + real(r8) :: CI ! CI output from CrownFireBehaveFM10 [m/min] + + ! allocate arrays + num_CBD = int((CBD_max - CBD_min) / CBD_inc + 1) + num_wind = int((wind_max - wind_min) / wind_inc + 1) + allocate(CBD(num_CBD)) + allocate(wind(num_wind)) + allocate(drying_ratio(size(dratio_vals))) + allocate(ROS_active_FM10(num_wind, size(dratio_vals))) + allocate(CI_FM10(num_CBD, size(dratio_vals))) + + do i = 1, num_CBD + CBD(i) = CBD_min + CBD_inc * (i-1) + + do j = 1, size(dratio_vals) + drying_ratio(j) = dratio_vals(j) + call CrownFireBehaveFM10(drying_ratio(j), fire_weather_index, SF_val_miner_total, & + SF_val_part_dens, wind_ref, CBD(i), ROS_active, CI) + CI_FM10(i, j) = CI + + end do + end do + + do i = 1, num_wind + wind(i) = (wind_min + wind_inc * (i-1)) * kmhr_to_mmin ! convert wind speed to m/min + + do j = 1, size(dratio_vals) + drying_ratio(j) = dratio_vals(j) + call CrownFireBehaveFM10(drying_ratio(j), fire_weather_index, SF_val_miner_total, & + SF_val_part_dens, wind(i), cbd_ref, ROS_active, CI) + ROS_active_FM10(i, j) = ROS_active + + end do + end do + +end subroutine TestCrownFireFM10 + +!========================================================================================= + +subroutine WriteCrownFireData(out_file, CBH, CWC, passive_crown_fi, & + smp, smp_alpha, LFMC, CBD, wind, drying_ratio, ROS_active_FM10, CI_FM10) + ! + ! DESCRIPTION: + ! write out data from the test + ! + use FatesConstantsMod, only : r8 => fates_r8 + use FatesUnitTestIOMod, only : OpenNCFile, CloseNCFile, RegisterNCDims + use FatesUnitTestIOMod, only : RegisterVar, EndNCDef, WriteVar + use FatesUnitTestIOMod, only : type_double + + implicit none + + ! ARGUMENTS: + character(len=*), intent(in) :: out_file + real(r8), intent(in) :: CBH(:) + real(r8), intent(in) :: CWC(:) + real(r8), intent(in) :: passive_crown_fi(:,:) + real(r8), intent(in) :: smp(:) + real(r8), intent(in) :: smp_alpha(:) + real(r8), intent(in) :: LFMC(:,:) + real(r8), intent(in) :: CBD(:) + real(r8), intent(in) :: wind(:) + real(r8), intent(in) :: drying_ratio(:) + real(r8), intent(in) :: ROS_active_FM10(:,:) + real(r8), intent(in) :: CI_FM10(:,:) + + ! LOCALS: + integer :: ncid ! netcdf id + character(len=20) :: dim_names(7) ! dimension names + integer :: dimIDs(7) ! dimension IDs + integer :: CBHID + integer :: CWCID + integer :: pasv_crwn_ID + integer :: smpID + integer :: smp_alpha_ID + integer :: LFMCID + integer :: CBDID + integer :: windID + integer :: dratioID + integer :: ROSACT_FM10_ID + integer :: CI_FM10_ID + + ! dimension names + dim_names = [character(len=20) :: 'CBH', 'CWC', 'smp', 'smp_alpha', 'CBD', 'wind', 'dratio'] + + ! open file + call OpenNCFile(trim(out_file), ncid, 'readwrite') + + ! register dimensions + call RegisterNCDims(ncid, dim_names, (/size(CBH), size(CWC), size(smp), & + size(smp_alpha), size(CBD), size(wind), size(drying_ratio)/), 7, dimIDs) + + ! first register dimension variables + ! register CBH + call RegisterVar(ncid, 'CBH', dimIDs(1:1), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: 'm', 'canopy base height'], 2, CBHID) + + ! register canopy water content + call RegisterVar(ncid, 'CWC', dimIDs(2:2), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: '%', 'canopy water content'], 2, CWCID) + + ! register soil matric potential + call RegisterVar(ncid, 'smp', dimIDs(3:3), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: 'MPa', 'soil matric potential'], 2, smpID) + + ! register LFMC model coeff associated with smp + call RegisterVar(ncid, 'smp_alpha', dimIDs(4:4), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: '', 'live fuel moisture model coefficient for soil matric potential'], 2, smp_alpha_ID) + + ! register canopy bulk density + call RegisterVar(ncid, 'CBD', dimIDs(5:5), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: 'kg/m3', 'canopy bulk density in biomass'], 2, CBDID) + + ! register wind speed + call RegisterVar(ncid, 'wind', dimIDs(6:6), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: 'm/min', 'open wind speed'], 2, windID) + + ! register drying ratio + call RegisterVar(ncid, 'dratio', dimIDs(7:7), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: '', 'drying ratio'], 2, dratioID) + + ! then register actual variables + + ! register min surface fire intensity for igniting crown fuels + call RegisterVar(ncid, 'passive_crown_fi', dimIDs(1:2), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'CBH CWC', 'kW/m', 'min fire intensity to ignite crown fuels'], & + 3, pasv_crwn_ID) + + ! register live fuel moisture content + call RegisterVar(ncid, 'LFMC', dimIDs(3:4), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'smp smp_alpha', '%', 'live fuel moisture content'], & + 3, LFMCID) + + ! register active crown fire spread rate using fuel model 10 + call RegisterVar(ncid, 'ROSACT_FM10', dimIDs(6:7), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind dratio', 'm/min', 'active crown fire ROS using fuel model 10'], & + 3, ROSACT_FM10_ID) + + ! register crowning index + call RegisterVar(ncid, 'CI_FM10', (/dimIDs(5),dimIDs(7)/), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'CBD dratio', 'm/min', 'wind speed at which a active crown fire is sustained'], & + 3, CI_FM10_ID) + + + ! finish defining variables + call EndNCDef(ncid) + + ! write out data + call WriteVar(ncid, CBHID, CBH(:)) + call WriteVar(ncid, CWCID, CWC(:)) + call WriteVar(ncid, pasv_crwn_ID, passive_crown_fi(:,:)) + call WriteVar(ncid, smpID, smp(:)) + call WriteVar(ncid, smp_alpha_ID, smp_alpha(:)) + call WriteVar(ncid, LFMCID, LFMC(:,:)) + call WriteVar(ncid, CBDID, CBD(:)) + call WriteVar(ncid, windID, wind(:)) + call WriteVar(ncid, dratioID, drying_ratio(:)) + call WriteVar(ncid, ROSACT_FM10_ID, ROS_active_FM10(:,:)) + call WriteVar(ncid, CI_FM10_ID, CI_FM10(:,:)) + + ! close file + call CloseNCFile(ncid) + + +end subroutine WriteCrownFireData \ No newline at end of file diff --git a/testing/functional_testing/fire/crownfire/crownfire_test.py b/testing/functional_testing/fire/crownfire/crownfire_test.py new file mode 100644 index 0000000000..2fa2fff426 --- /dev/null +++ b/testing/functional_testing/fire/crownfire/crownfire_test.py @@ -0,0 +1,239 @@ +""" +Concrete class for running the crownfire functional test for FATES. +""" +import os +import numpy as np +import xarray as xr +import pandas as pd +import matplotlib.pyplot as plt + +from functional_class import FunctionalTest +from utils_plotting import blank_plot + +COLORS = ["#793922", "#99291F", "#CC9728", "#6B8939", "#2C778A", "#2C378A"] +MPERMIN_TO_KMPERHOUR = 0.06 + + +class CrownFireTest(FunctionalTest): + """CrownFire test class""" + + name = "crownfire" + + def __init__(self, test_dict): + super().__init__( + CrownFireTest.name, + test_dict["test_dir"], + test_dict["test_exe"], + test_dict["out_file"], + test_dict["use_param_file"], + test_dict["other_args"], + ) + self.plot = True + + def plot_output(self, run_dir: str, save_figs: bool, plot_dir: str): + """Plot output associated with crown fire tests + + Args: + run_dir (str): run directory + out_file (str): output file + save_figs (bool): whether or not to save the figures + plot_dir (str): plot directory + """ + + # read in crown fire data + cfire_dat = xr.open_dataset(os.path.join(run_dir, self.out_file)) + + self.plot_passcrwn_fi(cfire_dat, save_figs, plot_dir) + self.plot_lfmc(cfire_dat, save_figs, plot_dir) + self.plot_rosact_fm10(cfire_dat, save_figs, plot_dir) + self.plot_ci_fm10(cfire_dat, save_figs, plot_dir) + + @staticmethod + def plot_passcrwn_fi(data: xr.Dataset, save_fig: bool, plot_dir: str = None): + """Plot min fire intensity required to ignite crown fuel + + Args: + data (xarray DataSet): the data set + save_fig (bool): whether or not to write out plot + plot_dir (str): if saving figure, where to write to + """ + data_frame = pd.DataFrame( + { + "CBH": np.tile(data.CBH, len(data.CWC)), + "CWC": np.repeat(data.CWC, len(data.CBH)), + "passive_crown_fi": data.passive_crown_fi.values.flatten(), + } + ) + + + max_CBH = data_frame["CBH"].max() + + plt.figure(figsize=(7, 5)) + axis = plt.subplot(111) + axis.spines["top"].set_visible(False) + axis.spines["bottom"].set_visible(False) + axis.spines["right"].set_visible(False) + axis.spines["left"].set_visible(False) + + axis.get_xaxis().tick_bottom() + axis.get_yaxis().tick_left() + + plt.xlim(0.0, max_CBH) + + plt.yticks(fontsize=10) + plt.xticks(fontsize=10) + plt.grid(True) + + CWC = np.unique(data_frame.CWC.values) + + for i, cwc in enumerate(CWC): + dat = data_frame[data_frame.CWC == cwc] + plt.plot( + dat.CBH.values, + dat["passive_crown_fi"].values, + lw=2, + color=COLORS[i], + label=cwc, + ) + + plt.xlabel("Canopy base height (m)", fontsize=11) + plt.ylabel("Threshold fire intensity to ignite crown fuels (kW m$^{-1}$)", fontsize=11) + plt.legend(loc="upper left", title="Canopy water content (%)") + + if save_fig: + fig_name = os.path.join(plot_dir, "passcrown_fi_plot.png") + plt.savefig(fig_name) + + @staticmethod + def plot_lfmc(data: xr.Dataset, save_fig: bool, plot_dir: str = None): + """Plot live fuel moisture content + + Args: + data (xarray DataSet): the data set + save_fig (bool): whether or not to write out plot + plot_dir (str): if saving figure, where to write to + """ + data_frame = pd.DataFrame( + { + "smp": np.tile(data.smp, len(data.smp_alpha)), + "smp_alpha": np.repeat(data.smp_alpha, len(data.smp)), + "LFMC": data.LFMC.values.flatten(), + } + ) + + min_smp = float(data_frame["smp"].min()) + max_lfmc = 150.0 + + blank_plot(0.0, -10.0, 150.0, 0.0, draw_horizontal_lines=True) + + smp_alpha_vals = np.unique(data_frame.smp_alpha.values) + + for i, alpha in enumerate(smp_alpha_vals): + dat = data_frame[data_frame["smp_alpha"] == alpha] + dat = dat.sort_values("smp") + + plt.plot( + dat.smp.values, + dat["LFMC"].values, + lw=2, + color=COLORS[i], + label=alpha, + ) + + plt.xlabel("Soil matric potential (MPa)", fontsize=11) + plt.ylabel("Live fuel moisture (%)", fontsize=11) + plt.legend(loc="upper right", title="Soil matric potential coeff") + plt.xlim(0.0,-10.0) + + if save_fig: + fig_name = os.path.join(plot_dir, "lfmc_plot.png") + plt.savefig(fig_name) + + @staticmethod + def plot_rosact_fm10(data: xr.Dataset, save_fig: bool, plot_dir: str = None): + """Plot active crown fire spread rate + + Args: + data (xarray DataSet): the data set + save_fig (bool): whether or not to write out plot + plot_dir (str): if saving figure, where to write to + """ + data_frame = pd.DataFrame( + { + "wind": np.tile(data.wind, len(data.dratio)), + "drying_ratio": np.repeat(data.dratio, len(data.wind)), + "ros_active_fm10": data.ROSACT_FM10.values.flatten(), + } + ) + + data_frame["wind_kmhr"] = data_frame.wind * MPERMIN_TO_KMPERHOUR # match usual wind speed unit in crown fire plots + + max_wind = data_frame["wind_kmhr"].max() + max_ros_active = data_frame["ros_active_fm10"].max() + + blank_plot(max_wind, 0.0, max_ros_active, 0.0, draw_horizontal_lines=True) + + dr_vals = np.unique(data_frame.drying_ratio.values) + for i, dr in enumerate(dr_vals): + dat = data_frame[data_frame["drying_ratio"] == dr] + + plt.plot( + dat.wind_kmhr.values, + dat["ros_active_fm10"].values, + lw=2, + color=COLORS[i], + label = dr, + ) + + plt.xlabel("Wind speed (km hr$^{-1}$)", fontsize=11) + plt.ylabel("Active crown fire ROS (m min$^{-1}$)", fontsize=11) + plt.legend(loc="upper right", title="Drying ratio") + + if save_fig: + fig_name = os.path.join(plot_dir, "ros_active_plot.png") + plt.savefig(fig_name) + + @staticmethod + def plot_ci_fm10(data: xr.Dataset, save_fig: bool, plot_dir: str = None): + """Plot crowning index + + Args: + data (xarray DataSet): the data set + save_fig (bool): whether or not to write out plot + plot_dir (str): if saving figure, where to write to + """ + data_frame = pd.DataFrame( + { + "CBD": np.tile(data.CBD, len(data.dratio)), + "drying_ratio": np.repeat(data.dratio, len(data.CBD)), + "CI_FM10": data.CI_FM10.values.flatten(), + } + ) + + data_frame["CI_kmhr"] = data_frame.CI_FM10* MPERMIN_TO_KMPERHOUR + + max_CBD = data_frame["CBD"].max() + max_CI = data_frame["CI_kmhr"].max() + + blank_plot(max_CBD, 0.0, max_CI, 0.0, draw_horizontal_lines=True) + + dr_vals = np.unique(data_frame.drying_ratio.values) + + for i, dr in enumerate(dr_vals): + dat = data_frame[data_frame["drying_ratio"] == dr] + + plt.plot( + dat.CBD.values, + dat["CI_kmhr"].values, + lw=2, + color=COLORS[i], + label = dr + ) + + plt.xlabel("Canopy bulk density (kg m$^{-3}$)", fontsize=11) + plt.ylabel("Crowning index (km hr$^{-1}$)", fontsize=11) + plt.legend(loc = "upper right", title="Drying ratio") + + if save_fig: + fig_name = os.path.join(plot_dir, "ci_plot.png") + plt.savefig(fig_name) \ No newline at end of file diff --git a/testing/functional_testing/fire/shr/CMakeLists.txt b/testing/functional_testing/fire/shr/CMakeLists.txt index b3ccabb546..53df61e470 100644 --- a/testing/functional_testing/fire/shr/CMakeLists.txt +++ b/testing/functional_testing/fire/shr/CMakeLists.txt @@ -1,6 +1,7 @@ list(APPEND fates_sources FatesTestFireMod.F90 SyntheticFuelModels.F90 + FatesTestCrownFireMod.F90 ) sourcelist_to_parent(fates_sources) \ No newline at end of file diff --git a/testing/functional_testing/fire/shr/FatesTestCrownFireMod.F90 b/testing/functional_testing/fire/shr/FatesTestCrownFireMod.F90 new file mode 100644 index 0000000000..4f66fd9b82 --- /dev/null +++ b/testing/functional_testing/fire/shr/FatesTestCrownFireMod.F90 @@ -0,0 +1,315 @@ + +module FatesTestCrownFireMod + ! + ! DESCRIPTION + ! modules to support testing crown fire + ! + + use FatesConstantsMod, only : r8 => fates_r8, nearzero + use FatesUnitTestIOMod, only : OpenNCFile, GetVar, CloseNCFile, RegisterNCDims + use FatesUnitTestIOMod, only : RegisterVar, EndNCDef, WriteVar + use FatesUnitTestIOMod, only : type_double, type_char, type_int + use SFParamsMod, only : SF_val_miner_total, SF_val_part_dens + use SFEquationsMod, only : OptimumPackingRatio, ReactionIntensity + use SFEquationsMod, only : HeatofPreignition, EffectiveHeatingNumber + use SFEquationsMod, only : WindFactor, PropagatingFlux + use SFEquationsMod, only : ForwardRateOfSpread + + implicit none + private + + public :: ROSWrapper + public :: EffectiveWindWrapper + public :: WriteCanopyFuelData + + ! ====================================================================================== + +contains + + subroutine ROSWrapper(bd, SAV, fuel_load, fmc, mef, fire_weather, effect_wind, ROS, i_r) + ! + ! DESCRIPTION + ! Calls a sequence of functions to calculate surface fire rate of spread + ! + ! ARGUMENTS: + real(r8), intent(in) :: bd ! fuel bulk density no trunk [kg/m3] + real(r8), intent(in) :: SAV ! fuel surface area to volume ratio [/cm] + real(r8), intent(in) :: fuel_load ! fuel amount excluding trunks [kgC/m2] + real(r8), intent(in) :: fmc ! weighted average fuel moisture content across non-trunk fuel classes [m3/m3] + real(r8), intent(in) :: mef ! weighted average of moisture of extinction across non-trunk fuel classes [m3/m3] + real(r8), intent(in) :: fire_weather ! fire weather index [unitless] + real(r8), intent(in) :: effect_wind ! effective wind speed [m/min] + real(r8), intent(out) :: ROS ! forward rate of spread [m/min] + real(r8), intent(out) :: i_r ! reaction intensity [kJ/m2/min], required for calculating HPA + ! + ! LOCALS: + real(r8) :: beta ! packing ratio [unitless] + real(r8) :: beta_op ! optimum packing ratio [unitless] + real(r8) :: beta_ratio ! relative packing ratio [unitless] + real(r8) :: xi ! propagating flux ratio [unitless] + real(r8) :: eps ! effective heating number [unitless] + real(r8) :: phi_wind ! wind factor [unitless] + real(r8) :: q_ig ! heat of pre-ignition [kJ/kg] + real(r8) :: fuel_net ! fuel load excluding mineral content [kgC/m2] + + ! fraction of fuel array occupied by fuels and optimum packing ratio + beta = bd / SF_val_part_dens + beta_op = OptimumPackingRatio(SAV) + ! relative packing ratio + if (beta_op < nearzero) then + beta_ratio = 0.0_r8 + else + beta_ratio = beta/beta_op + end if + + ! remove mineral content + fuel_net = fuel_load*(1.0_r8 - SF_val_miner_total) + + ! reaction intensity + i_r = ReactionIntensity(fuel_net/0.45_r8, SAV, beta_ratio, & + fmc, mef) + + ! heat of preignition + q_ig = HeatofPreignition(fmc) + + ! effective heating number + eps = EffectiveHeatingNumber(SAV) + + ! wind factor + phi_wind = WindFactor(effect_wind, beta_ratio, SAV) + + ! propogation flux + xi = PropagatingFlux(beta, SAV) + + ! forward rate of spread [m/min] + ROS = ForwardRateOfSpread(bd, eps, q_ig, i_r, xi, phi_wind) + + end subroutine ROSWrapper + + ! ====================================================================================== + + + subroutine EffectiveWindWrapper(tree_area, grass_area, area, wind, effect_wind) + ! + ! DESCRIPTION + ! Call a sequence of calculations to get effective wind speed + ! + ! ARGUMENTS + real(r8), intent(in) :: tree_area ! total tree area of the test patch [m2] + real(r8), intent(in) :: grass_area ! total grass area of the test patch [m2] + real(r8), intent(in) :: area ! patch area [m2] + real(r8), intent(in) :: wind ! wind speed [m/min] + real(r8), intent(out) :: effect_wind ! returned effective wind speed [m/min] + ! + ! LOCALS: + real(r8) :: tree_frac ! tree fraction of the test patch [0-1] + real(r8) :: grass_frac ! grass fraction of the test patch [0-1] + real(r8) :: bare_frac ! bare ground fraction of the test patch [0-1] + ! + ! CONSTANTS: + real(r8), parameter :: wind_atten_tree = 0.4_r8 ! wind attenuation factor for tree fraction + real(r8), parameter :: wind_atten_grass = 0.6_r8 ! wind attenuation factor for grass fraction + + tree_frac = tree_area / area + grass_frac = grass_area / area + grass_frac = min(grass_frac, 1.0_r8 - tree_frac) + bare_frac = 1.0_r8 - tree_frac - grass_frac + + effect_wind = wind * (tree_frac*wind_atten_tree + & + (grass_frac + bare_frac)*wind_atten_grass ) + + + end subroutine EffectiveWindWrapper + + ! ====================================================================================== + + subroutine WriteCanopyFuelData(out_file, nfuelmods, nstands, nwind, Wind, NI, CWC, & + CBD, CBH, canopy_fuel_load, ROS_front, FI, FI_init, ROS_actfm10, ROS_critical, ROS_actCI, & + CFB, ROS_final, FI_final, fuel_models, patch_types) + ! + ! DESCRIPTION: + ! write out data from canopy fuel functional test + ! + + ! ARGUMENTS: + character(len=*), intent(in) :: out_file + integer, intent(in) :: nfuelmods + integer, intent(in) :: nstands + integer, intent(in) :: nwind + real(r8), intent(in) :: Wind(:) + real(r8), intent(in) :: NI(:) + real(r8), intent(in) :: CWC(:) + real(r8), intent(in) :: CBD(:) + real(r8), intent(in) :: CBH(:) + real(r8), intent(in) :: canopy_fuel_load(:) + real(r8), intent(in) :: ROS_front(:,:,:,:) + real(r8), intent(in) :: FI(:,:,:,:) + real(r8), intent(in) :: FI_init(:,:) + real(r8), intent(in) :: ROS_actfm10(:,:,:,:,:) + real(r8), intent(in) :: ROS_critical(:) + real(r8), intent(in) :: ROS_actCI(:,:,:,:,:) + real(r8), intent(in) :: CFB(:,:,:,:,:) + real(r8), intent(in) :: ROS_final(:,:,:,:,:) + real(r8), intent(in) :: FI_final(:,:,:,:,:) + integer, intent(in) :: fuel_models(:) + integer, intent(in) :: patch_types(:) + + ! LOCALS: + integer :: ncid ! netcdf id + character(len=20) :: dim_names(5) ! dimension names + integer :: dimIDs(5) ! dimension IDs + integer :: modID, patchID, windID, niID, cwcID + integer :: CBDID, CBHID, cflID + integer :: ROS_frontID, ROS_actID,ROS_ciID, ROS_minID, ROS_finID + integer :: FIID, FI_initID, FI_finID + integer :: CFBID + + ! dimension names + dim_names = [character(len=20) :: 'wind','fire_weather','canopy_water', & + 'patch_type', 'fuel_model'] + ! open file + call OpenNCFile(trim(out_file), ncid, 'readwrite') + + ! register dimensions + call RegisterNCDims(ncid, dim_names, (/nwind, size(NI), size(CWC), & + nstands, nfuelmods/), 5, dimIDs) + + ! first register dimension variables + + ! register wind speed + call RegisterVar(ncid, 'wind', dimIDs(1:1), type_int, & + [character(len=20) :: 'units', 'long_name' ], & + [character(len=150) :: 'm min-1', 'wind speed index'], 2, windID) + + ! register fire weather index + call RegisterVar(ncid, 'fire_weather', dimIDs(2:2), type_int, & + [character(len=20) :: 'units', 'long_name' ], & + [character(len=150) :: '', 'fire weather index'], 2, niID) + + ! register canopy water content + call RegisterVar(ncid, 'canopy_water', dimIDs(3:3), type_int, & + [character(len=20) :: 'units', 'long_name' ], & + [character(len=150) :: '%', 'canopy water index'], 2, cwcID) + + ! register patch types + call RegisterVar(ncid, 'patch_type', dimIDs(4:4), type_int, & + [character(len=20) :: 'units', 'long_name' ], & + [character(len=150) :: '', 'patch type index'], 2, patchID) + + ! register fuel models + call RegisterVar(ncid, 'fuel_model', dimIDs(5:5), type_int, & + [character(len=20) :: 'units', 'long_name' ], & + [character(len=150) :: '', 'fuel model index'], 2, modID) + + ! register variables + + ! register actual variables + + ! register CBD + call RegisterVar(ncid, 'CBD', dimIDs(4:4), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'patch_type', 'kg m-3', 'canopy bulk density'], & + 3, CBDID) + + ! register CBH + call RegisterVar(ncid, 'CBH', dimIDs(4:4), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'patch_type', 'm', 'canopy base height'], & + 3, CBHID) + + ! register canopy fuel load + call RegisterVar(ncid, 'Canopy_fuel', dimIDs(4:4), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'patch_type', 'kg m-2', 'canopy fuel density'], & + 3, cflID) + + ! register surface fire spread rate + call RegisterVar(ncid, 'ROS_front', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type fuel_model', 'm min-1', & + 'surface forward rate of spread'], 3, ROS_frontID) + + ! register surface fire intensity + call RegisterVar(ncid, 'FI', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type fuel_model', 'kW m-1', & + 'surface fire intensity'], 3, FIID) + + ! register min FI to initiate crown fire + call RegisterVar(ncid, 'FI_init', (/dimIDs(4), dimIDs(3)/), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'patch_type canopy_water', 'kW m-1', & + 'fire intensity to ignite crown fuel'], 3, FI_initID) + + ! register active crown fire ROS + call RegisterVar(ncid, 'ROS_active', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(3), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type canopy_water fuel_model', & + 'm min-1', 'active crown fire ROS assuming FM10'], 3, ROS_actID) + + ! register active crown fire ROS at CI wind speed + call RegisterVar(ncid, 'ROS_act_ci', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(3), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type canopy_water fuel_model', & + 'm min-1', 'active crown fire ROS assuming FM10 at crowning index ws'], 3, ROS_ciID) + + ! register critical ROS to sustain active crown fire + call RegisterVar(ncid, 'ROS_min', dimIDs(4:4), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'patch_type', 'm min-1', 'min ROS to sustain active crown fire'], & + 3, ROS_minID) + + ! register crown fraction burnt + call RegisterVar(ncid, 'CFB', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(3), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type canopy_water fuel_model', '', & + 'crown fraction burnt by crown fire'], 3, CFBID) + + ! register final ROS + call RegisterVar(ncid, 'ROS_final', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(3), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type canopy_water fuel_model', 'm min-1', & + 'ROS after crown fire update'], 3, ROS_finID) + + ! register final FI + call RegisterVar(ncid, 'FI_final', (/dimIDs(1), dimIDs(2), dimIDs(4), dimIDs(3), dimIDs(5)/), & + type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'wind fire_weather patch_type canopy_water fuel_model', & + 'kW m-1', 'FI after crown fire update'], 3, FI_finID) + + ! finish defining variables + call EndNCDef(ncid) + + ! write out data + call WriteVar(ncid, patchID, patch_types(:)) + call WriteVar(ncid, modID, fuel_models(:)) + call WriteVar(ncid, windID, Wind(:)) + call WriteVar(ncid, niID, NI(:)) + call WriteVar(ncid, cwcID, CWC(:)) + call WriteVar(ncid, CBDID, CBD(:)) + call WriteVar(ncid, CBHID, CBH(:)) + call WriteVar(ncid, cflID, canopy_fuel_load(:)) + call WriteVar(ncid, ROS_frontID, ROS_front(:,:,:,:)) + call WriteVar(ncid, FIID, FI(:,:,:,:)) + call WriteVar(ncid, FI_initID, FI_init(:,:)) + call WriteVar(ncid, ROS_actID, ROS_actfm10(:,:,:,:,:)) + call WriteVar(ncid, ROS_minID, ROS_critical(:)) + call WriteVar(ncid, ROS_ciID, ROS_actCI(:,:,:,:,:)) + call WriteVar(ncid, CFBID, CFB(:,:,:,:,:)) + call WriteVar(ncid, ROS_finID, ROS_final(:,:,:,:,:)) + call WriteVar(ncid, FI_finID, FI_final(:,:,:,:,:)) + + call CloseNCFile(ncid) + + end subroutine WriteCanopyFuelData + + + +end module FatesTestCrownFireMod diff --git a/testing/functional_tests.cfg b/testing/functional_tests.cfg index fec8e075d6..c8f764487b 100644 --- a/testing/functional_tests.cfg +++ b/testing/functional_tests.cfg @@ -40,3 +40,17 @@ test_exe = FATES_firemort_exe out_file = fire_mortality_out.nc use_param_file = True other_args = [] + +[crownfire] +test_dir = fates_crownfire_ftest +test_exe = FATES_crownfire_exe +out_file = crownfire_out.nc +use_param_file = True +other_args = [] + +[canopyfuel] +test_dir = fates_canopyfuel_ftest +test_exe = FATES_canopyfuel_exe +out_file = canopyfuel_out.nc +use_param_file = True +other_args = [] diff --git a/testing/load_functional_tests.py b/testing/load_functional_tests.py index 3b72dcb2b7..6e575aa3b3 100644 --- a/testing/load_functional_tests.py +++ b/testing/load_functional_tests.py @@ -8,3 +8,5 @@ from functional_testing.fire.ros.ros_test import ROSTest from functional_testing.patch.patch_test import PatchTest from functional_testing.fire.mortality.fire_mortality_test import FireMortTest +from functional_testing.fire.crownfire.crownfire_test import CrownFireTest +from functional_testing.fire.canopyfuel.canopyfuel_test import CanopyFuelTest diff --git a/testing/testing_shr/FatesUnitTestIOMod.F90 b/testing/testing_shr/FatesUnitTestIOMod.F90 index ef261f1d0f..8150265532 100644 --- a/testing/testing_shr/FatesUnitTestIOMod.F90 +++ b/testing/testing_shr/FatesUnitTestIOMod.F90 @@ -26,6 +26,8 @@ module FatesUnitTestIOMod module procedure WriteVar1DReal module procedure WriteVar2DReal module procedure WriteVar3DReal + module procedure WriteVar4DReal + module procedure WriteVar5DReal module procedure WriteVar1DInt module procedure WriteVar2DInt module procedure WriteVar1DChar @@ -561,6 +563,40 @@ end subroutine WriteVar3DReal ! ===================================================================================== + subroutine WriteVar4DReal(ncid, varID, data) + ! + ! DESCRIPTION: + ! Write 3D real data + ! + + ! ARGUMENTS: + integer, intent(in) :: ncid ! netcdf file id + integer, intent(in) :: varID ! variable ID + real(r8), intent(in) :: data(:,:,:,:) ! data to write + + call Check(nf90_put_var(ncid, varID, data(:,:,:,:))) + + end subroutine WriteVar4DReal + + ! ===================================================================================== + + subroutine WriteVar5DReal(ncid, varID, data) + ! + ! DESCRIPTION: + ! Write 5d real data + + ! ARGUMENTS: + integer, intent(in) :: ncid ! netcdf file id + integer, intent(in) :: varID ! variable ID + real(r8), intent(in) :: data(:,:,:,:,:) ! data to write + + call Check(nf90_put_var(ncid, varID, data(:,:,:,:,:))) + + end subroutine WriteVar5DReal + + ! ===================================================================================== + + subroutine WriteVar1DInt(ncid, varID, data) ! ! DESCRIPTION: @@ -629,4 +665,4 @@ end subroutine WriteVar2DChar ! ===================================================================================== -end module FatesUnitTestIOMod \ No newline at end of file +end module FatesUnitTestIOMod diff --git a/testing/testing_shr/SyntheticPatchTypes.F90 b/testing/testing_shr/SyntheticPatchTypes.F90 index 64747e5b6d..e7f3877255 100644 --- a/testing/testing_shr/SyntheticPatchTypes.F90 +++ b/testing/testing_shr/SyntheticPatchTypes.F90 @@ -234,6 +234,57 @@ subroutine GetSyntheticPatchData(this) pft_ids=(/6, 2, 2, 9/), & canopy_layers=(/1, 1, 2, 2/)) + call this%AddPatch(patch_id=6, patch_name='low density', area=500.0_r8, & + ages=(/100.0_r8, 80.0_r8, 80.0_r8, 50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/90.0_r8, 70.0_r8, 70.0_r8, 40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 10.0_r8/), & + densities=(/0.0001_r8, 0.002_r8, 0.002_r8, 0.005_r8, 0.01_r8, 0.0001_r8, 0.0001_r8, 0.0001_r8/), & + pft_ids=(/1, 1, 2, 3, 2, 3, 4, 4/), & + canopy_layers=(/1, 1, 1, 1, 1, 2, 2, 2/)) + + call this%AddPatch(patch_id=7, patch_name='high density', area=500.0_r8, & + ages=(/100.0_r8, 80.0_r8, 80.0_r8, 50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/90.0_r8, 70.0_r8, 70.0_r8, 40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 10.0_r8/), & + densities=(/0.0001_r8, 0.002_r8, 0.002_r8, 0.005_r8, 0.01_r8, 0.01_r8, 0.01_r8, 0.01_r8/), & + pft_ids=(/1, 1, 2, 3, 2, 3, 4, 4/), & + canopy_layers=(/1, 1, 1, 1, 1, 2, 2, 2/)) + + call this%AddPatch(patch_id=8, patch_name='young stand', area=500.0_r8, & + ages=(/50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 30.0_r8/), & + dbhs=(/40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 15.0_r8/), & + densities=(/0.006_r8, 0.01_r8, 0.006_r8, 0.01_r8, 0.01_r8/), & + pft_ids=(/1, 2, 3, 3, 4/), & + canopy_layers=(/1, 1, 1, 2, 2/)) + + call this%AddPatch(patch_id=9, patch_name='low pine', area=500.0_r8, & + ages=(/100.0_r8, 80.0_r8, 80.0_r8, 50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/90.0_r8, 70.0_r8, 70.0_r8, 40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 10.0_r8/), & + densities=(/0.0001_r8, 0.002_r8, 0.002_r8, 0.005_r8, 0.01_r8, 0.0001_r8, 0.0001_r8, 0.0001_r8/), & + pft_ids=(/1, 1, 1, 1, 1, 1, 4, 4/), & + canopy_layers=(/1, 1, 1, 1, 1, 2, 2, 2/)) + + call this%AddPatch(patch_id=10, patch_name='high pine', area=500.0_r8, & + ages=(/100.0_r8, 80.0_r8, 80.0_r8, 50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/90.0_r8, 70.0_r8, 70.0_r8, 40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 10.0_r8/), & + densities=(/0.0001_r8, 0.002_r8, 0.002_r8, 0.005_r8, 0.01_r8, 0.01_r8, 0.01_r8, 0.01_r8/), & + pft_ids=(/1, 1, 1, 1, 1, 1, 4, 4/), & + canopy_layers=(/1, 1, 1, 1, 1, 2, 2, 2/)) + + call this%AddPatch(patch_id=11, patch_name='low fir', area=500.0_r8, & + ages=(/100.0_r8, 80.0_r8, 80.0_r8, 50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/90.0_r8, 70.0_r8, 70.0_r8, 40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 10.0_r8/), & + densities=(/0.0001_r8, 0.002_r8, 0.002_r8, 0.005_r8, 0.01_r8, 0.0001_r8, 0.0001_r8, 0.0001_r8/), & + pft_ids=(/3, 3, 3, 3, 3, 3, 4, 4/), & + canopy_layers=(/1, 1, 1, 1, 1, 2, 2, 2/)) + + call this%AddPatch(patch_id=12, patch_name='high fir', area=500.0_r8, & + ages=(/100.0_r8, 80.0_r8, 80.0_r8, 50.0_r8, 50.0_r8, 20.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/90.0_r8, 70.0_r8, 70.0_r8, 40.0_r8, 30.0_r8, 20.0_r8, 15.0_r8, 10.0_r8/), & + densities=(/0.0001_r8, 0.002_r8, 0.002_r8, 0.005_r8, 0.01_r8, 0.01_r8, 0.01_r8, 0.01_r8/), & + pft_ids=(/3, 3, 3, 3, 3, 3, 4, 4/), & + canopy_layers=(/1, 1, 1, 1, 1, 2, 2, 2/)) + + + end subroutine GetSyntheticPatchData ! --------------------------------------------------------------------------------------