Diagnostically track forest edge

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/FatesPatchMod.F90

@@ -32,6 +32,7 @@ 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 shr_infnan_mod,         only : nan => shr_infnan_nan, assignment(=)
   use shr_log_mod,            only : errMsg => shr_log_errMsg
 
@@ -74,6 +75,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]
@@ -279,6 +286,7 @@ 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))
 
       ! initialize all values to nan
       call this%NanValues()
@@ -453,6 +461,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
@@ -898,6 +910,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 +1240,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/SFMainMod.F90

@@ -81,7 +81,7 @@ subroutine UpdateFireWeather(currentSite, bc_in)
 
     use FatesConstantsMod,      only : tfrz => t_water_freeze_k_1atm
     use FatesConstantsMod,      only : sec_per_day, sec_per_min
-    use EDTypesMod,             only : CalculateTreeGrassAreaSite
+    use FatesEdgeForestMod,     only : CalculateTreeGrassAreaSite
     use FatesInterfaceTypesMod, only : hlm_use_managed_fire
     use SFParamsMod,            only : SF_val_rxfire_tpup, SF_val_rxfire_tplw, SF_val_rxfire_rhup, &
                                        SF_val_rxfire_rhlw, SF_val_rxfire_wdup, SF_val_rxfire_wdlw

main/CMakeLists.txt

@@ -6,6 +6,9 @@ list(APPEND clm_sources
   EDTypesMod.F90
   EDPftvarcon.F90
   FatesConstantsMod.F90
+  FatesEcotypesMod.F90
+  FatesEdgeForestMod.F90
+  FatesEdgeForestParamsMod.F90
   FatesHydraulicsMemMod.F90
   FatesParametersInterface.F90
   FatesUtilsMod.F90
@@ -27,6 +30,10 @@ list(APPEND fates_sources
   FatesSizeAgeTypeIndicesMod.F90
   FatesIntegratorsMod.F90
   FatesUtilsMod.F90
-  FatesSynchronizedParamsMod.F90)
+  FatesSynchronizedParamsMod.F90
+  FatesEcotypesMod.F90
+  FatesEdgeForestMod.F90
+  FatesEdgeForestParamsMod.F90
+  )
 
 sourcelist_to_parent(fates_sources)

main/EDParamsMod.F90

@@ -213,7 +213,9 @@ module EDParamsMod
    integer, protected, public :: max_cohort_per_patch
    character(len=param_string_length), parameter, public :: maxcohort_name = "fates_maxcohort"
 
-
+   ! Ecotypes parameters
+   real(r8),protected,public :: forest_tree_fraction_threshold  ! Tree fraction above which a patch is "forest"
+   character(len=param_string_length),parameter,public :: forest_tree_fraction_threshold_name = "fates_forest_tree_fraction_threshold"
 
    ! Logging Control Parameters (ONLY RELEVANT WHEN USE_FATES_LOGGING = TRUE)
    ! ----------------------------------------------------------------------------------------------
@@ -266,6 +268,7 @@ module EDParamsMod
    character(len=param_string_length),parameter,public :: eca_name_plant_escalar = "fates_cnp_eca_plant_escalar"
 
    public :: FatesParamsInit
+   public :: FatesParamsInitForFactory
    public :: FatesRegisterParams
    public :: FatesReceiveParams
    public :: FatesReportParams
@@ -342,6 +345,7 @@ subroutine FatesParamsInit()
     dev_arbitrary                         = nan
     damage_event_code                     = -9
     damage_canopy_layer_code              = -9
+    forest_tree_fraction_threshold        = nan
     landuse_grazing_carbon_use_eff        = nan
     landuse_grazing_nitrogen_use_eff      = nan
     landuse_grazing_phosphorus_use_eff    = nan
@@ -350,6 +354,21 @@ subroutine FatesParamsInit()
 
   end subroutine FatesParamsInit
 
+  !-----------------------------------------------------------------------
+
+  subroutine FatesParamsInitForFactory()
+    ! Initialize some parameters that are needed for unit-testing factories
+
+    allocate(ED_val_history_ageclass_bin_edges(7))
+    ED_val_history_ageclass_bin_edges = [0, 1, 2, 5, 10, 20, 50]
+
+    allocate(ED_val_history_sizeclass_bin_edges(13))
+    ED_val_history_sizeclass_bin_edges = [0, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100]
+
+    allocate(ED_val_history_coageclass_bin_edges(13))
+    ED_val_history_coageclass_bin_edges = [0, 5]
+  end subroutine FatesParamsInitForFactory
+
   !-----------------------------------------------------------------------
   subroutine FatesRegisterParams(fates_params)
     ! Register the parameters we want the host to provide, and
@@ -518,6 +537,9 @@ subroutine FatesRegisterParams(fates_params)
 
     call fates_params%RegisterParameter(name=damage_name_canopy_layer_code, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=forest_tree_fraction_threshold_name, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
 
     call fates_params%RegisterParameter(name=name_landuse_grazing_carbon_use_eff, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
@@ -739,6 +761,9 @@ subroutine FatesReceiveParams(fates_params)
          data=tmpreal)
     damage_canopy_layer_code = nint(tmpreal)
 
+    call fates_params%RetrieveParameter(name=forest_tree_fraction_threshold_name, &
+         data=forest_tree_fraction_threshold)
+
     ! parameters that are arrays of size defined within the params file and thus need allocating as well
     call fates_params%RetrieveParameterAllocate(name=ED_name_history_sizeclass_bin_edges, &
           data=ED_val_history_sizeclass_bin_edges)
@@ -864,6 +889,7 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),'(a,L2)') 'active_crown_fire = ',active_crown_fire
         write(fates_log(),fmt0) 'damage_event_code = ',damage_event_code
         write(fates_log(),fmt0) 'damage_canopy_layer_code = ', damage_canopy_layer_code
+        write(fates_log(),fmt0) 'forest_tree_fraction_threshold = ', forest_tree_fraction_threshold
 	write(fates_log(),fmt0) 'landuse_grazing_carbon_use_eff = ', landuse_grazing_carbon_use_eff
         write(fates_log(),fmt0) 'name_landuse_grazing_nitrogen_use_eff = ', name_landuse_grazing_nitrogen_use_eff
         write(fates_log(),fmt0) 'name_landuse_grazing_phosphorus_use_eff = ', name_landuse_grazing_phosphorus_use_eff

main/EDTypesMod.F90

@@ -605,12 +605,12 @@ module EDTypesMod
 
        procedure, public :: get_current_landuse_statevector
        procedure, public :: get_secondary_young_fraction
+       procedure, public :: GetNumberOfPatches
 
   end type ed_site_type
 
   ! Make public necessary subroutines and functions
   public :: dump_site
-  public :: CalculateTreeGrassAreaSite
   public :: set_patchno
 
 contains
@@ -737,39 +737,6 @@ end subroutine ZeroMassBalFlux
 
   ! =====================================================================================
 
-  subroutine CalculateTreeGrassAreaSite(csite, tree_fraction, grass_fraction, bare_fraction)
-    !
-    !  DESCRIPTION:
-    !  Calculates total grass, tree, and bare fractions for a site
-
-    ! ARGUMENTS:
-    type(ed_site_type), intent(inout) :: csite          ! site object
-    real(r8),           intent(out)   :: tree_fraction  ! total site tree fraction
-    real(r8),           intent(out)   :: grass_fraction ! total site grass fraction
-    real(r8),           intent(out)   :: bare_fraction  ! total site bare fraction
-
-    ! LOCALS:
-    type(fates_patch_type), pointer :: currentPatch ! patch object
-
-    tree_fraction = 0.0_r8
-    grass_fraction = 0.0_r8
-
-    currentPatch => csite%oldest_patch
-    do while(associated(currentPatch))
-      if (currentPatch%nocomp_pft_label /= nocomp_bareground) then
-        call currentPatch%UpdateTreeGrassArea()
-        tree_fraction = tree_fraction + currentPatch%total_tree_area/AREA
-        grass_fraction = grass_fraction + currentPatch%total_grass_area/AREA
-      end if 
-      currentPatch => currentPatch%younger
-    end do
-
-    ! if cover > 1.0, grasses are under the trees
-    grass_fraction = min(grass_fraction, 1.0_r8 - tree_fraction)
-    bare_fraction = 1.0_r8 - tree_fraction - grass_fraction
-
-  end subroutine CalculateTreeGrassAreaSite
-
   !---------------------------------------------------------------------------------------
 
   subroutine dump_site(csite) 
@@ -865,4 +832,26 @@ function get_secondary_young_fraction(this) result(secondary_young_fraction)
 
    end function get_secondary_young_fraction
 
+   function GetNumberOfPatches(this) result(n_patches)
+     ! DESCRIPTION
+     ! Returns number of patches at site
+     !
+     ! ARGUMENTS:
+     class(ed_site_type) :: this
+     !
+     ! RETURN VALUE:
+     integer :: n_patches
+     !
+     ! LOCAL VARIABLES:
+     type(fates_patch_type), pointer :: currentPatch
+
+     n_patches = 0
+     currentPatch => this%youngest_patch
+     do while(associated(currentPatch))
+        n_patches = n_patches + 1
+        currentPatch => currentPatch%older
+     enddo
+
+  end function GetNumberOfPatches
+
 end module EDTypesMod

main/FatesEcotypesMod.F90

@@ -0,0 +1,81 @@
+module FatesEcotypesMod
+
+  use FatesConstantsMod, only : r8 => fates_r8
+  use EDTypesMod, only : ed_site_type
+  use FatesPatchMod, only : fates_patch_type
+
+  implicit none
+  private  ! By default everything is private
+
+  ! Make public necessary subroutines and functions
+  public :: IsPatchForest
+  ! For unit testing
+  public :: DoesPatchHaveForest_TreeCover
+  public :: DoesPatchHaveForest_GrassBiomass
+
+contains
+
+  ! =====================================================================================
+
+  function DoesPatchHaveForest_TreeCover(patchptr, forest_tree_fraction_threshold)
+  ! DESCRIPTION:
+  ! Return boolean: Is this patch "forest"?
+  !
+  ! ARGUMENTS:
+  type(fates_patch_type), intent(in), pointer :: patchptr  ! pointer to patch object
+  real(r8), intent(in) :: forest_tree_fraction_threshold ! Tree fraction above which a patch is "forest"
+  !
+  ! RETURN VALUE
+  logical :: DoesPatchHaveForest_TreeCover
+  !
+  ! LOCAL VARIABLES
+  real(r8) :: tree_fraction = 0._r8
+
+  if (patchptr%area > 0._r8) then
+      tree_fraction = patchptr%total_tree_area / patchptr%area
+  else
+      tree_fraction = 0._r8
+  end if
+
+  DoesPatchHaveForest_TreeCover = tree_fraction > forest_tree_fraction_threshold
+
+  end function DoesPatchHaveForest_TreeCover
+
+
+  function DoesPatchHaveForest_GrassBiomass(patchptr, grass_biomass_threshold)
+  ! DESCRIPTION:
+  ! Return boolean: Does this patch have grass biomass above a threshold?
+  !
+  ! ARGUMENTS:
+  type(fates_patch_type), intent(in), pointer :: patchptr  ! pointer to patch object
+  real(r8), intent(in) :: grass_biomass_threshold ! Live grass biomass (kgC/m2) above which a patch is considered to "have grass"
+  !
+  ! RETURN VALUE
+  logical :: DoesPatchHaveForest_GrassBiomass
+
+  DoesPatchHaveForest_GrassBiomass = patchptr%livegrass > grass_biomass_threshold
+
+  end function DoesPatchHaveForest_GrassBiomass
+
+
+  function IsPatchForest(patchptr, forest_tree_fraction_threshold, grass_biomass_threshold)
+  ! DESCRIPTION:
+  ! Return boolean: Is this patch forest according to tree cover and, optionally, grass biomass?
+  !
+  ! ARGUMENTS:
+  type(fates_patch_type), intent(in), pointer :: patchptr  ! pointer to patch object
+  real(r8), intent(in) :: forest_tree_fraction_threshold ! Tree fraction above which a patch is "forest"
+  real(r8), intent(in), optional :: grass_biomass_threshold ! Live grass biomass (kgC/m2) above which a patch is considered to "have grass"
+  !
+  ! RETURN VALUE
+  logical :: IsPatchForest
+
+  IsPatchForest = DoesPatchHaveForest_TreeCover(patchptr, forest_tree_fraction_threshold)
+  if (IsPatchForest .and. present(grass_biomass_threshold)) then
+     IsPatchForest = .not. DoesPatchHaveForest_GrassBiomass(patchptr, grass_biomass_threshold)
+  end if
+
+  end function IsPatchForest
+
+
+end module FatesEcotypesMod