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 : is_patch_forest
+    use EDParamsMod          , only : forest_tree_fraction_threshold
+    use FatesEdgeForestMod   , only : calculate_edgeforest_area
 
     ! !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 = is_patch_forest(currentPatch, forest_tree_fraction_threshold)
+
           currentPatch => currentPatch%younger
        end do !patch loop
 
+       if (hlm_use_edge_forest == itrue) then
+          call calculate_edgeforest_area(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]
@@ -269,6 +276,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()
@@ -443,6 +451,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
@@ -878,6 +890,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
@@ -1207,6 +1220,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

@@ -78,7 +78,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
 
     ! ARGUMENTS:
     type(ed_site_type), intent(inout), target :: currentSite

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

@@ -58,6 +58,7 @@ module EDInitMod
   use FatesInterfaceTypesMod         , only : hlm_use_inventory_init
   use FatesInterfaceTypesMod         , only : hlm_use_fixed_biogeog
   use FatesInterfaceTypesMod         , only : hlm_use_tree_damage
+  use FatesInterfaceTypesMod         , only : hlm_use_edge_forest
   use FatesInterfaceTypesMod         , only : hlm_use_sp
   use FatesInterfaceTypesMod         , only : hlm_use_luh
   use FatesInterfaceTypesMod         , only : numpft

main/EDParamsMod.F90

@@ -212,7 +212,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)
    ! ----------------------------------------------------------------------------------------------
@@ -341,6 +343,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
@@ -517,6 +520,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)
@@ -738,6 +744,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)
@@ -863,6 +872,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

@@ -574,7 +574,6 @@ module EDTypesMod
 
   ! Make public necessary subroutines and functions
   public :: dump_site
-  public :: CalculateTreeGrassAreaSite
   public :: set_patchno
 
 contains
@@ -701,39 +700,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) 

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 :: is_patch_forest
+  ! For unit testing
+  public :: does_patch_have_forest_tcthresh
+  public :: does_patch_have_grass_bmthresh
+
+contains
+
+  ! =====================================================================================
+
+  function does_patch_have_forest_tcthresh(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 :: does_patch_have_forest_tcthresh
+  !
+  ! 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
+
+  does_patch_have_forest_tcthresh = tree_fraction > forest_tree_fraction_threshold
+
+  end function does_patch_have_forest_tcthresh
+
+
+  function does_patch_have_grass_bmthresh(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 :: does_patch_have_grass_bmthresh
+
+  does_patch_have_grass_bmthresh = patchptr%livegrass > grass_biomass_threshold
+
+  end function does_patch_have_grass_bmthresh
+
+
+  function is_patch_forest(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 :: is_patch_forest
+
+  is_patch_forest = does_patch_have_forest_tcthresh(patchptr, forest_tree_fraction_threshold)
+  if (is_patch_forest .and. present(grass_biomass_threshold)) then
+     is_patch_forest = .not. does_patch_have_grass_bmthresh(patchptr, grass_biomass_threshold)
+  end if
+
+  end function is_patch_forest
+
+
+end module FatesEcotypesMod