Add NVP (Non-Vascular Plant) PFT in CTSM–FATES

biogeochem/EDCanopyStructureMod.F90

@@ -19,6 +19,8 @@ module EDCanopyStructureMod
   use EDCohortDynamicsMod   , only : terminate_cohorts, terminate_cohort, fuse_cohorts
   use EDCohortDynamicsMod   , only : InitPRTObject
   use FatesAllometryMod     , only : tree_lai_sai
+  ! [PORTED by Hui Tang: NVP allometry routine for LAI/thickness update]
+  use FatesAllometryMod     , only : NVP_allom
   use EDTypesMod            , only : ed_site_type
   use EDTypesMod            , only : set_patchno
   use FatesAllometryMod     , only : VegAreaLayer
@@ -36,6 +38,9 @@ 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_nvp              ! [PORTED by Hui Tang: NVP master flag]
+  use FatesInterfaceTypesMod     , only : hlm_nvp_rad_model_ground ! [PORTED by Hui Tang: NVP radiation model switch]
+  use LeafBiophysicsMod          , only : lb_params                ! [PORTED by Hui Tang: NVP PFT identification]
   use FatesInterfaceTypesMod     , only : numpft
   use FatesInterfaceTypesMod, only : bc_in_type
   use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
@@ -1143,19 +1148,26 @@ subroutine leaf_area_profile( currentSite )
                       remainder = dinc_vai(iv)
                    end if
 
-                   cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + &
-                        remainder * fleaf * currentCohort%c_area/cpatch%total_canopy_area
-
-                   cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + &
-                        remainder * fleaf * currentCohort%c_area/cpatch%total_canopy_area * &
-                        fraction_exposed
-
-                   cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + &
-                        remainder * (1._r8 - fleaf) * currentCohort%c_area/cpatch%total_canopy_area
-
-                   cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + &
-                        remainder * (1._r8 - fleaf) * currentCohort%c_area/cpatch%total_canopy_area * &
-                        fraction_exposed
+                   ! [PORTED by Hui Tang: Approach A excludes NVP from Norman solver layer profiles]
+                   ! canopy_area_profile is always accumulated (needed for photosynthesis leaf_area weighting).
+                   ! LAI/SAI are excluded for NVP so the Norman solver uses soil albedo as its
+                   ! lower boundary without NVP absorbing as a canopy layer (avoids double-counting with R3).
+                   if (.not. (hlm_nvp_rad_model_ground == itrue .and. &
+                               lb_params%stomatal_intercept(ft) <= nearzero)) then
+                      cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + &
+                           remainder * fleaf * currentCohort%c_area/cpatch%total_canopy_area
+
+                      cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + &
+                           remainder * fleaf * currentCohort%c_area/cpatch%total_canopy_area * &
+                           fraction_exposed
+
+                      cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + &
+                           remainder * (1._r8 - fleaf) * currentCohort%c_area/cpatch%total_canopy_area
+
+                      cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + &
+                           remainder * (1._r8 - fleaf) * currentCohort%c_area/cpatch%total_canopy_area * &
+                           fraction_exposed
+                   end if
 
                    cpatch%canopy_area_profile(cl,ft,iv) = cpatch%canopy_area_profile(cl,ft,iv) + &
                         currentCohort%c_area/cpatch%total_canopy_area
@@ -1176,17 +1188,22 @@ subroutine leaf_area_profile( currentSite )
                         elai_layer,esai_layer,tlai_layer,tsai_layer)
 
 
-                   cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + &
-                        tlai_layer * currentCohort%c_area/cpatch%total_canopy_area
+                   ! [PORTED by Hui Tang: Approach A excludes NVP from Norman solver layer profiles]
+                   ! canopy_area_profile is always accumulated (needed for photosynthesis leaf_area weighting).
+                   if (.not. (hlm_nvp_rad_model_ground == itrue .and. &
+                               lb_params%stomatal_intercept(ft) <= nearzero)) then
+                      cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + &
+                           tlai_layer * currentCohort%c_area/cpatch%total_canopy_area
 
-                   cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + &
-                        elai_layer * currentCohort%c_area/cpatch%total_canopy_area
+                      cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + &
+                           elai_layer * currentCohort%c_area/cpatch%total_canopy_area
 
-                   cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + &
-                        tsai_layer * currentCohort%c_area/cpatch%total_canopy_area
+                      cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + &
+                           tsai_layer * currentCohort%c_area/cpatch%total_canopy_area
 
-                   cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + &
-                        esai_layer * currentCohort%c_area/cpatch%total_canopy_area
+                      cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + &
+                           esai_layer * currentCohort%c_area/cpatch%total_canopy_area
+                   end if
 
                    cpatch%canopy_area_profile(cl,ft,iv) = cpatch%canopy_area_profile(cl,ft,iv) + &
                         currentCohort%c_area/cpatch%total_canopy_area
@@ -1349,12 +1366,18 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
     ! Locals
     type (fates_cohort_type) , pointer :: currentCohort
     integer :: s, ifp, c, p
+    ! [PORTED by Hui Tang: loop indices for NVP LAI subtraction loop]                                                                                                                                       
+    integer :: cl, ft, iv  
+
     type (fates_patch_type)  , pointer :: currentPatch
     real(r8) :: bare_frac_area
     real(r8) :: total_patch_area
     real(r8) :: total_canopy_area
     real(r8) :: total_patch_leaf_stem_area
-    real(r8) :: weight  ! Weighting for cohort variables in patch
+    real(r8) :: weight           ! Weighting for cohort variables in patch
+    ! [PORTED by Hui Tang: NVP aggregation locals]
+    real(r8) :: nvp_crown_area   ! sum of c_area over all NVP cohorts in patch [m2/ha]
+                                 ! used to compute nvp_frac_pa and to normalize nvp_dz_pa
 
     do s = 1,nsites
 
@@ -1364,6 +1387,10 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        bc_out(s)%dleaf_pa(:) = 0._r8
        bc_out(s)%z0m_pa(:) = 0._r8
        bc_out(s)%displa_pa(:) = 0._r8
+       ! [PORTED by Hui Tang: zero NVP patch-level fields before aggregation]
+       bc_out(s)%nvp_dz_pa(:)   = 0._r8
+       bc_out(s)%nvp_frac_pa(:) = 0._r8
+       bc_out(s)%lai_nvp_pa(:)  = 0._r8
 
        currentPatch => sites(s)%oldest_patch
        c = fcolumn(s)
@@ -1470,6 +1497,54 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
              bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai')
              bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai')
 
+             ! [PORTED by Hui Tang: remove NVP thallus area contribution from patch LAI sent to CLM]
+             ! NVP is a ground-level mat, not part of the vascular canopy.
+             ! In Approach A (elai_profile=0 for NVP) this subtracts zero; in Approach B it removes
+             ! the NVP leaf area that was accumulated by calc_areaindex.
+             if (hlm_use_nvp == itrue) then
+                do cl = 1, currentPatch%NCL_p
+                   do ft = 1, numpft
+                      if (lb_params%stomatal_intercept(ft) <= nearzero) then
+                         do iv = 1, currentPatch%nrad(cl,ft)
+                            bc_out(s)%elai_pa(ifp) = bc_out(s)%elai_pa(ifp) - &
+                                 currentPatch%canopy_area_profile(cl,ft,iv) * currentPatch%elai_profile(cl,ft,iv)
+                            bc_out(s)%tlai_pa(ifp) = bc_out(s)%tlai_pa(ifp) - &
+                                 currentPatch%canopy_area_profile(cl,ft,iv) * currentPatch%tlai_profile(cl,ft,iv)
+                         end do
+                      end if
+                   end do
+                end do
+                bc_out(s)%elai_pa(ifp) = max(0._r8, bc_out(s)%elai_pa(ifp))
+                bc_out(s)%tlai_pa(ifp) = max(0._r8, bc_out(s)%tlai_pa(ifp))
+             end if
+
+             ! [PORTED by Hui Tang: cohort-to-patch NVP aggregation]
+             ! Accumulate nvp_dz_pa (area-weighted mean NVP thickness where present)
+             ! and nvp_frac_pa (fraction of patch covered by NVP).
+             if (hlm_use_nvp == itrue) then
+                nvp_crown_area = 0._r8
+                currentCohort => currentPatch%shortest
+                do while (associated(currentCohort))
+                   if (currentCohort%nvp_dz > nearzero) then
+                      bc_out(s)%nvp_dz_pa(ifp) = bc_out(s)%nvp_dz_pa(ifp) + &
+                           currentCohort%nvp_dz * currentCohort%c_area
+                      ! [PORTED by Hui Tang: accumulate area-weighted NVP LAI for Beer's law]
+                      bc_out(s)%lai_nvp_pa(ifp) = bc_out(s)%lai_nvp_pa(ifp) + &
+                           currentCohort%treelai * currentCohort%c_area
+                      nvp_crown_area = nvp_crown_area + currentCohort%c_area
+                   end if
+                   currentCohort => currentCohort%taller
+                end do
+                if (nvp_crown_area > nearzero) then
+                   ! Normalize dz and LAI to get area-weighted means within NVP-covered area
+                   bc_out(s)%nvp_dz_pa(ifp)   = bc_out(s)%nvp_dz_pa(ifp)  / nvp_crown_area
+                   ! [PORTED by Hui Tang: normalize lai to area-weighted mean within NVP-covered area]
+                   bc_out(s)%lai_nvp_pa(ifp)   = bc_out(s)%lai_nvp_pa(ifp) / nvp_crown_area
+                   ! Coverage fraction relative to full patch area
+                   bc_out(s)%nvp_frac_pa(ifp) = nvp_crown_area / AREA
+                end if
+             end if
+
              ! Fraction of vegetation free of snow. This is used to flag those
              ! patches which shall under-go photosynthesis
              ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let
@@ -1713,19 +1788,30 @@ 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 )
+   ! [PORTED by Hui Tang: NVP cohorts use NVP_allom; trees use tree_lai_sai]
+   if (hlm_use_nvp == itrue .and. &
+       lb_params%stomatal_intercept(currentCohort%pft) <= nearzero) then
+      ! NVP: derive treelai and nvp_dz from leaf carbon
+      call NVP_allom(currentCohort%treelai, leaf_c, currentCohort%nvp_dz, &
+                            inverse=.false.)
+      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
 
-   ! Do not update stem area index of SP vegetation
-   if (hlm_use_sp .eq. ifalse) then
-      currentCohort%treesai = treesai
+      ! Number of actual vegetation layers in this cohort's crown
+      currentCohort%nv = GetNVegLayers(currentCohort%treelai+currentCohort%treesai)
    end if
-
-   ! Number of actual vegetation layers in this cohort's crown
-   currentCohort%nv = GetNVegLayers(currentCohort%treelai+currentCohort%treesai)
-
+
   end subroutine UpdateCohortLAI
 
   ! ===============================================================================================

biogeochem/FatesAllometryMod.F90

@@ -99,6 +99,8 @@ module FatesAllometryMod
   use FatesGlobals     , only : endrun => fates_endrun
   use FatesGlobals     , only : FatesWarn,N2S,A2S,I2S
   use EDParamsMod      , only : nlevleaf,dinc_vai,dlower_vai
+  ! [PORTED by Hui Tang: NVP allometry parameters from FATES parameter file]
+  use EDParamsMod      , only : nvp_bulk_density, nvp_sla
   use DamageMainMod    , only : GetCrownReduction
   use LeafBiophysicsMod, only : DecayCoeffVcmax
 
@@ -124,6 +126,8 @@ module FatesAllometryMod
   public :: set_root_fraction  ! Generic wrapper to calculate normalized
                                ! root profiles
   public :: leafc_from_treelai ! Calculate target leaf carbon for a given treelai for SP mode
+  ! [PORTED by Hui Tang: NVP (moss/lichen) bidirectional LAI<->biomass and LAI->thickness conversion]
+  public :: NVP_allom
   public :: VegAreaLayer
 
   public :: tree_lai_sai       ! LAI and SAI calculations must work together, thus they
@@ -2800,6 +2804,8 @@ subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot)
     integer, parameter :: jackson_beta_profile_type   = 1
     integer, parameter :: exponential_1p_profile_type = 2
     integer, parameter :: exponential_2p_profile_type = 3
+    ! [PORTED by Hui Tang: mode 4 = NVP (no soil roots); root_fraction stays all-zero]
+    integer, parameter :: no_root_profile_type         = 4
 
     integer :: root_profile_type
     integer :: corr_id(1)        ! This is the bin with largest fraction
@@ -2835,9 +2841,12 @@ subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot)
        call exponential_1p_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), prt_params%fnrt_prof_a(ft)) 
     case ( jackson_beta_profile_type )
        call jackson_beta_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), prt_params%fnrt_prof_a(ft))
-    case ( exponential_2p_profile_type ) 
-       call exponential_2p_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), & 
+    case ( exponential_2p_profile_type )
+       call exponential_2p_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), &
              prt_params%fnrt_prof_a(ft),prt_params%fnrt_prof_b(ft))
+    ! [PORTED by Hui Tang: NVP has no soil roots; root_fraction already zeroed above]
+    case ( no_root_profile_type )
+       ! No root profile: root_fraction(:) = 0._r8 (initialized at line 2827, no further action needed)
 
     case default
        write(fates_log(),*) 'An undefined root profile type was specified'
@@ -2846,9 +2855,12 @@ subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot)
     end select
 
 
-    correction = 1._r8 - sum(root_fraction)
-    corr_id = maxloc(root_fraction)
-    root_fraction(corr_id(1)) = root_fraction(corr_id(1)) + correction
+    ! [PORTED by Hui Tang: skip normalization for no_root_profile_type; all zeros is intentional]
+    if (nint(prt_params%fnrt_prof_mode(ft)) /= no_root_profile_type) then
+       correction = 1._r8 - sum(root_fraction)
+       corr_id = maxloc(root_fraction)
+       root_fraction(corr_id(1)) = root_fraction(corr_id(1)) + correction
+    end if
 
 
 
@@ -3544,6 +3556,57 @@ subroutine size2dbh(size,ipft,dbh,dbh_maxh)
    end subroutine size2dbh
    ! ============================================================================
 
+   ! [PORTED by Hui Tang: NVP (moss/lichen) allometry - LAI<->biomass and LAI->layer thickness]
+   subroutine NVP_allom(lai, leaf_c, nvp_dz, inverse)
+      !
+      ! !DESCRIPTION:
+      ! Bidirectional conversion between NVP leaf carbon and LAI, and diagnosis
+      ! of the NVP layer thickness from LAI.
+      !
+      ! Parameters (from Porada et al. 2013 and related literature, read from FATES parameter file):
+      !   nvp_sla          [m2 g-1]  specific leaf area of NVP
+      !   nvp_bulk_density [kg m-3]  NVP bulk density
+      !   nv_c2b   = 2.32            carbon-to-biomass ratio of NVP (kgC / kgDM)
+      !   nv_t2l   = 1.0             total-to-leaf biomass ratio (all biomass is leaf for NVP)
+      !
+      ! When inverse = .false.: leaf_c  ->  lai, nvp_dz
+      ! When inverse = .true. : lai     ->  leaf_c  (nvp_dz still diagnosed from lai)
+      !
+      ! nvp_dz = lai / (nvp_sla [m2 g-1] * 1e3 [g kg-1] * nvp_bulk_density [kg m-3])
+      !        = lai / (nvp_sla * g_per_kg * nvp_bulk_density)
+      !
+      ! !ARGUMENTS:
+      real(r8), intent(inout) :: lai     ! NVP leaf area index [m2 leaf / m2 crown]
+      real(r8), intent(inout) :: leaf_c  ! NVP leaf carbon per unit crown area [kgC m-2]
+      real(r8), intent(out)   :: nvp_dz  ! NVP layer thickness [m]
+      logical,  intent(in)    :: inverse ! .false. = leaf_c->lai; .true. = lai->leaf_c
+      !
+      ! !LOCAL VARIABLES:
+      real(r8), parameter :: nv_c2b  = 2.32_r8 ! carbon to dry-mass ratio [kgC kgDM-1]
+      real(r8), parameter :: nv_t2l  = 1.0_r8  ! total-to-leaf biomass ratio [-]
+      !------------------------------------------------------------------------
+
+      if (inverse) then
+         ! lai -> leaf_c
+         ! leaf_mass_per_area [kgDM m-2] = lai [m2 m-2] / (nvp_sla [m2 g-1] * g_per_kg [g kg-1])
+         ! leaf_c [kgC m-2] = leaf_mass_per_area * nv_c2b / nv_t2l
+         leaf_c = (lai / (nvp_sla * g_per_kg)) * nv_c2b / nv_t2l
+      else
+         ! leaf_c -> lai
+         ! leaf_mass [kgDM m-2] = leaf_c * nv_t2l / nv_c2b
+         ! lai = leaf_mass * nvp_sla * g_per_kg
+         lai = (leaf_c * nv_t2l / nv_c2b) * nvp_sla * g_per_kg
+      end if
 
+      ! Diagnose layer thickness from LAI regardless of direction
+      ! nvp_dz [m] = lai [m2 m-2] / (nvp_sla [m2 g-1] * g_per_kg [g kg-1] * nvp_bulk_density [kg m-3])
+      if (nvp_bulk_density > nearzero) then
+         nvp_dz = lai / (nvp_sla * g_per_kg * nvp_bulk_density)
+      else
+         nvp_dz = 0.0_r8
+      end if
+
+   end subroutine NVP_allom
+   ! ============================================================================
 
 end module FatesAllometryMod