[WIP] Edge forest flammability enhancements

biogeochem/EDCanopyStructureMod.F90

@@ -35,6 +35,7 @@ module EDCanopyStructureMod
   use FatesInterfaceTypesMod     , only : hlm_use_planthydro
   use FatesInterfaceTypesMod     , only : hlm_use_cohort_age_tracking
   use FatesInterfaceTypesMod     , only : hlm_use_sp
+  use FatesInterfaceTypesMod     , only : hlm_use_edge_forest
   use FatesInterfaceTypesMod     , only : numpft
   use FatesInterfaceTypesMod, only : bc_in_type
   use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
@@ -1319,6 +1320,9 @@ subroutine canopy_summarization( nsites, sites, bc_in )
     use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index
     use EDtypesMod                , only : area
     use FatesConstantsMod         , only : itrue
+    use FatesEcotypesMod     , only : IsPatchForest
+    use EDParamsMod          , only : forest_tree_fraction_threshold
+    use FatesEdgeForestMod   , only : CalculateEdgeForestArea
 
     ! !ARGUMENTS
     integer                 , intent(in)            :: nsites
@@ -1447,9 +1451,15 @@ subroutine canopy_summarization( nsites, sites, bc_in )
              currentPatch%total_canopy_area = currentPatch%area
           endif
 
+          currentPatch%is_forest = IsPatchForest(currentPatch, forest_tree_fraction_threshold)
+
           currentPatch => currentPatch%younger
        end do !patch loop
 
+       if (hlm_use_edge_forest == itrue) then
+          call CalculateEdgeForestArea(sites(s))
+       end if
+
        call leaf_area_profile(sites(s))
 
        if(hlm_radiation_model.eq.twostr_solver) then

biogeochem/EDPatchDynamicsMod.F90

@@ -680,7 +680,8 @@ subroutine spawn_patches( currentSite, bc_in )
 
                    call newPatch%Create(age, site_areadis, i_landusechange_receiverpatchlabel, i_nocomp_pft, &
                                          num_swb, numpft, currentSite%nlevsoil, hlm_current_tod,              &
-                                         hlm_regeneration_model)
+                                         hlm_regeneration_model, &
+                                         currentSite%youngest_patch)
 
                    ! Initialize the litter pools to zero, these
                    ! pools will be populated by looping over the existing patches
@@ -1442,7 +1443,7 @@ subroutine spawn_patches( currentSite, bc_in )
 
                 call buffer_patch%Create(0._r8, 0._r8, i_land_use_label, 0, &
                      num_swb, numpft, currentSite%nlevsoil, hlm_current_tod,              &
-                     hlm_regeneration_model)
+                     hlm_regeneration_model, currentSite%youngest_patch)
 
                 ! Initialize the litter pools to zero
                 do el=1,num_elements
@@ -1735,7 +1736,7 @@ subroutine split_patch(currentSite, currentPatch, new_patch, fraction_to_keep, a
     call new_patch%Create(0._r8, temp_area, &
          currentPatch%land_use_label, currentPatch%nocomp_pft_label, &
          num_swb, numpft, currentSite%nlevsoil, hlm_current_tod, &
-         hlm_regeneration_model)
+         hlm_regeneration_model, currentSite%youngest_patch)
 
     ! Initialize the litter pools to zero, these
     ! pools will be populated shortly
@@ -3272,6 +3273,14 @@ subroutine fuse_2_patches(csite, dp, rp)
     ! Radiation
     rp%rad_error(1)         = (dp%rad_error(1)*dp%area + rp%rad_error(1)*rp%area) * inv_sum_area
     rp%rad_error(2)         = (dp%rad_error(2)*dp%area + rp%rad_error(2)*rp%area) * inv_sum_area
+
+    ! Fire weather
+    rp%fireWeather%precip = (dp%fireWeather%precip*dp%area + rp%fireWeather%precip*rp%area) * inv_sum_area
+    rp%fireWeather%rh = (dp%fireWeather%rh*dp%area + rp%fireWeather%rh*rp%area) * inv_sum_area
+    rp%fireWeather%temp_C = (dp%fireWeather%temp_C*dp%area + rp%fireWeather%temp_C*rp%area) * inv_sum_area
+    rp%fireWeather%wind = (dp%fireWeather%wind*dp%area + rp%fireWeather%wind*rp%area) * inv_sum_area
+    rp%fireWeather%fire_weather_index = (dp%fireWeather%fire_weather_index*dp%area + rp%fireWeather%fire_weather_index*rp%area) * inv_sum_area
+    rp%fireWeather%effective_windspeed = (dp%fireWeather%effective_windspeed*dp%area + rp%fireWeather%effective_windspeed*rp%area) * inv_sum_area
 
     rp%area = rp%area + dp%area !THIS MUST COME AT THE END!
 

biogeochem/FatesPatchMod.F90

@@ -32,6 +32,9 @@ module FatesPatchMod
   use FatesRadiationMemMod,   only : num_rad_stream_types
   use FatesInterfaceTypesMod, only : hlm_hio_ignore_val
   use FatesInterfaceTypesMod, only : numpft
+  use FatesInterfaceTypesMod, only : nlevedgeforest
+  use SFFireWeatherMod,       only : fire_weather
+  use SFNesterovMod,          only : nesterov_index
   use shr_infnan_mod,         only : nan => shr_infnan_nan, assignment(=)
   use shr_log_mod,            only : errMsg => shr_log_errMsg
 
@@ -74,6 +77,12 @@ module FatesPatchMod
 
     !---------------------------------------------------------------------------
 
+    ! FOREST INFO
+    logical  :: is_forest                    ! whether the patch is "forest" according to FATES param file criteria
+    real(r8), dimension(:), allocatable :: area_in_edgeforest_bins
+
+    !---------------------------------------------------------------------------
+
     ! RUNNING MEANS
     !class(rmean_type),    pointer :: t2m                  ! place-holder for 2m air temperature (variable window-size)
     class(rmean_type),     pointer :: tveg24               ! 24-hour mean vegetation temperature [K]
@@ -201,9 +210,17 @@ module FatesPatchMod
     !---------------------------------------------------------------------------
 
     ! FUELS AND FIRE
+    ! fire weather
+    class(fire_weather), pointer :: fireWeather      ! fire weather object
+
     ! fuel characteristics
     real(r8)              :: livegrass               ! total aboveground grass biomass in patch [kgC/m2]
 
+    ! number of fires
+    real(r8) ::  fdi                                 ! daily probability an ignition event will start a fire
+    real(r8) ::  NF                                  ! daily ignitions in km2
+    real(r8) ::  NF_successful                       ! daily ignitions in km2 that actually lead to fire
+
     ! fire spread
     real(r8)              :: ros_front               ! rate of forward  spread of fire [m/min]
     real(r8)              :: ros_back                ! rate of backward spread of fire [m/min]
@@ -258,7 +275,7 @@ module FatesPatchMod
 
     !===========================================================================
 
-    subroutine Init(this, num_swb, num_levsoil)
+    subroutine Init(this, num_swb, num_levsoil, from_patch_for_fire)
       !
       !  DESCRIPTION:
       !  Initialize a new patch - allocate arrays and set values to nan and/or 0.0
@@ -268,6 +285,7 @@ subroutine Init(this, num_swb, num_levsoil)
       class(fates_patch_type), intent(inout) :: this        ! patch object
       integer,                 intent(in)    :: num_swb     ! number of shortwave broad-bands to track
       integer,                 intent(in)    :: num_levsoil ! number of soil layers
+      class(fates_patch_type), intent(in), optional :: from_patch_for_fire  ! patch to take fire info from
 
       ! allocate arrays 
       allocate(this%tr_soil_dir(num_swb))
@@ -279,6 +297,21 @@ subroutine Init(this, num_swb, num_levsoil)
       allocate(this%sabs_dir(num_swb))
       allocate(this%sabs_dif(num_swb))
       allocate(this%fragmentation_scaler(num_levsoil))
+      allocate(this%area_in_edgeforest_bins(nlevedgeforest))
+
+      ! allocate and copy or initialize patch-level fire info
+      allocate(nesterov_index :: this%fireWeather)
+      if (present(from_patch_for_fire)) then
+        call this%fireWeather%CopyFrom(from_patch_for_fire%fireWeather)
+        this%FDI              = from_patch_for_fire%FDI
+        this%NF               = from_patch_for_fire%NF
+        this%NF_successful    = from_patch_for_fire%NF_successful
+      else
+        call this%fireWeather%Init()
+        this%FDI              = 0.0_r8
+        this%NF               = 0.0_r8
+        this%NF_successful    = 0.0_r8
+      end if
 
       ! initialize all values to nan
       call this%NanValues()
@@ -453,6 +486,10 @@ subroutine NanValues(this)
       this%ncl_p                        = fates_unset_int
       this%land_use_label               = fates_unset_int
       this%age_since_anthro_disturbance = nan
+
+      ! FOREST INFO
+      this%is_forest                       = .false.
+      this%area_in_edgeforest_bins(:) = nan
 
       ! LEAF ORGANIZATION
       this%pft_agb_profile(:,:)         = nan
@@ -706,7 +743,7 @@ end subroutine InitLitter
     !===========================================================================
 
     subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft,    &
-      num_levsoil, current_tod, regeneration_model)
+      num_levsoil, current_tod, regeneration_model, from_patch_for_fire)
       !
       ! DESCRIPTION:
       ! create a new patch with input and default values
@@ -723,10 +760,15 @@ subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft,
       integer,                 intent(in)    :: num_levsoil        ! number of soil layers
       integer,                 intent(in)    :: current_tod        ! time of day [seconds past 0Z]
       integer,                 intent(in)    :: regeneration_model ! regeneration model version
+      class(fates_patch_type), intent(in), optional :: from_patch_for_fire  ! patch to take fire info from
 
       ! initialize patch
       ! sets all values to nan, then some values to zero
-      call this%Init(num_swb, num_levsoil)
+      if (present(from_patch_for_fire)) then
+        call this%Init(num_swb, num_levsoil, from_patch_for_fire)
+      else
+        call this%Init(num_swb, num_levsoil)
+      end if
 
       ! initialize running means for patch
       call this%InitRunningMeans(current_tod, regeneration_model, num_pft)
@@ -898,6 +940,7 @@ subroutine FreeMemory(this, regeneration_model, numpft)
                  this%sabs_dir,                 &
                  this%sabs_dif,                 &
                  this%fragmentation_scaler,     &
+                 this%area_in_edgeforest_bins,  &
                  stat=istat, errmsg=smsg)
 
       ! These arrays are allocated via a call from EDCanopyStructureMod
@@ -1227,6 +1270,7 @@ subroutine Dump(this)
       write(fates_log(),*) 'pa%ncl_p              = ',this%ncl_p
       write(fates_log(),*) 'pa%total_canopy_area  = ',this%total_canopy_area
       write(fates_log(),*) 'pa%total_tree_area    = ',this%total_tree_area
+      write(fates_log(),*) 'pa%is_forest          = ',this%is_forest
       write(fates_log(),*) 'pa%total_grass_area   = ',this%total_grass_area
       write(fates_log(),*) 'pa%zstar              = ',this%zstar
       write(fates_log(),*) 'pa%gnd_alb_dif        = ',this%gnd_alb_dif(:)

fire/SFEquationsMod.F90

@@ -442,10 +442,6 @@ real(r8) function AreaBurnt(fire_size, num_ignitions, FDI)
       ! 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
       !

fire/SFFireWeatherMod.F90

@@ -12,10 +12,17 @@ module SFFireWeatherMod
     real(r8) :: effective_windspeed  ! effective wind speed, corrected for by tree/grass cover [m/min]
     integer  :: rx_flag              ! prescribed fire burn window flag [1=burn window present; 0=no burn window]
 
+    real(r8) :: temp_C  ! daily averaged temperature [deg C]
+    real(r8) :: precip  ! daily precip [mm/day]
+    real(r8) :: rh      ! daily relative humidity [%]
+    real(r8) :: wind    ! wind speed [m/min]
+
     contains 
 
     procedure(initialize_fire_weather), public, deferred :: Init
-    procedure(update_fire_weather),     public, deferred :: UpdateIndex
+    procedure(copy_fire_weather),       public, deferred :: CopyFrom
+    procedure(update_fire_weather_index), public, deferred :: UpdateIndex
+    procedure,                          public           :: UpdateFireWeatherData
     procedure,                          public           :: UpdateEffectiveWindSpeed
     procedure,                          public           :: UpdateRxfireBurnWindow
 
@@ -30,25 +37,49 @@ subroutine initialize_fire_weather(this)
 
     end subroutine initialize_fire_weather
 
-    subroutine update_fire_weather(this, temp_C, precip, rh, wind)
+    subroutine copy_fire_weather(this, from)
+
+      import :: fire_weather
+
+      class(fire_weather), intent(inout) :: this
+      class(fire_weather), intent(in)    :: from
+
+    end subroutine copy_fire_weather
+
+    subroutine update_fire_weather_index(this)
 
       use FatesConstantsMod, only : r8 => fates_r8
 
       import :: fire_weather 
 
       class(fire_weather), intent(inout) :: this
-      real(r8),            intent(in)    :: temp_C
-      real(r8),            intent(in)    :: precip
-      real(r8),            intent(in)    :: rh
-      real(r8),            intent(in)    :: wind
 
-    end subroutine update_fire_weather
+    end subroutine update_fire_weather_index
 
    end interface
 
   contains 
 
-  subroutine UpdateEffectiveWindSpeed(this, wind_speed, tree_fraction, grass_fraction,   &
+  subroutine UpdateFireWeatherData(this, temp_C, precip, rh, wind)
+    !
+    !  DESCRIPTION:
+    !  Updates fire weather variables
+
+    ! ARGUMENTS
+    class(fire_weather), intent(inout) :: this   ! fire weather class
+    real(r8),            intent(in)    :: temp_C ! daily averaged temperature [degrees C]
+    real(r8),            intent(in)    :: precip ! daily precipitation [mm]
+    real(r8),            intent(in)    :: rh     ! daily relative humidity [%]
+    real(r8),            intent(in)    :: wind   ! daily wind speed [m/min]
+
+    this%temp_C = temp_C
+    this%precip = precip
+    this%rh = rh
+    this%wind = wind
+
+  end subroutine UpdateFireWeatherData
+
+  subroutine UpdateEffectiveWindSpeed(this, tree_fraction, grass_fraction,   &
       bare_fraction)
     !
     !  DESCRIPTION:
@@ -60,25 +91,23 @@ subroutine UpdateEffectiveWindSpeed(this, wind_speed, tree_fraction, grass_fract
 
     ! ARGUMENTS
     class(fire_weather), intent(inout) :: this           ! fire weather class
-    real(r8),            intent(in)    :: wind_speed     ! wind speed [m/min]
     real(r8),            intent(in)    :: tree_fraction  ! tree fraction [0-1]
     real(r8),            intent(in)    :: grass_fraction ! grass fraction [0-1]
     real(r8),            intent(in)    :: bare_fraction  ! bare ground fraction [0-1]
 
-    this%effective_windspeed = wind_speed*(tree_fraction*wind_atten_treed +              &
+    this%effective_windspeed = this%wind * (tree_fraction*wind_atten_treed +              &
       (grass_fraction + bare_fraction)*wind_atten_grass)
 
   end subroutine UpdateEffectiveWindSpeed
 
-  subroutine UpdateRxfireBurnWindow(this, rxfire_switch, temp_C, rh, wind, temp_up,      &
+  subroutine UpdateRxfireBurnWindow(this, rxfire_switch, temp_up, &
     temp_low,rh_up, rh_low, wind_up, wind_low)
 
+    use FatesConstantsMod, only : sec_per_min
+
     ! ARGUMENTS
     class(fire_weather), intent(inout) :: this           ! fire weather class
-    real(r8),            intent(in)    :: temp_C         ! daily averaged temperature [degrees C]
     integer,             intent(in)    :: rxfire_switch  ! whether prescribed fire is turned on  
-    real(r8),            intent(in)    :: rh             ! daily relative humidity [%]
-    real(r8),            intent(in)    :: wind           ! wind speed [m/min]
     real(r8),            intent(in)    :: temp_up        ! user defined upper bound for temp when define a burn window
     real(r8),            intent(in)    :: temp_low       ! user defined lower bound for temp when define a burn window
     real(r8),            intent(in)    :: rh_up          ! user defined upper bound for relative humidity
@@ -87,21 +116,25 @@ subroutine UpdateRxfireBurnWindow(this, rxfire_switch, temp_C, rh, wind, temp_up
     real(r8),            intent(in)    :: wind_low       ! user defined lower bound for wind speed
 
     ! LOCAL VARIABLES
+    real(r8) :: wind_m_per_s  ! wind speed (m/s)
     real(r8) :: t_check  ! intermediate value derived from temp condition check
     real(r8) :: rh_check ! intermediate value derived from RH condition check
     real(r8) :: ws_check ! intermediate value derived from wind speed condition check
 
     if (rxfire_switch .eq. ifalse) return   
+
+    ! Convert from m/min to m/s
+    wind_m_per_s = this%wind / sec_per_min
 
     ! check if ambient temperature, relative humidity, and wind speed
     ! are within user defined ranges by comparing current weather
     ! condition to the lower and upper bounds defined. when within range,
     ! it should result in negative value or zero (at the boundary condition)
     ! for each check below 
 
-    t_check = (temp_C - temp_low)*(temp_C - temp_up)
-    rh_check = (rh - rh_low)*(rh - rh_up)
-    ws_check = (wind - wind_low)*(wind - wind_up)
+    t_check = (this%temp_C - temp_low)*(this%temp_C - temp_up)
+    rh_check = (this%rh - rh_low)*(this%rh - rh_up)
+    ws_check = (wind_m_per_s - wind_low)*(wind_m_per_s - wind_up)
 
     if (t_check <= 0.0_r8 .and. rh_check <= 0.0_r8 .and. ws_check <= 0.0_r8) then
       this%rx_flag = 1