Merge branch 'refactor_cnallocation'

Refactor NutrientCompetition / CNAllocation to provide hooks for AgSys Major refactor of NutrientCompetition / CNAllocation to provide hooks for AgSys crop model: separates the NutrientCompetition modules into pieces based on (1) consolidating duplicate code between the Clm45 and FlexibleCN versions, and (2) separating pieces that will vs. won't be used for crop patches when running with the upcoming AgSys crop model. I have restored the old CNAllocationMod, with some of the responsibilities that it used to have. (I'm not sure it's appropriate to have the calculation of gpp and maint resp in CNAllocationMod, but I left it there because it has always been combined with the allocation code, including back when we had a separate CNAllocationMod.)
samsrabin · Apr 30, 2022 · e2932c9 · e2932c9
2 parents 82a63cc + 6ec5166
commit e2932c9
Show file tree

Hide file tree

Showing 15 changed files with 1,655 additions and 1,319 deletions.
diff --git a/doc/ChangeLog b/doc/ChangeLog
@@ -1,4 +1,82 @@
 ===============================================================
+Tag name: ctsm5.1.dev092
+Originator(s): sacks (Bill Sacks)
+Date: Fri Apr 29 18:31:48 MDT 2022
+One-line Summary: Refactor NutrientCompetition / CNAllocation to provide hooks for AgSys
+
+Purpose and description of changes
+----------------------------------
+
+Major refactor of NutrientCompetition / CNAllocation to provide hooks
+for AgSys crop model: separates the NutrientCompetition modules into
+pieces based on (1) consolidating duplicate code between the Clm45 and
+FlexibleCN versions, and (2) separating pieces that will vs. won't be
+used for crop patches when running with the upcoming AgSys crop model.
+
+I have restored the old CNAllocationMod, with some of the
+responsibilities that it used to have. (I'm not sure it's appropriate to
+have the calculation of gpp and maint resp in CNAllocationMod, but I
+left it there because it has always been combined with the allocation
+code, including back when we had a separate CNAllocationMod.)
+
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+(Details of any changes will be given in the "Answer changes" section below.)
+
+    [Put an [X] in the box for any configuration with significant answer changes.]
+
+[ ] clm5_1
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+
+Testing summary:
+----------------
+
+  regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing):
+
+    cheyenne ---- PASS
+    izumi ------- PASS
+
+
+Answer changes
+--------------
+
+Changes answers relative to baseline: NO - though potential for answer
+changes in unusual cases
+
+   One change in this tag has the potential for answer changes in
+   unusual cases, even though no answer changes were observed in
+   testing: Previously, only the FlexibleCN code applied some logic that
+   changed crop allocation fractions during the grainfill period if
+   peaklai had been reached. I have changed this so that this logic is
+   applied with or without FlexibleCN (because I have moved this block
+   of code to the CNAllocation module, which is shared between the
+   FlexibleCN and non-FlexibleCN versions). I thought this would change
+   answers, but it appears not to, at least based on the tests in the
+   test suite as well as an extra 5-year test I did at f19 resolution
+   (ERS_Ly5.f19_g17.IHistClm45BgcCrop.cheyenne_intel.clm-cropMonthOutput).
+   It's possible that this changes answers in rare cases or with an
+   unusual combination of options that we don't test: specifically, it
+   might change answers for the atypical situation where you are running
+   with FUN but not FlexibleCN. (In this case, I believe this change is
+   the correct thing to do.)
+
+
+Other details
+-------------
+Pull Requests that document the changes (include PR ids):
+https://github.com/ESCOMP/CTSM/pull/1705
+
+===============================================================
+===============================================================
 Tag name: ctsm5.1.dev091
 Originator(s): rgknox (Ryan Knox,rgknox@lbl.gov)
 Date: Fri Apr 22 14:11:50 EDT 2022

diff --git a/doc/ChangeSum b/doc/ChangeSum
@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+    ctsm5.1.dev092    sacks 04/29/2022 Refactor NutrientCompetition / CNAllocation to provide hooks for AgSys
     ctsm5.1.dev091   rgknox 04/22/2022 clm decomp method is now passed to fates to enabled mimics coupling
     ctsm5.1.dev090 samrabin 03/31/2022 Fix misleading name of "gddplant"
     ctsm5.1.dev089    sacks 03/31/2022 For CLM45 apply peaklai to aleaf in grainfill

diff --git a/src/biogeochem/CNAllocationMod.F90 b/src/biogeochem/CNAllocationMod.F90
diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90
@@ -11,7 +11,7 @@ module CNDriverMod
   use decompMod                       , only : bounds_type
   use perf_mod                        , only : t_startf, t_stopf
   use clm_varctl                      , only : use_nitrif_denitrif, use_nguardrail
-  use clm_varctl                      , only : iulog, use_crop
+  use clm_varctl                      , only : iulog, use_crop, use_crop_agsys
   use SoilBiogeochemDecompCascadeConType, only : mimics_decomp, century_decomp, decomp_method
   use CNSharedParamsMod               , only : use_fun
   use CNVegStateType                  , only : cnveg_state_type
@@ -42,7 +42,6 @@ module CNDriverMod
   use SaturatedExcessRunoffMod        , only : saturated_excess_runoff_type
   use ActiveLayerMod                  , only : active_layer_type
   use SoilWaterRetentionCurveMod      , only : soil_water_retention_curve_type
-  use CropReprPoolsMod                    , only : nrepr
   !
   ! !PUBLIC TYPES:
   implicit none
@@ -83,7 +82,8 @@ end subroutine CNDriverInit
 
   !-----------------------------------------------------------------------
   subroutine CNDriverNoLeaching(bounds,                                                    &
-       num_soilc, filter_soilc, num_soilp, filter_soilp, num_pcropp, filter_pcropp,        &
+       num_soilc, filter_soilc, num_soilp, filter_soilp,                                   &
+       num_pcropp, filter_pcropp, num_soilnopcropp, filter_soilnopcropp,                   &
        num_exposedvegp, filter_exposedvegp, num_noexposedvegp, filter_noexposedvegp,       &
        cnveg_state_inst,                                                                   &
        cnveg_carbonflux_inst, cnveg_carbonstate_inst,                                      &
@@ -114,6 +114,8 @@ subroutine CNDriverNoLeaching(bounds,
     use clm_varpar                        , only: nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools
     use subgridAveMod                     , only: p2c
     use CropType                          , only: crop_type
+    use CNAllocationMod                   , only: calc_gpp_mr_availc, calc_crop_allocation_fractions
+    use CNAllocationMod                   , only: calc_allometry
     use CNNDynamicsMod                    , only: CNNDeposition,CNNFixation, CNNFert, CNSoyfix,CNFreeLivingFixation
     use CNMRespMod                        , only: CNMResp
     use CNFUNMod                          , only: CNFUNInit  !, CNFUN 
@@ -151,6 +153,8 @@ subroutine CNDriverNoLeaching(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_soilnopcropp       ! number of non-prog. crop soil patches in filter
+    integer                                 , intent(in)    :: filter_soilnopcropp(:) ! filter for non-prog. crop soil patches
     integer                                 , intent(in)    :: num_exposedvegp        ! number of points in filter_exposedvegp
     integer                                 , intent(in)    :: filter_exposedvegp(:)  ! patch filter for non-snow-covered veg
     integer                                 , intent(in)    :: num_noexposedvegp       ! number of points in filter_noexposedvegp
@@ -203,8 +207,6 @@ subroutine CNDriverNoLeaching(bounds,
     real(r8):: pmnf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) !potential mineral N flux, from one pool to another
     real(r8):: p_decomp_npool_to_din(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions)  ! potential flux to dissolved inorganic N
     real(r8):: p_decomp_cn_gain(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools)  ! C:N ratio of the flux gained by the receiver pool
-    real(r8):: arepr(bounds%begp:bounds%endp,nrepr) ! reproductive allocation coefficient(s) (only used for use_crop)
-    real(r8):: aroot(bounds%begp:bounds%endp) ! root allocation coefficient (only used for use_crop)
     integer :: begp,endp
     integer :: begc,endc
 
@@ -386,16 +388,44 @@ subroutine CNDriverNoLeaching(bounds,
 
      end if
 
+     call t_startf('cnalloc')
+     call calc_gpp_mr_availc( &
+          bounds, num_soilp, filter_soilp, &
+          crop_inst, photosyns_inst, canopystate_inst, &
+          cnveg_carbonstate_inst, cnveg_carbonflux_inst, &
+          c13_cnveg_carbonflux_inst, c14_cnveg_carbonflux_inst)
+
+     if (.not. use_crop_agsys) then
+        call calc_crop_allocation_fractions(bounds, num_pcropp, filter_pcropp, &
+             crop_inst, cnveg_state_inst)
+     end if
+
+     call calc_allometry(num_soilp, filter_soilp, &
+          cnveg_carbonflux_inst, cnveg_state_inst)
+     call t_stopf('cnalloc')
+
      call t_startf('calc_plant_nutrient_demand')
+     ! We always call calc_plant_nutrient_demand for natural veg patches, but only call
+     ! it for crop patches if NOT running with AgSys (since AgSys calculates the relevant
+     ! output variables in its own way).
      call nutrient_competition_method%calc_plant_nutrient_demand ( &
-         bounds, num_soilp, filter_soilp,                                 &
-         photosyns_inst, crop_inst, canopystate_inst,                     &
-         cnveg_state_inst, cnveg_carbonstate_inst, cnveg_carbonflux_inst, &
-         c13_cnveg_carbonflux_inst, c14_cnveg_carbonflux_inst,            &
-         cnveg_nitrogenstate_inst, cnveg_nitrogenflux_inst,               &
-         soilbiogeochem_carbonflux_inst, soilbiogeochem_nitrogenstate_inst, &
-         energyflux_inst, &
-         aroot=aroot(begp:endp), arepr=arepr(begp:endp,:))
+          bounds,                                                          &
+          num_soilnopcropp, filter_soilnopcropp, .false.,                  &
+          crop_inst, canopystate_inst,                                     &
+          cnveg_state_inst, cnveg_carbonstate_inst, cnveg_carbonflux_inst, &
+          cnveg_nitrogenstate_inst, cnveg_nitrogenflux_inst,               &
+          soilbiogeochem_carbonflux_inst, soilbiogeochem_nitrogenstate_inst, &
+          energyflux_inst)
+     if (.not. use_crop_agsys) then
+        call nutrient_competition_method%calc_plant_nutrient_demand ( &
+             bounds,                                                          &
+             num_pcropp, filter_pcropp, .true.,                               &
+             crop_inst, canopystate_inst,                                     &
+             cnveg_state_inst, cnveg_carbonstate_inst, cnveg_carbonflux_inst, &
+             cnveg_nitrogenstate_inst, cnveg_nitrogenflux_inst,               &
+             soilbiogeochem_carbonflux_inst, soilbiogeochem_nitrogenstate_inst, &
+             energyflux_inst)
+     end if
 
      ! get the column-averaged plant_ndemand (needed for following call to SoilBiogeochemCompetition)
 
@@ -429,8 +459,6 @@ subroutine CNDriverNoLeaching(bounds,
          c13_cnveg_carbonflux_inst, c14_cnveg_carbonflux_inst, &
          cnveg_nitrogenstate_inst, cnveg_nitrogenflux_inst, &
          soilbiogeochem_nitrogenstate_inst, &
-         aroot=aroot(begp:endp), &
-         arepr=arepr(begp:endp,:), &
          fpg_col=soilbiogeochem_state_inst%fpg_col(begc:endc))
     call t_stopf('calc_plant_nutrient_competition')
 

diff --git a/src/biogeochem/CNPhenologyMod.F90 b/src/biogeochem/CNPhenologyMod.F90
@@ -28,6 +28,8 @@ module CNPhenologyMod
   use CNVegnitrogenstateType          , only : cnveg_nitrogenstate_type
   use CNVegnitrogenfluxType           , only : cnveg_nitrogenflux_type
   use CropType                        , only : crop_type
+  use CropType                        , only : cphase_planted, cphase_leafemerge
+  use CropType                        , only : cphase_grainfill, cphase_harvest
   use pftconMod                       , only : pftcon
   use SoilStateType                   , only : soilstate_type
   use TemperatureType                 , only : temperature_type
@@ -48,6 +50,7 @@ module CNPhenologyMod
   public :: CNPhenologyreadNML   ! Read namelist
   public :: CNPhenologyInit      ! Initialization
   public :: CNPhenology          ! Update
+  public :: CropPhase            ! Get the current phase of each crop patch
 
   ! !PUBLIC for unit testing
   public :: CNPhenologySetNML         ! Set the namelist setttings explicitly for unit tests
@@ -2012,7 +2015,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
          offset_flag(p) = 0._r8 ! carbon and nitrogen transfers
 
          if (croplive(p)) then
-            cphase(p) = 1._r8
+            cphase(p) = cphase_planted
 
             ! call vernalization if winter temperate cereal planted, living, and the
             ! vernalization factor is not 1;
@@ -2045,8 +2048,14 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
                hui(p) = max(hui(p),huigrain(p))
             endif
 
+            ! The following conditionals are similar to those in CropPhase. However, they
+            ! differ slightly because here we are potentially setting a new crop phase,
+            ! whereas CropPhase is just designed to get the current, already-determined
+            ! phase. However, despite these differences: if you make changes to the
+            ! following conditionals, you should also check to see if you should make
+            ! similar changes in CropPhase.
             if (leafout(p) >= huileaf(p) .and. hui(p) < huigrain(p) .and. idpp < mxmat(ivt(p))) then
-               cphase(p) = 2._r8
+               cphase(p) = cphase_leafemerge
                if (abs(onset_counter(p)) > 1.e-6_r8) then
                   onset_flag(p)    = 1._r8
                   onset_counter(p) = dt
@@ -2076,7 +2085,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
                harvest_count(p) = harvest_count(p) + 1
                crop_inst%hdates_thisyr(p, harvest_count(p)) = real(jday, r8)
                croplive(p) = .false.     ! no re-entry in greater if-block
-               cphase(p) = 4._r8
+               cphase(p) = cphase_harvest
                if (tlai(p) > 0._r8) then ! plant had emerged before harvest
                   offset_flag(p) = 1._r8
                   offset_counter(p) = dt
@@ -2106,7 +2115,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
                ! Use CN's simple formula at least as a place holder (slevis)
 
             else if (hui(p) >= huigrain(p)) then
-               cphase(p) = 3._r8
+               cphase(p) = cphase_grainfill
                bglfr(p) = 1._r8/(leaf_long(ivt(p))*avg_dayspyr*secspday)
             end if
 
@@ -2144,6 +2153,67 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
 
   end subroutine CropPhenology
 
+  !-----------------------------------------------------------------------
+  subroutine CropPhase(bounds, num_pcropp, filter_pcropp, &
+       crop_inst, cnveg_state_inst, crop_phase)
+    !
+    ! !DESCRIPTION:
+    ! Get the current phase of each crop patch.
+    !
+    ! The returned values (in crop_phase) are from the set of cphase_* values defined in
+    ! CropType. The returned values in crop_phase are only valid for patches where
+    ! croplive is true; the values are undefined where croplive is false and should not be
+    ! used there!
+    !
+    ! This has logic similar to that in CropPhenology. If you make changes here, you
+    ! should also check if similar changes need to be made in CropPhenology.
+    !
+    ! !ARGUMENTS:
+    type(bounds_type)      , intent(in)    :: bounds
+    integer                , intent(in)    :: num_pcropp       ! number of prog crop patches in filter
+    integer                , intent(in)    :: filter_pcropp(:) ! filter for prognostic crop patches
+    type(crop_type)        , intent(in)    :: crop_inst
+    type(cnveg_state_type) , intent(in)    :: cnveg_state_inst
+    real(r8)               , intent(inout) :: crop_phase(bounds%begp:)
+    !
+    ! !LOCAL VARIABLES:
+    integer :: p, fp
+
+    character(len=*), parameter :: subname = 'CropPhase'
+    !-----------------------------------------------------------------------
+    SHR_ASSERT_ALL_FL((ubound(crop_phase) == [bounds%endp]), sourcefile, __LINE__)
+
+    associate( &
+         croplive =>    crop_inst%croplive_patch        , & ! Input: [logical  (:) ]  Flag, true if planted, not harvested
+         hui      =>    crop_inst%hui_patch             , & ! Input: [real(r8) (:) ]  gdd since planting (gddplant)
+         leafout  =>    crop_inst%gddtsoi_patch         , & ! Input: [real(r8) (:) ]  gdd from top soil layer temperature
+         huileaf  =>    cnveg_state_inst%huileaf_patch  , & ! Input: [real(r8) (:) ]  heat unit index needed from planting to leaf emergence
+         huigrain =>    cnveg_state_inst%huigrain_patch   & ! Input: [real(r8) (:) ]  same to reach vegetative maturity
+         )
+
+    do fp = 1, num_pcropp
+       p = filter_pcropp(fp)
+
+       if (croplive(p)) then
+          ! Start with cphase_planted, but this might get changed in the later
+          ! conditional blocks.
+          crop_phase(p) = cphase_planted
+          if (leafout(p) >= huileaf(p) .and. hui(p) < huigrain(p)) then
+             crop_phase(p) = cphase_leafemerge
+          else if (hui(p) >= huigrain(p)) then
+             ! Since we know croplive is true, any hui greater than huigrain implies that
+             ! we're in the grainfill stage: if we were passt gddmaturity then croplive
+             ! would be false.
+             crop_phase(p) = cphase_grainfill
+          end if
+       end if
+    end do
+
+    end associate
+
+  end subroutine CropPhase
+
+
   !-----------------------------------------------------------------------
   subroutine CropPhenologyInit(bounds)
     !