Crown fire Scott & Reinhardt 2001

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
 

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

fire/CMakeLists.txt

@@ -6,6 +6,7 @@ list(APPEND fates_sources
   FatesFuelMod.F90
   FatesFuelClassesMod.F90
   SFEquationsMod.F90
+  CrownFireEquationsMod.F90
   )
 
 sourcelist_to_parent(fates_sources)

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