Toy moss model

biogeochem/EDCanopyStructureMod.F90

@@ -12,6 +12,7 @@ module EDCanopyStructureMod
   use FatesConstantsMod     , only : rsnbl_math_prec
   use FatesConstantsMod     , only : nocomp_bareground
   use FatesConstantsMod,      only : i_term_mort_type_canlev
+  use FatesConstantsMod,      only : m2_per_ha
   use FatesGlobals          , only : fates_log
   use EDPftvarcon           , only : EDPftvarcon_inst
   use PRTParametersMod      , only : prt_params
@@ -2232,7 +2233,9 @@ end subroutine UpdatePatchLAI
   subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area)
 
    ! Update LAI and related variables for a given cohort
-
+
+   use FatesMossMod, only : SLA_M2LEAF_PER_KGMOSS
+
    ! Arguments
    type(fates_cohort_type),intent(inout), target   :: currentCohort
    real(r8), intent(in) :: canopy_layer_tlai(nclmax)  ! total leaf area index of each canopy layer
@@ -2244,20 +2247,26 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area)
 
    ! Obtain the leaf carbon
    leaf_c = currentCohort%prt%GetState(leaf_organ,carbon12_element)
-
-   ! Note that tree_lai has an internal check on the canopy location
-   call  tree_lai_sai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n,           &
-          currentCohort%canopy_layer, canopy_layer_tlai, currentCohort%vcmax25top, currentCohort%dbh, currentCohort%crowndamage,          &
-          currentCohort%canopy_trim, currentCohort%efstem_coh, 4, currentCohort%treelai, treesai )
-
-   ! Do not update stem area index of SP vegetation
-   if (hlm_use_sp .eq. ifalse) then
-      currentCohort%treesai = treesai
-   end if
 
-   ! Number of actual vegetation layers in this cohort's crown
-   currentCohort%nv = GetNVegLayers(currentCohort%treelai+currentCohort%treesai)
+   if (prt_params%moss(currentCohort%pft) == itrue) then
+      currentCohort%treelai = leaf_c * currentCohort%n / m2_per_ha * SLA_M2LEAF_PER_KGMOSS
+      currentCohort%treesai = 0._r8
+      currentCohort%nv = 1
+   else
+
+      ! Note that tree_lai has an internal check on the canopy location
+      call  tree_lai_sai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n,           &
+             currentCohort%canopy_layer, canopy_layer_tlai, currentCohort%vcmax25top, currentCohort%dbh,    currentCohort%crowndamage,          &
+             currentCohort%canopy_trim, currentCohort%efstem_coh, 4, currentCohort%treelai, treesai )
+
+      ! Do not update stem area index of SP vegetation
+      if (hlm_use_sp .eq. ifalse) then
+         currentCohort%treesai = treesai
+      end if
 
+      ! Number of actual vegetation layers in this cohort's crown
+      currentCohort%nv = GetNVegLayers(currentCohort%treelai+currentCohort%treesai)
+   end if
   end subroutine UpdateCohortLAI
 
   ! ===============================================================================================

biogeochem/EDPatchDynamicsMod.F90

@@ -3191,7 +3191,28 @@ subroutine fuse_2_patches(csite, dp, rp)
     rp%frac_burnt           = (dp%frac_burnt*dp%area + rp%frac_burnt*rp%area) * inv_sum_area
     rp%btran_ft(:)          = (dp%btran_ft(:)*dp%area + rp%btran_ft(:)*rp%area) * inv_sum_area
     rp%zstar                = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area
-    rp%c_stomata            = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area
+    if (dp%IsAllMoss() .or. rp%IsAllMoss()) then
+       ! We set c_stomata to the unset real value in FatesPlantRespPhotosynthDrive() because mosses
+       ! are expected to have zero resistance. If exactly one fusing patch is all moss, take
+       ! c_stomata from that. This is not accurate, but it's only a diagnostic and for
+       ! one timestep at most (until the next call of FatesPlantRespPhotosynthDrive).
+       ! TODO:
+       !    * If keeping fates_unset_r8 kludge in FatesPlantRespPhotosynthDrive(), code up a way
+       !      to quickly recalculate c_stomata without having to call FatesPlantRespPhotosynthDrive()
+       !      again.
+       !    * Otherwise, c_stomata should always have meaningful values, so remove this special
+       !      handling and just calculate the new value as the area-weighted mean of the donor and
+       !      recipient patches.
+       if (.not. dp%IsAllMoss()) then
+          rp%c_stomata      = dp%c_stomata
+       else if (.not. rp%IsAllMoss()) then
+          rp%c_stomata      = rp%c_stomata
+       else
+          rp%c_stomata      = fates_unset_r8
+       end if
+    else
+       rp%c_stomata         = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area
+    end if
     rp%c_lblayer            = (dp%c_lblayer*dp%area + rp%c_lblayer*rp%area) * inv_sum_area
 
     ! Radiation

biogeochem/EDPhysiologyMod.F90

@@ -702,7 +702,12 @@ subroutine trim_canopy( currentSite )
                currentCohort%dbh, currentCohort%crowndamage, currentCohort%canopy_trim, &
                currentCohort%efstem_coh, 0, currentCohort%treelai, currentCohort%treesai )
 
-          currentCohort%nv = GetNVegLayers(currentCohort%treelai+currentCohort%treesai)
+          ! TODO: Would be safer to add this handling to GetNVegLayers()
+          if (prt_params%moss(ipft) == itrue) then
+             currentCohort%nv = 1
+          else
+             currentCohort%nv = GetNVegLayers(currentCohort%treelai+currentCohort%treesai)
+          end if
 
           leaf_veg_frac = currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
 
@@ -1885,6 +1890,11 @@ subroutine satellite_phenology(currentSite, bc_in)
              call endrun(msg=errMsg(sourcefile, __LINE__))
           end if
 
+          if (prt_params%moss(fates_pft) == itrue) then
+             write(fates_log(),*) 'Moss in SP mode'
+             call endrun(msg=errMsg(sourcefile, __LINE__))
+          end if
+
           ! Call routine to invert SP drivers into cohort properites.
           call assign_cohort_SP_properties(currentCohort, currentSite%sp_htop(fates_pft), currentSite%sp_tlai(fates_pft)     , currentSite%sp_tsai(fates_pft),currentPatch%area,ifalse,leaf_c)
 

biogeochem/FatesPatchMod.F90

@@ -241,6 +241,7 @@ module FatesPatchMod
       procedure :: FreeMemory
       procedure :: Dump
       procedure :: CheckVars
+      procedure :: IsAllMoss
 
   end type fates_patch_type
 
@@ -1286,6 +1287,33 @@ subroutine CheckVars(this, var_aliases, return_code)
 
     end subroutine CheckVars
 
+    !===========================================================================
+
+    function IsAllMoss(this)
+      !
+      ! DESCRIPTION:
+      !  Checks whether all cohorts on patch are moss
+
+      ! ARGUMENTS:
+      class(fates_patch_type), intent(inout) :: this
+
+      ! LOCALS:
+      type(fates_cohort_type), pointer :: currentCohort ! cohort object
+
+      ! RESULT:
+      logical :: IsAllMoss
+
+      currentCohort => this%tallest
+      IsAllMoss = associated(currentCohort)  ! If patch has no cohorts, it's not all moss
+      do while(associated(currentCohort) .and. IsAllMoss)
+        if (prt_params%moss(currentCohort%pft) == ifalse) then
+          IsAllMoss = .false.
+        end if
+        currentCohort => currentCohort%shorter
+      end do
+
+    end function IsAllMoss
+
     !===========================================================================  
 
 end module FatesPatchMod

biogeophys/CMakeLists.txt

@@ -1,5 +1,6 @@
 list(APPEND fates_sources
   FatesHydroWTFMod.F90
+  FatesMossMod.F90
   LeafBiophysicsMod.F90
   FatesPlantHydraulicsMod.F90)
 

biogeophys/FatesMossMod.F90

@@ -0,0 +1,232 @@
+module FatesMossMod
+
+! ==============================================================================================
+! This module contains the relevant code for moss.
+!
+! WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
+!
+!  MOSS IS AN EXPERIMENTAL OPTION THAT IS STILL UNDERGOING TESTING.
+!
+! WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
+!
+! ==============================================================================================
+
+  use FatesConstantsMod , only: r8 => fates_r8, nearzero
+  use FatesConstantsMod , only: umolC_to_kgC
+  use FatesConstantsMod , only: days_per_sec, days_per_year
+  use FatesGlobals, only      : endrun => fates_endrun
+  use FatesGlobals, only      : fates_log
+
+
+  implicit none
+
+  ! Constants from Bonan & Korzukhin (1989)
+  real(r8), parameter :: Q_KG_PER_KGMOSS = 0.12       ! Respiration? parameter
+  real(r8), parameter :: B_KG_PER_KGMOSS = 0.136      ! Mortality? parameter
+  real(r8), parameter :: EXT = 0.5      ! Light extinction coefficient
+  real(r8), parameter :: SLA_M2LEAF_PER_KGMOSS = 1.0       ! Specific leaf area (m2/kg)
+  real(r8), parameter :: FRAC_ASSIM_TO_SPORES = 0.001 ! Proportion of assimilation spent on reproduction. UNUSED but keeping for possible future work.
+  real(r8), parameter :: PMAX_KG_PER_M2_PER_YR = 0.35    ! Maximum moss production (kg/m2/yr)
+  real(r8), parameter :: LRMIN = 0.01   ! Light compensation point
+  real(r8), parameter :: LRMAX = 0.05   ! Light compensation point
+
+  ! Constants from UVAFME (maybe also from Bonan & Korzukhin, 1989, but unsure)
+  real(r8), parameter :: BULK_MOSS_KG_PER_M3 = 18       ! Moss bulk density (kg/m3)
+  real(r8), parameter :: FCGF_ALFF_THRESH = 0.75  ! Above this level of light on forest floor, moss assimilation starts to decrease
+  real(r8), parameter :: FCGF_INTERCEPT = 1.5625  ! Intercept parameter in fcgf equation
+  real(r8), parameter :: DLGF_DECLIT_THRESH = 0  ! Above this level of deciduous litter (kg/m2 plot), moss assimilation starts to decrease
+  real(r8), parameter :: DLGF_DECAY = 2.932  ! Decay parameter of moss assimilation with increasing deciduous leaf litter
+
+
+  ! PUBLIC MEMBER FUNCTIONS:
+  public :: light_growth_multiplier
+  public :: forest_cover_growth_multiplier
+  public :: litter_growth_multiplier
+  public :: moss
+
+contains
+
+!------------------------------------------------------------------------------
+function available_light_under_canopy_and_moss(canopy_lai, moss_biom_kg_per_m2plot_before) result(al)
+  real(r8), intent(in)    :: canopy_lai  ! Leaf area index of canopy (i.e., excluding moss) (m2 leaves / m2 plot)
+  real(r8), intent(in)    :: moss_biom_kg_per_m2plot_before  ! Moss biomass (kg/m2) before this timestep
+  real(r8) :: moss_lai  ! Leaf area index of moss (m2 leaves / m2 plot)
+  real(r8) :: lai       ! Total leaf area index, canopy + moss (m2 leaves / m2 plot)
+  real(r8) :: al        ! Available light
+
+  moss_lai = moss_biom_kg_per_m2plot_before * SLA_M2LEAF_PER_KGMOSS
+  lai = canopy_lai + moss_lai
+  al = exp(-1.0*EXT*lai)
+end function available_light_under_canopy_and_moss
+
+!------------------------------------------------------------------------------
+function light_growth_multiplier(canopy_lai, moss_biom_kg_per_m2plot_before) result(algf)
+  real(r8), intent(in)    :: canopy_lai  ! Leaf area index of canopy (i.e., excluding moss) (m2 leaves / m2 plot)
+  real(r8), intent(in)    :: moss_biom_kg_per_m2plot_before  ! Moss biomass (kg/m2) before this timestep
+  real(r8) :: al        ! Available light
+  real(r8) :: algf      ! Light growth multiplier
+
+  al = available_light_under_canopy_and_moss(canopy_lai, moss_biom_kg_per_m2plot_before)
+  algf = (al - LRMIN)/(LRMAX - LRMIN)
+  algf = max(0.0, algf)
+  algf = min(1.0, algf)
+end function light_growth_multiplier
+
+!------------------------------------------------------------------------------
+function forest_cover_growth_multiplier(alff) result(fcgf)
+  real(r8), intent(in)    :: alff    ! Available light on the forest floor (0-1)
+  real(r8) :: fcgf      ! Forest cover growth factor
+  if (alff > FCGF_ALFF_THRESH) then
+      fcgf = FCGF_INTERCEPT - alff**2
+  else
+      fcgf = 1.0
+  end if
+  if (fcgf > 1.0) fcgf = 1.0
+  if (fcgf < 0.0) fcgf = 0.0
+end function forest_cover_growth_multiplier
+
+!------------------------------------------------------------------------------
+function litter_growth_multiplier(decid_litter) result(dlgf)
+  real(r8), intent(in)    :: decid_litter    ! Fresh deciduous leaf litter (kg/m2)
+  real(r8) :: dlgf      ! Deciduous leaf litter growth multiplier
+  if (decid_litter > DLGF_DECLIT_THRESH) then
+      dlgf = exp(-DLGF_DECAY * decid_litter)
+      if (dlgf <= 0.0) dlgf = 0.0
+      if (dlgf >= 1.0) dlgf = 1.0
+  else
+      dlgf = 1.0
+  end if
+end function litter_growth_multiplier
+
+subroutine moss_biomass_change_kg_per_m2(q_kg_per_kg_moss_in, b_kg_per_kg_moss_in, assim_eff, moss_biom_before, moss_resp, moss_mort, moss_biom_after)
+  ! ALL UNITS KG / M2 PLOT UNLESS OTHERWISE SPECIFIED
+  real(r8), intent(in) :: q_kg_per_kg_moss_in  ! Respiration? parameter
+  real(r8), intent(in) :: b_kg_per_kg_moss_in  ! Mortality? parameter
+  real(r8), intent(in) :: assim_eff ! Assimilation (kg/m2)
+  real(r8), intent(in) :: moss_biom_before   ! Moss biomass (kg/m2) before this timestep
+  real(r8), intent(out) :: moss_resp  ! Moss respiration (kg/m2) during this timestep
+  real(r8), intent(out) :: moss_mort  ! Moss mortality (kg/m2) during this timestep
+  real(r8), intent(out) :: moss_biom_after   ! Moss biomass (kg/m2) after this timestep
+
+  ! Local variables
+  real(r8) :: moss_respmort  ! Sum of moss respiration and mortality (kg/m2) during this timestep
+  real(r8) :: qb_sum  ! Sum of Q and B
+  real(r8) :: moss_biom_change     ! Change in moss biomass (kg/m2) during this timestep
+
+  ! Total losses to respiration and mortality
+  moss_resp = q_kg_per_kg_moss_in * moss_biom_before
+  moss_mort = b_kg_per_kg_moss_in * moss_biom_before
+  moss_respmort = moss_resp + moss_mort
+
+  ! Not enough moss to account for mortality/respiration
+  ! Set moss loss to all of current biomass
+  if (moss_respmort > moss_biom_before) then
+      qb_sum = q_kg_per_kg_moss_in + b_kg_per_kg_moss_in
+      moss_resp = moss_biom_before * (q_kg_per_kg_moss_in / qb_sum)
+      moss_mort = moss_biom_before * (b_kg_per_kg_moss_in / qb_sum)
+      moss_respmort = moss_resp + moss_mort
+  end if
+
+  ! Net change in moss biomass
+  moss_biom_change = assim_eff - moss_respmort
+
+  ! Update biomass
+  moss_biom_after = moss_biom_before + moss_biom_change
+end subroutine moss_biomass_change_kg_per_m2
+
+!------------------------------------------------------------------------------
+subroutine moss(alff, canopy_lai, decid_litter, dtime, moss_biom_kg_per_m2plot_inout, moss_to_litter_flux_kg_per_m2plot, moss_to_atmos_flux_kg_per_m2plot, &
+                livemoss_depth_m, psn_z, anet_av_z, lmr_z)
+  !
+  !  Calculates annual moss growth and mortality
+  !  Adapted from Bonan and Korzukhin 1989 Vegetatio 84:31-44
+  !  Further adapted from Foster et al. (2019, Ecol. Mod., doi: 10.1016/j.ecolmodel.2019.108765)
+  !
+
+  ! Arguments
+  real(r8), intent(in)    :: alff    ! Available light on the forest floor (0-1)
+  real(r8), intent(in)    :: canopy_lai  ! Leaf area index of canopy (i.e., excluding moss) (m2 leaves / m2 plot)
+  real(r8), intent(in)    :: decid_litter  ! Fresh deciduous leaf litter (kg/m2)
+  real(r8), intent(in)    :: dtime  ! Time step length (s/timestep)
+  real(r8), intent(inout) :: moss_biom_kg_per_m2plot_inout  ! Moss biomass (kg/m2)
+  real(r8), intent(out)   :: moss_to_litter_flux_kg_per_m2plot  ! Flux from moss to litter (kg/m2)
+  real(r8), intent(out)   :: moss_to_atmos_flux_kg_per_m2plot  ! Flux from moss to atmosphere (kg/m2)
+  real(r8), intent(out)   :: livemoss_depth_m  ! Depth (m) of live moss layer
+  ! Outputs for FatesPlantRespPhotosynthDrive()
+  real(r8), intent(out)   :: psn_z  ! GPP [umolC/m2leaf/s]
+  real(r8), intent(out)   :: anet_av_z  ! "net leaf photosynthesis" [umol C/m2leaf/timestep]
+  real(r8), intent(out)   :: lmr_z  ! leaf maintenance (dark) respiration [umolC/m2leaf/s]
+
+  ! Local variables
+  real(r8) :: per_year_to_per_sec
+  real(r8) :: per_year_to_per_timestep
+  real(r8) :: moss_biom_kg_per_plot_before    ! Moss biomass (kg per plot) before this timestep
+  real(r8) :: moss_biom_kg_per_m2plot_before  ! Moss biomass (kg/m2) before this timestep
+  real(r8) :: moss_biom_kg_per_m2plot_after   ! Moss biomass (kg/m2) after this timestep
+  real(r8) :: algf      ! Available light growth factor
+  real(r8) :: fcgf      ! Forest cover growth factor
+  real(r8) :: dlgf      ! Deciduous leaf litter growth factor
+  real(r8) :: ddgf      ! Moisture growth factor
+  real(r8) :: assim_kg_per_m2leaf     ! Moss assimilation rate (kg/m2)
+  real(r8) :: assim_eff_kg_per_kgmoss ! Effective assimilation (kg/kg)
+  real(r8) :: assim_eff_kg_per_m2plot ! Assimilation (kg/m2)
+  real(r8) :: moss_biom_change_kg_per_m2plot     ! Change in moss biomass (kg/m2) during this timestep
+  real(r8) :: moss_resp  ! Moss respiration (kg/m2) during this timestep
+  real(r8) :: moss_mort  ! Moss mortality (kg/m2) during this timestep
+
+  real(r8) :: moss_to_litter_flux_kg_per_plot    ! Flux from moss to litter (kg)
+  real(r8) :: moss_respmort_kg_per_kgmoss  ! Moss respiration and mortality (kg / kg moss) during this timestep
+  real(r8) :: moss_respmort_kg_per_m2plot  ! Moss respiration and mortality (kg/m2) during this timestep
+
+  ! Save this for later
+  moss_biom_kg_per_m2plot_before = moss_biom_kg_per_m2plot_inout
+
+  ! Light growth multiplier
+  algf = light_growth_multiplier(canopy_lai, moss_biom_kg_per_m2plot_before)
+
+  ! Forest cover growth multiplier
+  fcgf = forest_cover_growth_multiplier(alff)
+
+  ! Deciduous leaf litter growth multiplier
+  dlgf = litter_growth_multiplier(decid_litter)
+
+  ! Moisture growth factor
+  ! TODO: Implement this
+  ddgf = 1.0
+
+  ! Moss assimilation
+  assim_kg_per_m2leaf = PMAX_KG_PER_M2_PER_YR*algf*fcgf*dlgf*ddgf
+  assim_eff_kg_per_kgmoss = SLA_M2LEAF_PER_KGMOSS * assim_kg_per_m2leaf
+  assim_eff_kg_per_m2plot = (assim_eff_kg_per_kgmoss * moss_biom_kg_per_m2plot_before)
+
+  ! Get fluxes from moss and change in moss biomass
+  call moss_biomass_change_kg_per_m2(Q_KG_PER_KGMOSS, B_KG_PER_KGMOSS, assim_eff_kg_per_m2plot, moss_biom_kg_per_m2plot_before, moss_resp, moss_mort, moss_biom_kg_per_m2plot_after)
+
+  ! Get flux from live moss to litter and atmosphere
+  moss_to_litter_flux_kg_per_m2plot = moss_mort
+  moss_to_atmos_flux_kg_per_m2plot = moss_resp
+
+  ! Thickness of live moss layer (m)
+  livemoss_depth_m = moss_biom_kg_per_m2plot_after / BULK_MOSS_KG_PER_M3
+
+  ! Convert certain variables to their UVAFME outputs
+  ! TODO: Change these to what FATES needs
+  moss_biom_kg_per_m2plot_inout = moss_biom_kg_per_m2plot_after
+
+  ! Outputs for FatesPlantRespPhotosynthDrive()
+  per_year_to_per_sec = days_per_sec / days_per_year
+  per_year_to_per_timestep = per_year_to_per_sec * dtime
+  anet_av_z = assim_kg_per_m2leaf / umolC_to_kgC * per_year_to_per_timestep ! "net leaf photosynthesis" [umol C/m2leaf/timestep]
+  psn_z = anet_av_z * dtime ! GPP [umolC/m2leaf/s]. Moss model doesn't calculate, so assume GPP = NPP.
+  if (moss_biom_kg_per_m2plot_before < nearzero) then
+     lmr_z = 0._r8
+  else
+     lmr_z = moss_resp / umolC_to_kgC / moss_biom_kg_per_m2plot_before / SLA_M2LEAF_PER_KGMOSS * per_year_to_per_sec ! leaf maintenance (dark) respiration [umolC/m2leaf/s]
+  end if
+
+end subroutine moss
+
+! =====================================================================================
+
+
+end module FatesMossMod