Merge pull request #1326 from adrifoster/fire_refactor_lightning

Refactor SPITFIRE - overhaul of fire behavior sections (but bit4bit)

Author: glemieux

biogeochem/EDPatchDynamicsMod.F90

@@ -434,7 +434,8 @@ subroutine disturbance_rates( site_in, bc_in)
        endif
 
        ! Fire Disturbance Rate
-       currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt
+       currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt 
+
 
        ! Fires can't burn the whole patch, as this causes /0 errors. 
        if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then

fire/FatesFuelMod.F90

@@ -34,6 +34,8 @@ module FatesFuelMod
       procedure :: UpdateFuelMoisture
       procedure :: AverageBulkDensity_NoTrunks
       procedure :: AverageSAV_NoTrunks
+      procedure :: CalculateFuelBurnt
+      procedure :: CalculateResidenceTime
 
   end type fuel_type
 
@@ -354,8 +356,8 @@ subroutine AverageSAV_NoTrunks(this, sav_fuel)
       !    Pyne et al., 1996 (Introduction to wildland fire)
 
       ! ARGUMENTS:
-      class(fuel_type),   intent(inout) :: this                       ! fuel class
-      real(r8),           intent(in)    :: sav_fuel(num_fuel_classes) ! surface area to volume ratio of all fuel types [/cm]
+      class(fuel_type), intent(inout) :: this                       ! fuel class
+      real(r8),         intent(in)    :: sav_fuel(num_fuel_classes) ! surface area to volume ratio of all fuel types [/cm]
 
       ! LOCALS:
       integer :: i ! looping index
@@ -376,4 +378,96 @@ end subroutine AverageSAV_NoTrunks
 
   !---------------------------------------------------------------------------------------
 
+    subroutine CalculateFuelBurnt(this, fuel_consumed)
+      ! DESCRIPTION:
+      !   Calculates the fraction and total amount of fuel burnt
+      !
+      
+      use SFParamsMod, only : SF_val_mid_moisture, SF_val_mid_moisture_Coeff
+      use SFParamsMod, only : SF_val_mid_moisture_Slope, SF_val_min_moisture
+      use SFParamsMod, only : SF_val_low_moisture_Coeff, SF_val_low_moisture_Slope
+      use SFParamsMod, only : SF_val_miner_total
+
+      ! ARGUMENTS:
+      class(fuel_type), intent(inout) :: this                            ! fuel class
+      real(r8),         intent(out)   :: fuel_consumed(num_fuel_classes) ! fuel consumed [kgC/m2]
+      
+      ! LOCALS:
+      real(r8) :: rel_moisture  ! relative moisture of fuel (moist/moisture of extinction) [unitless]
+      integer  :: i             ! looping index
+      
+      ! CONSTANTS:
+      real(r8), parameter :: max_grass_frac = 0.8_r8 ! maximum fraction burnt for live grass fuels
+      
+      this%frac_burnt(:) = 1.0_r8
+        
+      ! Calculate fraction of litter is burnt for all classes. 
+      ! Equation B1 in Thonicke et al. 2010
+      do i = 1, num_fuel_classes        
+        
+        rel_moisture = this%effective_moisture(i)                  
+        
+        if (rel_moisture <= SF_val_min_moisture(i)) then
+          ! very dry litter
+          this%frac_burnt(i) = 1.0_r8 
+        else if (rel_moisture > SF_val_min_moisture(i) .and. rel_moisture <= SF_val_mid_moisture(i)) then
+          ! low to medium moisture
+          this%frac_burnt(i) = max(0.0_r8, min(1.0_r8, SF_val_low_moisture_Coeff(i) - &
+            SF_val_low_moisture_Slope(i)*rel_moisture))
+        else if (rel_moisture > SF_val_mid_moisture(i) .and. rel_moisture <= 1.0_r8) then
+          ! medium to high moisture
+          this%frac_burnt(i) = max(0.0_r8, min(1.0_r8, SF_val_mid_moisture_Coeff(i) - &
+            SF_val_mid_moisture_Slope(i)*rel_moisture))
+        else 
+          ! very wet litter
+          this%frac_burnt(i) = 0.0_r8  
+        endif
+        
+        ! we can't ever kill all of the grass
+        if (i == fuel_classes%live_grass()) then
+          this%frac_burnt(i) = min(max_grass_frac, this%frac_burnt(i))
+        end if
+        
+        ! reduce fraction burnt based on mineral content
+        this%frac_burnt(i) = this%frac_burnt(i)*(1.0_r8 - SF_val_miner_total)
+        
+        ! calculate fuel consumed
+        fuel_consumed(i) = this%frac_burnt(i)*this%loading(i)
+      end do
+
+    end subroutine CalculateFuelBurnt
+    
+    !-------------------------------------------------------------------------------------
+    
+    subroutine CalculateResidenceTime(this, tau_l)
+      !
+      !  DESCRIPTION:
+      !  Calculates fire residence time, duration of lethal bole heating [min]
+      !  This is used for determining cambial kill of woody cohorts
+      !
+      !  From Peterson & Ryan (1986)
+      !
+      
+      ! ARGUMENTS:
+      class(fuel_type), intent(in)  :: this  ! fuel class
+      real(r8),         intent(out) :: tau_l ! duration of lethal bole heating [min]
+      
+      ! LOCALS:
+      integer :: i ! looping index
+      
+      tau_l = 0.0_r8 
+      do i = 1, num_fuel_classes
+        if (i /= fuel_classes%trunks()) then 
+          ! don't include 1000-hr fuels
+          ! convert loading from kgC/m2 to g/cm2
+          tau_l = tau_l + 39.4_r8*(this%frac_loading(i)*this%non_trunk_loading/0.45_r8/10.0_r8)* &
+            (1.0_r8 - ((1.0_r8 - this%frac_burnt(i))**0.5_r8))
+        end if 
+      end do
+      
+      ! cap the residence time to 8mins, as suggested by literature survey by P&R (1986)
+      tau_l = min(8.0_r8, tau_l) 
+
+    end subroutine CalculateResidenceTime
+    
 end module FatesFuelMod

fire/SFEquationsMod.F90

@@ -16,13 +16,19 @@ module SFEquationsMod
   public :: MaximumReactionVelocity
   public :: OptimumReactionVelocity
   public :: OptimumPackingRatio
+  public :: MoistureCoefficient
   public :: ReactionIntensity
   public :: HeatofPreignition
   public :: EffectiveHeatingNumber
   public :: WindFactor
   public :: PropagatingFlux
   public :: ForwardRateOfSpread
   public :: BackwardRateOfSpread
+  public :: FireDuration
+  public :: LengthToBreadth
+  public :: FireSize
+  public :: AreaBurnt
+  public :: FireIntensity
 
   contains 
 
@@ -62,7 +68,11 @@ real(r8) function MaximumReactionVelocity(SAV)
       ! ARGUMENTS:
       real(r8), intent(in) :: SAV ! fuel surface area to volume ratio [/cm]
 
-      MaximumReactionVelocity = 1.0_r8/(0.0591_r8 + 2.926_r8*(SAV**(-1.5_r8)))
+      if (SAV < nearzero) then 
+        MaximumReactionVelocity = 0.0_r8
+      else 
+        MaximumReactionVelocity = 1.0_r8/(0.0591_r8 + 2.926_r8*(SAV**(-1.5_r8)))
+      end if
 
     end function MaximumReactionVelocity
 
@@ -112,6 +122,13 @@ real(r8) function MoistureCoefficient(moisture, MEF)
       ! LOCALS:
       real(r8) :: mw_weight ! relative fuel moisture/fuel moisture of extinction
 
+      if (MEF < nearzero) then 
+        ! this really should never happen - essentially this means fuel can never burn
+        ! but we are putting this here to avoid divide by zeros 
+        MoistureCoefficient = 0.0_r8
+        return
+      end if 
+      
       ! average values for litter pools (dead leaves, twigs, small and large branches), plus grass
       mw_weight = moisture/MEF
 
@@ -314,5 +331,137 @@ real(r8) function BackwardRateOfSpread(ros_front, wind_speed)
     end function BackwardRateOfSpread
 
     !-------------------------------------------------------------------------------------
+    
+    real(r8) function FireDuration(FDI)
+      !
+      !  DESCRIPTION:
+      !  Calculates fire duration [min]
+      !
+      !  Equation 14 in Thonicke et al. 2010
+      !
+      use SFParamsMod, only : SF_val_max_durat, SF_val_durat_slope
+
+      ! ARGUMENTS:
+      real(r8), intent(in) :: FDI  ! fire danger index [0-1]
+
+      FireDuration = (SF_val_max_durat + 1.0_r8)/(1.0_r8 + SF_val_max_durat*             &
+        exp(SF_val_durat_slope*FDI))
+
+    end function FireDuration
+
+    !-------------------------------------------------------------------------------------
+  
+    real(r8) function LengthToBreadth(effective_windspeed, tree_fraction)
+      !
+      !  DESCRIPTION:
+      !  Calculates length to breadth ratio of fire ellipse [unitless], used for calculating area burnt 
+      !
+      !     Canadian Forest Fire Behavior Prediction System Ont.Int.Rep. ST-X-3, 1992
+      !     Information Report GLC-X-10, Wotten et al. 2009
+      !
+      use FatesConstantsMod, only : m_per_km, min_per_hr
+
+      ! ARGUMENTS:
+      real(r8), intent(in) :: effective_windspeed ! effective windspeed [m/min]
+      real(r8), intent(in) :: tree_fraction       ! tree fraction [0-1]
+      
+      ! LOCALS:
+      real(r8) :: windspeed_km_hr ! effective windspeed, converted to correct units [km/hr]
+            
+      ! CONSTANTS:
+      real(r8), parameter :: lb_threshold = 0.55_r8 ! tree canopy fraction below which to use grassland length-to-breadth eqn
+
+      windspeed_km_hr = effective_windspeed/m_per_km*min_per_hr
+
+      if (windspeed_km_hr < 1.0_r8) then 
+        LengthToBreadth = 1.0_r8
+      else
+        if (tree_fraction > lb_threshold) then 
+          LengthToBreadth = 1.0_r8 + 8.729_r8*((1.0_r8 - exp(-0.03_r8*windspeed_km_hr))**2.155_r8)
+        else  
+          LengthToBreadth = 1.1_r8*(windspeed_km_hr**0.464_r8)
+        endif
+      endif
+
+    end function LengthToBreadth
+
+    !-------------------------------------------------------------------------------------
+    
+    real(r8) function FireSize(length_to_breadth, ros_back, ros_forward, fire_duration)
+      !
+      !  DESCRIPTION:
+      !  Calculates fire size [m2]
+      !
+      !  Eq 14 Arora and Boer JGR 2005 (area of an ellipse)
+      !
+      use FatesConstantsMod, only : pi_const
+
+      ! ARGUMENTS:
+      real(r8), intent(in) :: length_to_breadth ! length to breadth ratio of fire ellipse [unitless]
+      real(r8), intent(in) :: ros_back          ! backwards rate of spread [m/min]
+      real(r8), intent(in) :: ros_forward       ! forward rate of spread [m/min]
+      real(r8), intent(in) :: fire_duration     ! fire duration [min]
+      
+      ! LOCALS:
+      real(r8) :: dist_back    ! distance fire has travelled backwards [m]
+      real(r8) :: dist_forward ! distance fire has travelled forward [m]
+      real(r8) :: fire_size    ! area of fire [m2]
+      
+      dist_back = ros_back*fire_duration
+      dist_forward = ros_forward*fire_duration
+      
+      ! Eq 14 Arora and Boer JGR 2005 (area of an ellipse)
+      if (length_to_breadth < nearzero) then 
+        FireSize = 0.0_r8
+      else 
+        FireSize = (pi_const/(4.0_r8*length_to_breadth))*((dist_forward + dist_back)**2.0_r8)
+      end if
+      
+    end function FireSize
+    
+    !-------------------------------------------------------------------------------------
+    
+    real(r8) function AreaBurnt(fire_size, num_ignitions, FDI)
+      !
+      !  DESCRIPTION:
+      !  Calculates area burnt [m2/km2/day]
+      !
+      ! daily area burnt = size fires in m2 * num ignitions per day per km2 * prob ignition starts fire
+      ! Thonicke 2010 Eq. 1
+      !
+      ! the denominator in the units of currentSite%NF is total gridcell area, but since we assume that ignitions 
+      ! are equally probable across patches, currentSite%NF is equivalently per area of a given patch
+      ! thus AreaBurnt has units of m2 burned area per km2 patch area per day
+      !
+      ! TO DO: Connect here with the Li & Levis GDP fire suppression algorithm. 
+      !     Equation 16 in arora and boer model JGR 2005
+      !
+
+      ! ARGUMENTS:
+      real(r8), intent(in) :: fire_size     ! fire size [m2]
+      real(r8), intent(in) :: num_ignitions ! number of ignitions [/km2/day]
+      real(r8), intent(in) :: FDI           ! fire danger index [0-1]
+      
+      AreaBurnt = fire_size*num_ignitions*FDI
+
+    end function AreaBurnt
+    
+    !-------------------------------------------------------------------------------------
+    
+    real(r8) function FireIntensity(fuel_consumed, ros)
+      !
+      !  DESCRIPTION:
+      !  Calculates fire intensity [kW/m]
+      !  Eq 15 Thonicke et al 2010
+
+      use SFParamsMod, only : SF_val_fuel_energy
+
+      ! ARGUMENTS:
+      real(r8), intent(in) :: fuel_consumed ! fuel consumed [kg/m2]
+      real(r8), intent(in) :: ros           ! rate of spread [m/s]
+      
+      FireIntensity = SF_val_fuel_energy*fuel_consumed*ros
+
+    end function FireIntensity
 
 end module SFEquationsMod