Add grazing of storage organ

biogeochem/EDPhysiologyMod.F90

@@ -2878,7 +2878,7 @@ 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  ! leaf that is eaten by grazers [kg]
+    real(r8) :: herbivory  ! leaf that is eaten by grazers [kg]
     real(r8) :: herbivory_element_use_efficiency   ! fraction of grazed biomass that is returned to litter pool versus atmosphere
     !----------------------------------------------------------------------
 
@@ -2909,6 +2909,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
 
     currentCohort => currentPatch%shortest
     do while(associated(currentCohort))
+       herbivory = 0.0_r8
 
        pft = currentCohort%pft
        call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, &
@@ -2924,7 +2925,9 @@ 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)
+       herbivory  = currentCohort%prt%GetHerbivory(leaf_organ, element_id)
+       ! we also graze storage now.
+       herbivory  = herbivory + currentCohort%prt%GetHerbivory(store_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
@@ -2971,7 +2974,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
           dcmpy_frac = GetDecompyFrac(pft,leaf_organ,dcmpy)
           litt%leaf_fines_in(dcmpy) = litt%leaf_fines_in(dcmpy) + &
                (leaf_m_turnover+repro_m_turnover + &
-               leaf_herbivory * herbivory_element_use_efficiency) * &
+               herbivory * herbivory_element_use_efficiency) * &
                plant_dens * dcmpy_frac
 
           dcmpy_frac = GetDecompyFrac(pft,fnrt_organ,dcmpy)
@@ -2989,7 +2992,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
 
        site_mass%herbivory_flux_out = &
             site_mass%herbivory_flux_out + &
-            leaf_herbivory * (1._r8 - herbivory_element_use_efficiency) * currentCohort%n
+            herbivory * (1._r8 - herbivory_element_use_efficiency) * currentCohort%n
 
 
        ! Assumption: turnover from deadwood and sapwood are lumped together in CWD pool

biogeochem/FatesLandUseChangeMod.F90

@@ -451,9 +451,9 @@ end subroutine GetInitLanduseTransitionRates
 
   !----------------------------------------------------------------------------------------------------
 
-  subroutine FatesGrazing(prt, ft, land_use_label, height,npp_acc, lai)
+  subroutine FatesGrazing(prt, ft, land_use_label, height,npp_acc, lai, leaf_store_frac)
 
-    use PRTGenericMod,    only : leaf_organ
+    use PRTGenericMod,    only : leaf_organ, store_organ
     use PRTGenericMod,    only : prt_vartypes
     use PRTLossFluxesMod, only : PRTHerbivoryLosses
     use EDParamsMod     , only : landuse_grazing_rate
@@ -470,6 +470,7 @@ subroutine FatesGrazing(prt, ft, land_use_label, height,npp_acc, lai)
     real(r8), intent(in) :: height
     real(r8), intent(in) :: npp_acc
     real(r8), intent(in) :: lai
+    real(r8), intent(in) :: leaf_store_frac ! fraction of leaf storage to total
 
     real(r8) :: grazing_rate    ! rate of grazing (or browsing) of leaf tissue [day -1]
     real(r8) :: crown_depth
@@ -493,7 +494,8 @@ subroutine FatesGrazing(prt, ft, land_use_label, height,npp_acc, lai)
        endif
 
        call PRTHerbivoryLosses(prt, leaf_organ, grazing_rate)
-
+       call PRTHerbivoryLosses(prt, store_organ, grazing_rate * leaf_store_frac)
+
     end if
 
   end subroutine FatesGrazing

main/EDMainMod.F90

@@ -348,6 +348,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
     use FatesInterfaceTypesMod, only : hlm_use_cohort_age_tracking
     use FatesConstantsMod, only : itrue
     use FatesConstantsMod     , only : nearzero
+    use FatesConstantsMod        , only : rsnbl_math_prec
     use EDCanopyStructureMod  , only : canopy_structure
 
 
@@ -401,6 +402,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
     real(r8) :: repro_c
     real(r8) :: total_c
     real(r8) :: store_c
+    real(r8) :: alive_c
+    real(r8) :: leaf_c_frac
 
     real(r8) :: cc_leaf_c
     real(r8) :: cc_fnrt_c
@@ -554,9 +557,16 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
 
 
              ! allow herbivores to graze
+             call currentCohort%prt%GetBiomass(carbon12_element, &
+                  sapw_c, struct_c, leaf_c, fnrt_c, store_c, repro_c, alive_c, total_c)
+             if (total_c > rsnbl_math_prec) then
+                leaf_c_frac = leaf_c/total_c
+             else
+                leaf_c_frac = 0._r8
+             endif
 
              call FatesGrazing(currentCohort%prt, ft, currentPatch%land_use_label, currentCohort%height,currentCohort%npp_acc, &
-                  currentCohort%treelai)
+                  currentCohort%treelai, leaf_c_frac)
 
              ! Conduct Maintenance Turnover (parteh)
              if(debug) call currentCohort%prt%CheckMassConservation(ft,3)

main/FatesHistoryInterfaceMod.F90

@@ -2440,7 +2440,7 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites)
     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) :: herbivory     ! mass of leaves eaten by herbivores [kg/yr]
     real(r8) :: n_perm2            ! abundance per m2
     real(r8) :: patch_fracarea  ! Fraction of area for this patch
     real(r8) :: crown_area_covered ! accumulator variable for patch crown area
@@ -2978,8 +2978,9 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites)
                   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
-                  hio_grazing_si(io_si) = hio_grazing_si(io_si) + leaf_herbivory * n_perm2 / days_per_year / sec_per_day
+                  herbivory   = ccohort%prt%GetHerbivory(leaf_organ, carbon12_element) * days_per_year
+                  herbivory   = herbivory + ccohort%prt%GetHerbivory(store_organ, carbon12_element) * days_per_year
+                  hio_grazing_si(io_si) = hio_grazing_si(io_si) + herbivory * n_perm2 / days_per_year / sec_per_day
 
                   ! Woody State Variables (basal area growth increment)
                   if ( prt_params%woody(ft) == itrue) then

parameter_files/fates_params_noresm.json

@@ -1703,7 +1703,7 @@
       "dims": ["fates_landuseclass"],
       "long_name": "fraction of leaf biomass consumed by grazers per day",
       "units": "1/day",
-      "data": [0.007, 0.0, 0.02, 0.02, 0.04]
+      "data": [0.0035, 0.0, 0.01, 0.01, 0.02]
     },
     "fates_landuseclass_name":{
       "dtype": "string",