Merge branch 'evasinha/lnd/add_winter_wheat_ylu' (PR #6986)

This feature adds new parameterizations to represent winter wheat crop in ELM.

The implementation is based on Yaqiong Lu’s implementation of winter wheat in CLM4.5. https://doi.org/10.5194/gmd-10-1873-2017

The model was calibrated using the methodology used in Sinha et al. 2022.
https://doi.org/10.1029/2022MS003171

[non BFB] for Crops

Author: peterdschwartz

cime_config/tests.py

@@ -72,7 +72,6 @@
             "ERS_D.f19_f19.IELM.elm-ic_f19_f19_ielm",
             "ERS_D.f09_g16.I1850ELMCN",
             "ERS_D.ne4pg2_oQU480.I20TRELM.elm-disableDynpftCheck",
-            "SMS_Ly2_P1x1_D.1x1_smallvilleIA.IELMCNCROP.elm-lulcc_sville",
             "ERS_D.f19_g16.I1850GSWCNPRDCTCBC.elm-ctc_f19_g16_I1850GSWCNPRDCTCBC",
             "ERS_D.f09_f09.IELM.elm-solar_rad",
             "ERS_D.f09_f09.IELM.elm-koch_snowflake",
@@ -87,12 +86,13 @@
             "ERS.f19_f19.I1850ELMCN",
             "ERS.f19_f19.I20TRELMCN",
             "SMS_Ld1.hcru_hcru.I1850CRUELMCN",
-            "SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-force_netcdf_pio",
+            "SMS_Ly1_P1x1.1x1_smallvilleIA.I20TRGSWCNPCROP.elm-lulcc_sville",
+            "SMS_Ly5_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-force_netcdf_pio",
             "ERS.f19_g16.I1850ELM.elm-betr",
             "ERS.f19_g16.I1850ELM.elm-vst",
             "ERS.f09_g16.I1850ELMCN.elm-bgcinterface",
             "SMS.r05_r05.I1850ELMCN.elm-qian_1948",
-            "SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-per_crop",
+            "SMS_Ly5_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-per_crop",
             "SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-fan",
             "SMS.r05_r05.IELM.elm-topounit",
             "SMS.r05_r05.IELM.elm-topounit_im2",

components/elm/cime_config/testdefs/testmods_dirs/elm/force_netcdf_pio/user_nl_elm

@@ -0,0 +1,4 @@
+hist_mfilt = 1, 1
+hist_nhtfrq = 0, 0
+hist_dov2xy = .true., .false.
+hist_fincl2 = 'GPP', 'ER', 'NEE', 'NPP', 'NBP', 'TOTSOMC', 'DMYIELD', 'PLANTDAY', 'HARVESTDAY'

components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/shell_commands

@@ -1,2 +1,3 @@
 ./xmlchange ELM_BLDNML_OPTS="-irrig .true." -append
 if [ `./xmlquery --value MACH` == chrysalis ]; then ./xmlchange FORCE_BUILD_SMP=TRUE; fi
+./xmlchange JOB_WALLCLOCK_TIME=2:00:00

components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm

@@ -2,3 +2,7 @@ flanduse_timeseries = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/landuse.timeseries_1x
 fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp50_simyr2000_c230905.nc'
 
 check_dynpft_consistency = .false.
+hist_mfilt = 1, 1
+hist_nhtfrq = 0, 0
+hist_dov2xy = .true., .false.
+hist_fincl2 = 'GPP', 'ER', 'NEE', 'NPP', 'NBP', 'TOTSOMC', 'DMYIELD', 'PLANTDAY', 'HARVESTDAY'

components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm

@@ -1,2 +1,6 @@
 fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp50_simyr2000_miscanthus_c230905.nc'
 do_budgets = .false.
+hist_mfilt = 1, 1
+hist_nhtfrq = 0, 0
+hist_dov2xy = .true., .false.
+hist_fincl2 = 'GPP', 'ER', 'NEE', 'NPP', 'NBP', 'TOTSOMC', 'DMYIELD', 'PLANTDAY', 'HARVESTDAY'

components/elm/src/biogeochem/AllocationMod.F90

@@ -386,7 +386,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
     use elm_varctl      , only : carbonphosphorus_only!
     use pftvarcon        , only: npcropmin, declfact, bfact, aleaff, arootf, astemf, noveg
     use pftvarcon        , only: arooti, fleafi, allconsl, allconss, grperc, grpnow, nsoybean
-    use pftvarcon        , only: iscft, percrop
+    use pftvarcon        , only: iscft, percrop, nwcereal, nwcerealirrig
     use elm_varpar       , only: nlevdecomp
     use elm_varcon       , only: nitrif_n2o_loss_frac, secspday
     !
@@ -423,6 +423,8 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
     !! Local P variables
     real(r8):: cpl,cpfr,cplw,cpdw,cpg                                    !C:N ratios for leaf, fine root, and wood
     real(r8):: puptake_prof(bounds%begc:bounds%endc, 1:nlevdecomp)
+    integer, parameter :: cphase_gf = 3                                  !Crop phenology phase grain fill
+    integer, parameter :: max_lai   = 1                                  !Maximum allowed lai
 
 
   !-----------------------------------------------------------------------
@@ -456,6 +458,8 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
 
          hui                          => crop_vars%gddplant_patch                              , & ! Input:  [real(r8) (:)   ]  =gdd since planting (gddplant)
          leafout                      => crop_vars%gddtsoi_patch                               , & ! Input:  [real(r8) (:)   ]  =gdd from top soil layer temperature
+         vf                           => crop_vars%vf_patch                                    , & ! Output: [real(r8) (:)   ]  vernalization factor
+         cphase                       => crop_vars%cphase_patch                                , & ! Output: [real(r8) (:)   ]  phenology phase
 
          xsmrpool                     => veg_cs%xsmrpool                       , & ! Input:  [real(r8) (:)   ]  (gC/m2) temporary photosynthate C pool
          leafc                        => veg_cs%leafc                          , & ! Input:  [real(r8) (:)   ]
@@ -689,7 +693,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
                   ! allocation rules for crops based on maturity and linear decrease
                   ! of amount allocated to roots over course of the growing season
 
-                  if (peaklai(p) == 1) then ! lai at maximum allowed
+                  if (peaklai(p) == max_lai) then ! lai at maximum allowed
                      arepr(p) = 0._r8
                      aleaf(p) = 1.e-5_r8
                      aroot(p) = max(0._r8, min(1._r8, arooti(ivt(p)) -   &
@@ -717,6 +721,13 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
                   aleafi(p) = aleaf(p) ! to reproductive phenology stage begins
                   grain_flag(p) = 0._r8 ! setting to 0 while in phase 2
 
+                  ! Added based on Yaqiong Lu et al., 2017 in Geosci. Model Dev.
+                  ! when peaklai==1, astem=0 and then astemi=0, so the astem in phase 3 will
+                  ! equal to 0 and therefore resulted a very large arepr and grainc
+                  if(peaklai(p)==max_lai .and. (ivt(p) == nwcereal .or. ivt(p) == nwcerealirrig)) then
+                    astemi(p)=0.8_r8
+                  end if
+
                   ! Phase 2 completed:
                   ! ==================
                   ! shift allocation either when enough gdd are accumulated or maximum number
@@ -757,7 +768,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
                   !would be bypassed altogether, not the intended outcome. I checked several of my output files and
                   !they all seemed to be going through the retranslocation loop for soybean - good news.
 
-                  if (ivt(p) /= nsoybean .or. astem(p) == astemf(ivt(p)) .or. peaklai(p) == 1._r8) then
+                  if (ivt(p) /= nsoybean .or. astem(p) == astemf(ivt(p)) .or. peaklai(p) == max_lai) then
                      if (grain_flag(p) == 0._r8) then
                         t1 = 1 / dt
                         leafn_to_retransn(p) = t1 * ((leafc(p) / leafcn(ivt(p))) - (leafc(p) / &
@@ -775,6 +786,12 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp
 
                   arepr(p) = 1._r8 - aroot(p) - astem(p) - aleaf(p)
 
+                  ! Added based on Yaqiong Lu et al., 2017 in Geosci. Model Dev.
+                  if(cphase(p) == cphase_gf .and. (ivt(p) == nwcereal .or. ivt(p) == nwcerealirrig)) then
+                     arepr(p) = arepr(p) * vf(p)
+                     aroot(p) = 1._r8 - aleaf(p) - astem(p) - arepr(p)
+                  end if
+
                else                   ! pre emergence
                   aleaf(p) = 1.e-5_r8 ! allocation coefficients should be irrelevant
                   astem(p) = 0._r8    ! because crops have no live carbon pools;

components/elm/src/biogeochem/CropType.F90

@@ -38,6 +38,16 @@ module CropType
      logical , pointer :: croplive_patch          (:)   ! patch Flag, true if planted, not harvested
      logical , pointer :: cropplant_patch         (:)   ! patch Flag, true if planted
      integer , pointer :: harvdate_patch          (:)   ! patch harvest date
+     real(r8), pointer :: rateh_patch             (:)   ! increase of tolerance caused by cold hardening index
+     real(r8), pointer :: rated_patch             (:)   ! loss of tolerance caused by dehardening
+     real(r8), pointer :: rates_patch             (:)   ! loss of tolerance caused by low temperature
+     real(r8), pointer :: rater_patch             (:)   ! loss of tolerance caused by respiration under snow
+     real(r8), pointer :: lt50_patch              (:)   ! the lethal temperature at which 50% of the individuals are damaged
+     real(r8), pointer :: fsurv_patch             (:)   ! winter wheat survival rate
+     real(r8), pointer :: accfsurv_patch          (:)   ! accumulated winter wheat survival rate
+     real(r8), pointer :: countfsurv_patch        (:)   ! count of accumulated winter wheat survival rate
+     real(r8), pointer :: wdd_patch               (:)   ! winter wheat weighted cumulated degree days
+     real(r8), pointer :: tcrown_patch            (:)   ! crown temperature
      real(r8), pointer :: fertnitro_patch         (:)   ! patch fertilizer nitrogen
      real(r8), pointer :: fertphosp_patch         (:)   ! patch fertilizer phosphorus
      real(r8), pointer :: gddplant_patch          (:)   ! patch accum gdd past planting date for crop       (ddays)
@@ -128,6 +138,16 @@ subroutine InitAllocate(this, bounds)
     allocate(this%croplive_patch         (begp:endp)) ; this%croplive_patch         (:) = .false.
     allocate(this%cropplant_patch        (begp:endp)) ; this%cropplant_patch        (:) = .false.
     allocate(this%harvdate_patch         (begp:endp)) ; this%harvdate_patch         (:) = huge(1) 
+    allocate(this%rateh_patch            (begp:endp)) ; this%rateh_patch            (:) = spval
+    allocate(this%rated_patch            (begp:endp)) ; this%rated_patch            (:) = spval
+    allocate(this%rates_patch            (begp:endp)) ; this%rates_patch            (:) = spval
+    allocate(this%rater_patch            (begp:endp)) ; this%rater_patch            (:) = spval
+    allocate(this%lt50_patch             (begp:endp)) ; this%lt50_patch             (:) = spval
+    allocate(this%fsurv_patch            (begp:endp)) ; this%fsurv_patch            (:) = spval
+    allocate(this%accfsurv_patch         (begp:endp)) ; this%accfsurv_patch         (:) = spval
+    allocate(this%countfsurv_patch       (begp:endp)) ; this%countfsurv_patch       (:) = spval
+    allocate(this%wdd_patch              (begp:endp)) ; this%wdd_patch              (:) = spval
+    allocate(this%tcrown_patch           (begp:endp)) ; this%tcrown_patch           (:) = spval
     allocate(this%fertnitro_patch        (begp:endp)) ; this%fertnitro_patch        (:) = spval
     allocate(this%fertphosp_patch        (begp:endp)) ; this%fertphosp_patch        (:) = spval
     allocate(this%gddplant_patch         (begp:endp)) ; this%gddplant_patch         (:) = spval
@@ -669,7 +689,8 @@ subroutine UpdateAccVars(this, bounds, temperature_vars)
           rbufslp(p) = max(0._r8, min(mxtmp(ivt), &
                veg_es%t_ref2m(p)-(SHR_CONST_TKFRZ + baset(ivt)))) &
                * dtime/SHR_CONST_CDAY
-          if (ivt == nwcereal .or. ivt == nwcerealirrig) then
+          ! Modified based on Yaqiong Lu et al., 2017 in Geosci. Model Dev.
+          if ((ivt == nwcereal .or. ivt == nwcerealirrig) .and. this%cphase_patch(p) > 1) then
              rbufslp(p) = rbufslp(p)*this%vf_patch(p)
           end if
        else
@@ -691,9 +712,9 @@ subroutine UpdateAccVars(this, bounds, temperature_vars)
                ((col_es%t_soisno(c,1)*col_pp%dz(c,1) + &
                col_es%t_soisno(c,2)*col_pp%dz(c,2))/(col_pp%dz(c,1)+col_pp%dz(c,2))) - &
                (SHR_CONST_TKFRZ + baset(ivt)))) * dtime/SHR_CONST_CDAY
-          if (ivt == nwcereal .or. ivt == nwcerealirrig) then
-             rbufslp(p) = rbufslp(p)*this%vf_patch(p)
-          end if
+          ! Removed rbufslp modification based on Yaqiong Lu et al., 2017 in Geosci. Model Dev.
+          ! Removed the vf control on gddtsoil, because the vernalization
+          ! occurs after leaf emerge and end at flowering
        else
           rbufslp(p) = accumResetVal
        end if