Grazing

biogeochem/EDPhysiologyMod.F90

@@ -2785,7 +2785,9 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
     ! and turnover in dying trees.
     !
     ! !USES:
-
+    use EDParamsMod           , only : landuse_grazing_carbon_use_eff
+    use EDParamsMod           , only : landuse_grazing_nitrogen_use_eff
+    use EDParamsMod           , only : landuse_grazing_phosphorus_use_eff
     !
     ! !ARGUMENTS
     type(ed_site_type), intent(inout), target :: currentSite
@@ -2832,6 +2834,8 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
     integer  :: numlevsoil        ! Actual number of soil layers
 
     real(r8) :: SF_val_CWD_frac_adj(4) !SF_val_CWD_frac adjusted based on cohort dbh
+    real(r8) :: leaf_herbivory
+    real(r8) :: herbivory_element_use_efficiency
     !----------------------------------------------------------------------
 
     ! -----------------------------------------------------------------------------------
@@ -2842,6 +2846,15 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
 
     element_id = litt%element_id
 
+    select case(element_id)
+       case (carbon12_element)
+          herbivory_element_use_efficiency = landuse_grazing_carbon_use_eff
+       case (nitrogen_element)
+          herbivory_element_use_efficiency = landuse_grazing_nitrogen_use_eff
+       case (phosphorus_element)
+          herbivory_element_use_efficiency = landuse_grazing_phosphorus_use_eff
+    end select
+
     ! Object tracking flux diagnostics for each element
     flux_diags => currentSite%flux_diags(element_pos(element_id))
 
@@ -2863,6 +2876,8 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
        fnrt_m          = currentCohort%prt%GetState(fnrt_organ,element_id)
        repro_m         = currentCohort%prt%GetState(repro_organ,element_id)
 
+       leaf_herbivory  = currentCohort%prt%GetHerbivory(leaf_organ, element_id)
+
        if (prt_params%woody(currentCohort%pft) == itrue) then
           ! Assumption: for woody plants fluxes from deadwood and sapwood go together in CWD pool
           leaf_m_turnover   = currentCohort%prt%GetTurnover(leaf_organ,element_id)
@@ -2907,7 +2922,9 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
        do dcmpy=1,ndcmpy
           dcmpy_frac = GetDecompyFrac(pft,leaf_organ,dcmpy)
           litt%leaf_fines_in(dcmpy) = litt%leaf_fines_in(dcmpy) + &
-               (leaf_m_turnover+repro_m_turnover) * plant_dens * dcmpy_frac
+               (leaf_m_turnover+repro_m_turnover + &
+               leaf_herbivory * herbivory_element_use_efficiency) * &
+               plant_dens * dcmpy_frac
 
           dcmpy_frac = GetDecompyFrac(pft,fnrt_organ,dcmpy)
           do ilyr = 1, numlevsoil
@@ -2920,6 +2937,11 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
             flux_diags%root_litter_input(pft) +  &
             (fnrt_m_turnover + store_m_turnover ) * currentCohort%n
 
+       ! send the part of the herbivory flux that doesn't go to litter to the atmosphere
+
+       currentSite%mass_balance(element_id)%herbivory_flux_out = &
+            currentSite%mass_balance(element_id)%herbivory_flux_out + &
+            leaf_herbivory * (1._r8 - herbivory_element_use_efficiency) * currentCohort%n
 
        ! Assumption: turnover from deadwood and sapwood are lumped together in CWD pool
 

biogeochem/FatesLandUseChangeMod.F90

@@ -33,6 +33,7 @@ module FatesLandUseChangeMod
   public :: GetLUHStatedata
   public :: GetInitLanduseTransitionRates
   public :: GetInitLanduseHarvestRate
+  public :: fates_grazing
 
   ! module data
   integer, parameter :: max_luh2_types_per_fates_lu_type = 5
@@ -424,4 +425,44 @@ subroutine GetInitLanduseTransitionRates(bc_in, min_allowed_landuse_fraction, &
 
   end subroutine GetInitLanduseTransitionRates
 
+  !----------------------------------------------------------------------------------------------------
+
+  subroutine fates_grazing(prt, ft, land_use_label, height)
+
+    use PRTGenericMod,    only : leaf_organ
+    use PRTGenericMod,    only : prt_vartypes
+    use PRTLossFluxesMod, only : PRTHerbivoryLosses
+    use EDParamsMod     , only : landuse_grazing_rate
+    use EDParamsMod     , only : landuse_grazing_maxheight
+    use EDPftvarcon     , only : EDPftvarcon_inst
+    use PRTParametersMod, only : prt_params
+    use FatesAllometryMod,only : CrownDepth
+
+    ! apply grazing and browsing to plants as a function of PFT, height (for woody plants), and land use label.
+
+    class(prt_vartypes), intent(inout), pointer :: prt
+    integer, intent(in)  :: ft
+    integer, intent(in)  :: land_use_label
+    real(r8), intent(in) :: height
+
+    real(r8) :: grazing_rate    ! rate of grazing (or browsing) of leaf tissue [day -1]
+    real(r8) :: crown_depth
+
+    grazing_rate = landuse_grazing_rate(land_use_label) * EDPftvarcon_inst%landuse_grazing_palatability(ft)
+
+    if ( grazing_rate .gt. 0._r8) then
+       if (prt_params%woody(ft)) then
+
+          call CrownDepth(height,ft,crown_depth)
+
+          grazing_rate = grazing_rate * &
+               max(0._r8, min(1._r8, &
+               (landuse_grazing_maxheight - (height - crown_depth )) / crown_depth))
+       endif
+
+       call PRTHerbivoryLosses(prt, leaf_organ, grazing_rate)
+    end if
+
+  end subroutine fates_grazing
+
 end module FatesLandUseChangeMod

main/EDMainMod.F90

@@ -104,6 +104,7 @@ module EDMainMod
   use PRTLossFluxesMod,       only : PRTReproRelease
   use EDPftvarcon,            only : EDPftvarcon_inst
   use FatesHistoryInterfaceMod, only : fates_hist
+  use FatesLandUseChangeMod,  only: fates_grazing
 
   ! CIME Globals
   use shr_log_mod         , only : errMsg => shr_log_errMsg
@@ -524,6 +525,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
              bc_out%ar_site = bc_out%ar_site + currentCohort%resp_acc_hold * & 
                   AREA_INV * currentCohort%n / hlm_days_per_year / sec_per_day
 
+             ! allow herbivores to graze
+             call fates_grazing(currentCohort%prt, ft, currentPatch%land_use_label, currentCohort%height)
+
              ! Conduct Maintenance Turnover (parteh)
              if(debug) call currentCohort%prt%CheckMassConservation(ft,3)
              if(any(currentSite%dstatus(ft) == [phen_dstat_moiston,phen_dstat_timeon])) then
@@ -918,7 +922,8 @@ subroutine TotalBalanceCheck (currentSite, call_index )
                   site_mass%seed_out + &
                   site_mass%flux_generic_out + &
                   site_mass%frag_out + &
-                  site_mass%aresp_acc
+                  site_mass%aresp_acc + &
+                  site_mass%herbivory_flux_out
 
        net_flux        = flux_in - flux_out
        error           = abs(net_flux - change_in_stock)
@@ -951,6 +956,7 @@ subroutine TotalBalanceCheck (currentSite, call_index )
           write(fates_log(),*) 'flux_generic_out: ',site_mass%flux_generic_out
           write(fates_log(),*) 'frag_out: ',site_mass%frag_out
           write(fates_log(),*) 'aresp_acc: ',site_mass%aresp_acc
+          write(fates_log(),*) 'herbivory_flux_out: ',site_mass%herbivory_flux_out
           write(fates_log(),*) 'error=net_flux-dstock:', error
           write(fates_log(),*) 'biomass', biomass_stock
           write(fates_log(),*) 'litter',litter_stock

main/EDParamsMod.F90

@@ -117,6 +117,13 @@ module EDParamsMod
    real(r8),protected,public  :: q10_mr     ! Q10 for respiration rate (for soil fragmenation and plant respiration)    (unitless)
    real(r8),protected,public  :: q10_froz   ! Q10 for frozen-soil respiration rates (for soil fragmentation)            (unitless)
 
+   ! grazing parameters
+   real(r8),protected,public :: landuse_grazing_carbon_use_eff
+   real(r8),protected,public :: landuse_grazing_maxheight
+   real(r8),protected,public :: landuse_grazing_nitrogen_use_eff
+   real(r8),protected,public :: landuse_grazing_phosphorus_use_eff
+   real(r8),protected,public :: landuse_grazing_rate(n_landuse_cats)
+
    ! Unassociated pft dimensioned free parameter that developers can use for testing arbitrary new hypotheses
    ! (THIS PARAMETER IS UNUSED, FEEL FREE TO USE IT FOR WHATEVER PURPOSE YOU LIKE. WE CAN
    !  HELP MIGRATE YOUR USAGE OF THE PARMETER TO A PERMANENT HOME LATER)
@@ -186,6 +193,7 @@ module EDParamsMod
    character(len=param_string_length),parameter,public :: ED_name_maxpatches_by_landuse = "fates_maxpatches_by_landuse"
    character(len=param_string_length),parameter,public :: ED_name_max_nocomp_pfts_by_landuse = "fates_max_nocomp_pfts_by_landuse"
 
+
    ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE)
    ! ----------------------------------------------------------------------------------------------
    real(r8),protected,public :: hydr_kmax_rsurf1         !  maximum conducitivity for unit root surface 
@@ -285,6 +293,13 @@ module EDParamsMod
                                                     ! leftovers will be left onsite as large CWD
    character(len=param_string_length),parameter,public :: logging_name_export_frac ="fates_landuse_logging_export_frac"   
 
+   ! grazing-related parameters
+   character(len=param_string_length),parameter,public :: name_landuse_grazing_rate               = "fates_landuse_grazing_rate"
+   character(len=param_string_length),parameter,public :: name_landuse_grazing_carbon_use_eff     = "fates_landuse_grazing_carbon_use_eff"
+   character(len=param_string_length),parameter,public :: name_landuse_grazing_maxheight          = "fates_landuse_grazing_maxheight"
+   character(len=param_string_length),parameter,public :: name_landuse_grazing_nitrogen_use_eff   = "fates_landuse_grazing_nitrogen_use_eff"
+   character(len=param_string_length),parameter,public :: name_landuse_grazing_phosphorus_use_eff = "fates_landuse_grazing_phosphorus_use_eff"
+
    real(r8),protected,public :: eca_plant_escalar  ! scaling factor for plant fine root biomass to 
                                                ! calculate nutrient carrier enzyme abundance (ECA)
 
@@ -368,6 +383,12 @@ subroutine FatesParamsInit()
     dev_arbitrary                         = nan
     damage_event_code                     = -9
     damage_canopy_layer_code              = -9
+    landuse_grazing_carbon_use_eff        = nan
+    landuse_grazing_nitrogen_use_eff      = nan
+    landuse_grazing_phosphorus_use_eff    = nan
+    landuse_grazing_maxheight             = nan
+    landuse_grazing_rate(:)               = nan
+
   end subroutine FatesParamsInit
 
   !-----------------------------------------------------------------------
@@ -571,7 +592,19 @@ 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=name_landuse_grazing_carbon_use_eff, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=name_landuse_grazing_maxheight, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=name_landuse_grazing_nitrogen_use_eff, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=name_landuse_grazing_phosphorus_use_eff, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
     ! non-scalar parameters
 
     call fates_params%RegisterParameter(name=ED_name_hydr_htftype_node, dimension_shape=dimension_shape_1d, &
@@ -607,6 +640,8 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_max_nocomp_pfts_by_landuse, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names_landuse)
 
+    call fates_params%RegisterParameter(name=name_landuse_grazing_rate, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names_landuse)
   end subroutine FatesRegisterParams
 
 
@@ -625,6 +660,7 @@ subroutine FatesReceiveParams(fates_params)
     real(r8), allocatable :: tmp_vector_by_landuse1(:)  ! local real vector for changing type on read
     real(r8), allocatable :: tmp_vector_by_landuse2(:)  ! local real vector for changing type on read
     real(r8), allocatable :: tmp_vector_by_landuse3(:)  ! local real vector for changing type on read
+    real(r8), allocatable :: tmp_vector_by_landuse4(:)  ! local real vector for changing type on read
 
     call fates_params%RetrieveParameter(name=ED_name_photo_temp_acclim_timescale, &
          data=photo_temp_acclim_timescale)
@@ -860,6 +896,25 @@ subroutine FatesReceiveParams(fates_params)
     hydr_htftype_node(:) = nint(hydr_htftype_real(:))
     deallocate(hydr_htftype_real)
 
+    call fates_params%RetrieveParameter(name=name_landuse_grazing_carbon_use_eff, &
+         data=landuse_grazing_carbon_use_eff)
+
+    call fates_params%RetrieveParameter(name=name_landuse_grazing_nitrogen_use_eff, &
+         data=landuse_grazing_nitrogen_use_eff)
+
+    call fates_params%RetrieveParameter(name=name_landuse_grazing_phosphorus_use_eff, &
+         data=landuse_grazing_phosphorus_use_eff)
+
+    call fates_params%RetrieveParameter(name=name_landuse_grazing_maxheight, &
+         data=landuse_grazing_maxheight)
+
+    call fates_params%RetrieveParameterAllocate(name=name_landuse_grazing_rate, &
+         data=tmp_vector_by_landuse4)
+
+    landuse_grazing_rate(:) = tmp_vector_by_landuse4(:)
+
+    deallocate(tmp_vector_by_landuse4)
+
   end subroutine FatesReceiveParams
 
   ! =====================================================================================
@@ -930,6 +985,11 @@ 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) '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
+        write(fates_log(),fmt0) 'name_landuse_grazing_maxheight = ', name_landuse_grazing_maxheight
+        write(fates_log(),fmt0) 'name_landuse_grazing_rate(:) = ', name_landuse_grazing_rate(:)
         write(fates_log(),*) '------------------------------------------------------'
 
      end if

main/EDPftvarcon.F90

@@ -286,6 +286,9 @@ module EDPftvarcon
      real(r8), allocatable :: landusechange_frac_exported(:)  ! fraction of land use change-generated wood material that is exported to wood product (the remainder is either burned or goes to litter)
      real(r8), allocatable :: landusechange_pprod10(:)        ! fraction of land use change wood product that goes to 10-year product pool (remainder goes to 100-year pool)
 
+     ! Grazing
+     real(r8), allocatable :: landuse_grazing_palatability(:) ! Relative intensity of leaf grazing/browsing per PFT (unitless 0-1)
+
    contains
      procedure, public :: Init => EDpftconInit
      procedure, public :: Register
@@ -818,7 +821,11 @@ subroutine Register_PFT(this, fates_params)
     name = 'fates_landuse_luc_pprod10'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
-
+
+    name = 'fates_landuse_grazing_palatability'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_dev_arbitrary_pft'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -1303,6 +1310,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%landusechange_pprod10)
 
+    name = 'fates_landuse_grazing_palatability'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%landuse_grazing_palatability)
+
   end subroutine Receive_PFT
 
   !-----------------------------------------------------------------------

main/EDTypesMod.F90

@@ -235,6 +235,8 @@ module EDTypesMod
                                        ! due to re-sizing patches when area math starts to lose
                                        ! precision
 
+     real(r8) :: herbivory_flux_out    ! loss of element due to grazing (and/or browsing) by herbivores
+
    contains
 
      procedure :: ZeroMassBalState
@@ -531,6 +533,7 @@ subroutine ZeroMassBalFlux(this)
       this%flux_generic_in   = 0._r8
       this%flux_generic_out  = 0._r8
       this%patch_resize_err  = 0._r8
+      this%herbivory_flux_out= 0._r8
 
       return
   end subroutine ZeroMassBalFlux

main/FatesHistoryInterfaceMod.F90

@@ -470,6 +470,7 @@ module FatesHistoryInterfaceMod
   integer :: ih_npp_agsw_si_scpf
   integer :: ih_npp_agdw_si_scpf
   integer :: ih_npp_stor_si_scpf
+  integer :: ih_grazing_si
 
   integer :: ih_mortality_canopy_si_scpf
   integer :: ih_mortality_understory_si_scpf
@@ -2388,6 +2389,7 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in)
     real(r8) :: fnrt_m_net_alloc   ! mass allocated to fine-root [kg/yr]
     real(r8) :: struct_m_net_alloc ! mass allocated to structure [kg/yr]
     real(r8) :: repro_m_net_alloc  ! mass allocated to reproduction [kg/yr]
+    real(r8) :: leaf_herbivory     ! mass of leaves eaten by herbivores [kg/yr]
     real(r8) :: n_perm2            ! abundance per m2
     real(r8) :: area_frac  ! Fraction of area for this patch
 
@@ -2440,6 +2442,7 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in)
          hio_npp_froot_si        => this%hvars(ih_npp_froot_si)%r81d, &
          hio_npp_croot_si        => this%hvars(ih_npp_croot_si)%r81d, &
          hio_npp_stor_si         => this%hvars(ih_npp_stor_si)%r81d, &
+         hio_grazing_si          => this%hvars(ih_grazing_si)%r81d, &
          hio_canopy_mortality_crownarea_si     => this%hvars(ih_canopy_mortality_crownarea_si)%r81d, &
          hio_ustory_mortality_crownarea_si => this%hvars(ih_understory_mortality_crownarea_si)%r81d, &
          hio_fire_c_to_atm_si  => this%hvars(ih_fire_c_to_atm_si)%r81d, &
@@ -2891,6 +2894,9 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in)
                   hio_npp_stor_si(io_si) = hio_npp_stor_si(io_si) +               &
                        store_m_net_alloc * n_perm2 / days_per_year / sec_per_day
 
+                  leaf_herbivory   = ccohort%prt%GetHerbivory(leaf_organ, carbon12_element) * days_per_year  !cdkcdk
+                  hio_grazing_si(io_si) = hio_grazing_si(io_si) + leaf_herbivory * n_perm2 / days_per_year / sec_per_day
+
                   ! Woody State Variables (basal area growth increment)
                   if ( prt_params%woody(ft) == itrue) then
 
@@ -6753,6 +6759,12 @@ subroutine define_history_vars(this, initialize_variables)
             upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables,                &
             index = ih_npp_stor_si)
 
+       call this%set_history_var(vname='FATES_GRAZING', units='kg m-2 s-1',    &
+            long='grazing by herbivores of leaves in kg carbon per m2 per second',          &
+            use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM',    &
+            upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables,                &
+            index = ih_grazing_si)
+
        hydro_active_if: if(hlm_use_planthydro.eq.itrue) then
           call this%set_history_var(vname='FATES_VEGH2O_DEAD', units = 'kg m-2',  &
                long='cumulative water stored in dead biomass due to mortality',  &