Merge branch 'evasinha/lnd/add_sugarcane_based_on_ECOSMOS' (PR #7681)

This feature adds new parameterizations to represent the sugarcane crop in ELM.

The implementation is based on sugarcane implementation in the Agricultural version of the Integrated Biosphere Simulator (Agro-IBIS)
model (Cuadra et al., 2012 and Colmanetti et al., 2024).
It is important to note that sugarcane is modeled as a perennial crop,
and this requires modifying the crop and percrop parameters in the parameter file, along with accompanying changes in ELM.

The model was calibrated using the methodology used in Sinha et al. 2022.

[non BFB] for Crops since the parameter file is updated.

Author: peterdschwartz

components/elm/bld/namelist_files/namelist_defaults.xml

@@ -123,7 +123,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
      for the CLM2 data in the CESM distribution -->
 <!-- Plant function types (relative to {csmdata}) -->
 <paramfile use_cn=".false." use_fates=".false." >lnd/clm2/paramdata/clm_params_c180301.nc</paramfile>
-<paramfile use_crop=".true." >lnd/clm2/paramdata/clm_params_c220224.nc</paramfile>
+<paramfile use_crop=".true." >lnd/clm2/paramdata/clm_params_c251006.nc</paramfile>
 <paramfile nu_com="ECA" use_fates=".true." use_crop=".false." >lnd/clm2/paramdata/clm_params.cbgc.c07292018.nc</paramfile>
 <paramfile nu_com="ECA" use_cn=".true." use_crop=".false." >lnd/clm2/paramdata/clm_params.cbgc.c07292018.nc</paramfile>
 <paramfile use_crop=".false." use_cn=".true."  >lnd/clm2/paramdata/clm_params_c211124.nc</paramfile>

components/elm/src/biogeochem/AllocationMod.F90

@@ -108,6 +108,7 @@ module AllocationMod
   real(r8)              :: dayscrecover         !number of days to recover negative cpool
   real(r8), allocatable :: arepr(:)             !reproduction allocation coefficient
   real(r8), allocatable :: aroot(:)             !root allocation coefficient
+  real(r8), allocatable :: aerial(:)            !aboveground allocation coefficient
 
   !$acc declare create(bdnr                )
   !$acc declare create(dayscrecover        )
@@ -247,6 +248,7 @@ subroutine AllocationInit ( bounds, elm_fates)
     if ( crop_prog )then
        allocate(arepr(bounds%begp:bounds%endp)); arepr(bounds%begp : bounds%endp) = nan
        allocate(aroot(bounds%begp:bounds%endp)); aroot(bounds%begp : bounds%endp) = nan
+       allocate(aerial(bounds%begp:bounds%endp)); aerial(bounds%begp : bounds%endp) = nan
     end if
 
     ! Allocate scratch space for ECA and FATES/ECA
@@ -476,7 +478,7 @@ subroutine TotalNPDemand(num_soilp,filter_soilp, photosyns_vars, &
             canopystate_vars, crop_vars, cnstate_vars, dt)
     !! Assess the total plant N demand and P demand
     use pftvarcon   , only: npcropmin, declfact, bfact, aleaff, arootf, astemf, noveg
-     use pftvarcon  , only: iscft, percrop, nwcereal, nwcerealirrig
+     use pftvarcon  , only: iscft, percrop, nwcereal, nwcerealirrig, nsugarcane, nsugarcaneirrig
     use pftvarcon   , only: arooti, fleafi, allconsl, allconss, grperc, grpnow, nsoybean
     use elm_varpar  , only: nlevdecomp
     use elm_varcon  , only: nitrif_n2o_loss_frac, secspday
@@ -504,7 +506,21 @@ subroutine TotalNPDemand(num_soilp,filter_soilp, photosyns_vars, &
     integer, parameter :: cphase_gf = 3                                  !Crop phenology phase grain fill
     integer, parameter :: max_lai   = 1                                  !Maximum allowed lai
     integer :: ivt, fp, p
-    !-----------------------------------------------------------------------
+
+    ! LOCAL VARAIBLES for sugarcane allocation based on Colmanetti et al., 2024  10.1016/j.eja.2023.127061
+    real(r8), parameter :: rootd      = 0.3_r8   ! root decline parameter
+    real(r8), parameter :: sf1        = 0.02_r8  ! parameters for linear carbon allocation to the stalk
+    real(r8), parameter :: ipf1       = 1.0_r8   ! parameters for linear carbon allocation to the stalk
+    real(r8), parameter :: ecf2       = 0.065_r8 ! parameters for logarithmic carbon allocation to the stalk
+    real(r8), parameter :: ipf2       = 10.0_r8  ! parameters for logarithmic carbon allocation to the stalk
+    real(r8), parameter :: sf3        = 0.006_r8 ! parameters for linear carbon allocation to the stalk-sucrose
+    real(r8), parameter :: ipf3       = 5.0_r8   ! parameters for linear carbon allocation to the stalk-sucrose
+    real(r8), parameter :: ecf4       = 0.017_r8 ! parameters for logarithmic carbon allocation to the stalk-sucrose
+    real(r8), parameter :: ipf4       = 10.0_r8  ! parameters for logarithmic carbon allocation to the stalk-sucrose
+    real(r8) :: rmat, af1, af2, af3, af4, sipf3, sipf4
+
+  !-----------------------------------------------------------------------
+
     associate(                                                                                   &
          woody                        => veg_vp%woody                                      , & ! Input:  [real(r8) (:)   ]  woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub)
          froot_leaf                   => veg_vp%froot_leaf                                 , & ! Input:  [real(r8) (:)   ]  allocation parameter: new fine root C per new leaf C (gC/gC)
@@ -548,6 +564,8 @@ subroutine TotalNPDemand(num_soilp,filter_soilp, photosyns_vars, &
          huileaf                      => cnstate_vars%huileaf_patch                            , & ! Input:  [real(r8) (:)   ]  heat unit index needed from planting to leaf emergence
          huigrain                     => cnstate_vars%huigrain_patch                           , & ! Input:  [real(r8) (:)   ]  same to reach vegetative maturity
          croplive                     => crop_vars%croplive_patch                           , & ! Input:  [logical  (:)   ]  flag, true if planted, not harvested
+         nyrs_crop_active             => crop_vars%nyrs_crop_active_patch                      , & ! Input:  [integer (:)    ]  number of years this crop patch has been active
+         onset_gdd                    => cnstate_vars%onset_gdd_patch                          , & ! Output: [real(r8) (:)   ]  onset growing degree days
          peaklai                      => cnstate_vars%peaklai_patch                            , & ! Input:  [integer  (:)   ]  1: max allowed lai; 0: not at max
          !lgsf                        => cnstate_vars%lgsf_patch                               , & ! Input:  [real(r8) (:)   ]  long growing season factor [0-1]
          aleafi                       => cnstate_vars%aleafi_patch                             , & ! Output: [real(r8) (:)   ]  saved allocation coefficient from phase 2
@@ -879,15 +897,60 @@ subroutine TotalNPDemand(num_soilp,filter_soilp, photosyns_vars, &
                g1 = 0.25_r8
 
             else if (croplive(p) .and. percrop(ivt) == 1.0_r8) then
-               arepr(p) = 0._r8
-               aroot(p) = max(0._r8, min(1._r8, arooti(ivt) -   &
-                    (arooti(ivt) - arootf(ivt)) *  &
-                    min(1._r8, hui(p)/gddmaturity(p))))
-               fleaf = fleafi(ivt) * (exp(-bfact(ivt)) -         &
-                    exp(-bfact(ivt)*hui(p)/gddmaturity(p))) / &      ! replacing huigrain with gddmaturity since huigrain does not exist for perennial crops
-                    (exp(-bfact(ivt))-1) ! fraction alloc to leaf (from J Norman alloc curve)
-               aleaf(p) = max(1.e-5_r8, (1._r8 - aroot(p)) * fleaf)
-               astem(p) = 1._r8 - arepr(p) - aleaf(p) - aroot(p)
+               if (ivt==nsugarcane .or. ivt==nsugarcaneirrig) then
+
+                  rmat = 100.0_r8 * (onset_gdd(p)/gddmaturity(p))
+
+                  ! Aerial(p) ranges from 0 to 1 and is the allocation
+                  ! coefficient for aboveground carbon components, including
+                  ! stems and leaves
+
+                  ! Aboveground sugarcane plant growth decreases by 0.6 in the
+                  ! year when sugarcane is planted and by 1 in the following years
+                  if (nyrs_crop_active(p) == 0) then ! Year 1
+                     aerial(p) = (1.0_r8 - arootf(ivt)) * min(1.0_r8, (1 - exp(-(rootd * 0.6 * rmat))))
+                  else
+                     aerial(p) = (1.0_r8 - arootf(ivt)) * min(1.0_r8, (1 - exp(-(rootd * rmat))))
+                  end if
+
+                  aroot(p) = 1.0_r8 - aerial(p)
+
+                  ! Carbon allocation to stems is nearly linear at the start of
+                  ! the sugarcane growth cycle and then follows a logarithmic
+                  ! pattern until the cycle's end
+                  ! af1 - fraction of the aboveground carbon initially allocated to the stem
+                  ! af2 - fraction of the aboveground carbon allocated to the
+                  ! stem over most of the plant's lifespan
+                  af1 = max(0._r8, (rmat * sf1) - (ipf1 * sf1))
+                  af2 = max(0._r8, (1.0_r8 - exp(-((rmat * ecf2) - (ipf2 * ecf2)))) )
+
+                  ! Allocation coeff for aboveground carbon is split between stem and leaves
+                  astem(p) = min( (1.0_r8 - aleaff(ivt) - arootf(ivt)), (aerial(p) * max(af1, af2)) )
+                  aleaf(p) = max(1.e-5_r8, (aerial(p) - astem(p)))
+
+
+                  sipf3 = ipf1 + (100.0_r8 - ipf1) * (ipf3 / 100.0_r8)
+                  sipf4 = ipf2 + (100.0_r8 - ipf2) * (ipf4 / 100.0_r8)
+
+                  af3 = max(0._r8, (rmat * sf3) - (sipf3 * sf3))
+                  af4 = max(0._r8, (1.0_r8 - exp(-((rmat * ecf4) - (sipf4 * ecf4)))) )
+
+                  ! Carbon allocated to the stem is then divided between
+                  ! stem sucrose or reporductive part and stem structural components
+                  arepr(p) = astem(p) * max(0._r8, af3, af4)
+                  astem(p) = astem(p) - arepr(p)
+
+               else ! if perennial but not sugarcane
+                  arepr(p) = 0._r8
+                  aroot(p) = max(0._r8, min(1._r8, arooti(ivt) -   &
+                       (arooti(ivt) - arootf(ivt)) *  &
+                       min(1._r8, hui(p)/gddmaturity(p))))
+                  fleaf = fleafi(ivt) * (exp(-bfact(ivt)) -         &
+                       exp(-bfact(ivt)*hui(p)/gddmaturity(p))) / &      ! replacing huigrain with gddmaturity since huigrain does not exist for perennial crops
+                       (exp(-bfact(ivt))-1) ! fraction alloc to leaf (from J Norman alloc curve)
+                  aleaf(p) = max(1.e-5_r8, (1._r8 - aroot(p)) * fleaf)
+                  astem(p) = 1._r8 - arepr(p) - aleaf(p) - aroot(p)
+               end if
 
                f1 = aroot(p) / aleaf(p)
                f3 = astem(p) / aleaf(p)

components/elm/src/biogeochem/CNEcosystemDynBetrMod.F90

@@ -59,7 +59,8 @@ module CNEcosystemDynBetrMod
   !-----------------------------------------------------------------------
   subroutine CNEcosystemDynBetr(bounds,                             &
          num_soilc, filter_soilc,                                        &
-         num_soilp, filter_soilp, num_pcropp, filter_pcropp, doalb,      &
+         num_soilp, filter_soilp, num_pcropp, filter_pcropp,             &
+         num_ppercropp, filter_ppercropp, doalb,                         &
          cnstate_vars, carbonflux_vars, carbonstate_vars,                &
          c13_carbonflux_vars, c13_carbonstate_vars,                      &
          c14_carbonflux_vars, c14_carbonstate_vars,                      &
@@ -118,6 +119,8 @@ subroutine CNEcosystemDynBetr(bounds,                             &
     integer                          , intent(in)    :: filter_soilp(:)   ! filter for soil patches
     integer                          , intent(in)    :: num_pcropp        ! number of prog. crop patches in filter
     integer                          , intent(in)    :: filter_pcropp(:)  ! filter for prognostic crop patches
+    integer                          , intent(in)    :: num_ppercropp     ! number of prog perennial crop patches in filter
+    integer                          , intent(in)    :: filter_ppercropp(:) ! filter for prognostic perennial crop patches
     logical                          , intent(in)    :: doalb             ! true = surface albedo calculation time step
     type(cnstate_type)               , intent(inout) :: cnstate_vars
     type(carbonflux_type)            , intent(inout) :: carbonflux_vars
@@ -189,7 +192,8 @@ subroutine CNEcosystemDynBetr(bounds,                             &
 
        call t_startf('MaintenanceResp')
        if (crop_prog) then
-          call NitrogenFert(bounds, num_soilc,filter_soilc, num_pcropp, filter_pcropp )
+          call NitrogenFert(bounds, num_soilc,filter_soilc, num_pcropp, filter_pcropp, &
+                            num_ppercropp, filter_ppercropp )
 
        end if
        call MaintenanceResp(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &

components/elm/src/biogeochem/CropType.F90

@@ -727,7 +727,7 @@ subroutine UpdateAccVars(this, bounds, temperature_vars)
   end subroutine UpdateAccVars
 
   !-----------------------------------------------------------------------
-  subroutine CropIncrementYear (this, num_pcropp, filter_pcropp)
+  subroutine CropIncrementYear (this, num_pcropp, filter_pcropp, num_ppercropp, filter_ppercropp)
     !
     ! !DESCRIPTION:
     ! Increment the crop year, if appropriate
@@ -741,6 +741,8 @@ subroutine CropIncrementYear (this, num_pcropp, filter_pcropp)
     class(crop_type) :: this
     integer , intent(in) :: num_pcropp       ! number of prog. crop patches in filter
     integer , intent(in) :: filter_pcropp(:) ! filter for prognostic crop patches
+    integer , intent(in) :: num_ppercropp    ! number of prog perennial crop patches in filter
+    integer , intent(in) :: filter_ppercropp(:) ! filter for prognostic perennial crop patches
     !
     ! !LOCAL VARIABLES:
     integer kyr   ! current year
@@ -762,6 +764,12 @@ subroutine CropIncrementYear (this, num_pcropp, filter_pcropp)
 
           this%nyrs_crop_active_patch(p) = this%nyrs_crop_active_patch(p) + 1
        end do
+
+       do fp = 1, num_ppercropp
+          p = filter_ppercropp(fp)
+
+          this%nyrs_crop_active_patch(p) = this%nyrs_crop_active_patch(p) + 1
+       end do
     end if
 
   end subroutine CropIncrementYear

components/elm/src/biogeochem/EcosystemDynMod.F90

@@ -287,6 +287,7 @@ subroutine EcosystemDynNoLeaching1(bounds,                          &
        num_soilc, filter_soilc,                                         &
        num_soilp, filter_soilp,                                         &
        num_pcropp, filter_pcropp,                                       &
+       num_ppercropp, filter_ppercropp,                                 &
        cnstate_vars,  atm2lnd_vars,           &
        canopystate_vars, soilstate_vars, crop_vars,   &
        ch4_vars, photosyns_vars,   &
@@ -328,6 +329,8 @@ subroutine EcosystemDynNoLeaching1(bounds,                          &
     integer                  , intent(in)    :: filter_soilp(:)   ! filter for soil patches
     integer                  , intent(in)    :: num_pcropp        ! number of prog. crop patches in filter
     integer                  , intent(in)    :: filter_pcropp(:)  ! filter for prognostic crop patches
+    integer                  , intent(in)    :: num_ppercropp     ! number of prog perennial crop patches in filter
+    integer                  , intent(in)    :: filter_ppercropp(:) ! filter for prognostic perennial crop patches
     type(cnstate_type)       , intent(inout) :: cnstate_vars
     type(atm2lnd_type)       , intent(in)    :: atm2lnd_vars
     type(canopystate_type)   , intent(in)    :: canopystate_vars
@@ -405,7 +408,8 @@ subroutine EcosystemDynNoLeaching1(bounds,                          &
        event = 'MaintenanceResp'
        call t_start_lnd(event)
        if (crop_prog) then
-          call NitrogenFert(bounds, num_soilc,filter_soilc, num_pcropp, filter_pcropp )
+          call NitrogenFert(bounds, num_soilc,filter_soilc, num_pcropp, filter_pcropp, &
+                            num_ppercropp, filter_ppercropp)
           call PhosphorusFert(bounds, num_soilc, filter_soilc )
 
           call CNSoyfix(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &

components/elm/src/biogeochem/NitrogenDynamicsMod.F90

@@ -387,7 +387,8 @@ subroutine NitrogenLeaching(num_soilc, filter_soilc, dt )
  end subroutine NitrogenLeaching
 
   !-----------------------------------------------------------------------
-  subroutine NitrogenFert(bounds, num_soilc, filter_soilc, num_pcropp, filter_pcropp)
+  subroutine NitrogenFert(bounds, num_soilc, filter_soilc, &
+                          num_pcropp, filter_pcropp, num_ppercropp, filter_ppercropp)
     !
     ! !DESCRIPTION:
     ! On the radiation time step, update the nitrogen fertilizer for crops
@@ -403,6 +404,8 @@ subroutine NitrogenFert(bounds, num_soilc, filter_soilc, num_pcropp, filter_pcro
     integer                 , intent(in)    :: filter_soilc(:) ! filter for soil columns
     integer , intent(in) :: num_pcropp       ! number of prog. crop patches in filter
     integer , intent(in) :: filter_pcropp(:) ! filter for prognostic crop patches
+    integer , intent(in) :: num_ppercropp     ! number of prog perennial crop patches in filter
+    integer , intent(in) :: filter_ppercropp(:) ! filter for prognostic perennial crop patches
     !
     ! !LOCAL VARIABLES:
     integer :: c,fc,p,fp                 ! indices
@@ -444,6 +447,10 @@ subroutine NitrogenFert(bounds, num_soilc, filter_soilc, num_pcropp, filter_pcro
             p = filter_pcropp(fp)
             totalfert(p) = synthfert(p) + manure(p)
          end do
+         do fp = 1, num_ppercropp
+            p = filter_ppercropp(fp)
+            totalfert(p) = synthfert(p) + manure(p)
+         end do
       end if
 
       ! if fan_to_bgc_crop == .true., FAN fills in the fert_to_sminn and totalfert for