diff --git a/cime_config/tests.py b/cime_config/tests.py index 1b619092f2a2..3f74472549a1 100644 --- a/cime_config/tests.py +++ b/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", diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/force_netcdf_pio/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/force_netcdf_pio/user_nl_elm new file mode 100644 index 000000000000..521985653fc8 --- /dev/null +++ b/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' diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/shell_commands b/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/shell_commands index 98f1754a6b3c..f4b9806a903d 100755 --- a/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/shell_commands +++ b/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 diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm index 55b35406dcbf..43bde325e639 100644 --- a/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm +++ b/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' diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm index 6554c076a27a..f29e1ca7093b 100644 --- a/components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm +++ b/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' diff --git a/components/elm/src/biogeochem/AllocationMod.F90 b/components/elm/src/biogeochem/AllocationMod.F90 index ba14dbc352db..591a4bf8f692 100644 --- a/components/elm/src/biogeochem/AllocationMod.F90 +++ b/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; diff --git a/components/elm/src/biogeochem/CropType.F90 b/components/elm/src/biogeochem/CropType.F90 index 8394d21ca348..1c5e6b36ad44 100644 --- a/components/elm/src/biogeochem/CropType.F90 +++ b/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 diff --git a/components/elm/src/biogeochem/PhenologyMod.F90 b/components/elm/src/biogeochem/PhenologyMod.F90 index c46e251b60cd..70f747f06673 100644 --- a/components/elm/src/biogeochem/PhenologyMod.F90 +++ b/components/elm/src/biogeochem/PhenologyMod.F90 @@ -1469,6 +1469,15 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & huigrain => cnstate_vars%huigrain_patch , & ! Output: [real(r8) (:) ] same to reach vegetative maturity cumvd => cnstate_vars%cumvd_patch , & ! Output: [real(r8) (:) ] cumulative vernalization d?ependence? hdidx => cnstate_vars%hdidx_patch , & ! Output: [real(r8) (:) ] cold hardening index? + lt50 => crop_vars%lt50_patch , & ! Output: [integer (:) ] the lethal temperature at which 50% of the individuals are damaged + wdd => crop_vars%wdd_patch , & ! Output: [integer (:) ] winter wheat weighted cumulated degree days + rateh => crop_vars%rateh_patch , & ! Output: [integer (:) ] increase of tolerance caused by cold hardening index + rated => crop_vars%rated_patch , & ! Output: [integer (:) ] loss of tolerance caused by dehardening + rates => crop_vars%rates_patch , & ! Output: [integer (:) ] loss of tolerance caused by low temperature + rater => crop_vars%rater_patch , & ! Output: [integer (:) ] loss of tolerance caused by respiration under snow + fsurv => crop_vars%fsurv_patch , & ! Output: [integer (:) ] winter wheat survival rate + accfsurv => crop_vars%accfsurv_patch , & ! Output: [integer (:) ] accumulated winter wheat survival rate + countfsurv => crop_vars%countfsurv_patch , & ! Output: [integer (:) ] count of accumulated winter wheat survival rate vf => crop_vars%vf_patch , & ! Output: [real(r8) (:) ] vernalization factor bglfr_leaf => cnstate_vars%bglfr_leaf_patch , & ! Output: [real(r8) (:) ] background leaf litterfall rate (1/s) bglfr_froot => cnstate_vars%bglfr_froot_patch , & ! Output: [real(r8) (:) ] background fine root litterfall rate (1/s) @@ -1501,6 +1510,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & plantmonth => crop_vars%plantmonth_patch , & ! Output: [real(r8) ):)] plant month plantday => crop_vars%plantday_patch , & ! Output: [real(r8) ):)] plant day harvday => crop_vars%harvday_patch , & ! Ouptut: [real(r8) ):)] harvest day + cphase => crop_vars%cphase_patch , & ! Output: [real(r8) (:)] phenology phase forc_rain => top_af%rain , & ! Input: [real(r8) (:)] rainfall rate wf2 => col_ws%wf2 & ! Output: [real(r8) (:)] soil water as frac. of whc for top 0.17 m ) @@ -1615,6 +1625,15 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & cumvd(p) = 0._r8 hdidx(p) = 0._r8 + lt50(p) = -5._r8 + wdd(p) = 0._r8 + rateh(p) = 0._r8 + rated(p) = 0._r8 + rater(p) = 0._r8 + rates(p) = 0._r8 + fsurv(p) = 1._r8 + accfsurv(p) = 1._r8 + countfsurv(p) = 1._r8 vf(p) = 0._r8 croplive(p) = .true. cropplant(p) = .true. @@ -1640,6 +1659,15 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & cumvd(p) = 0._r8 hdidx(p) = 0._r8 + lt50(p) = -5._r8 + wdd(p) = 0._r8 + rateh(p) = 0._r8 + rated(p) = 0._r8 + rater(p) = 0._r8 + rates(p) = 0._r8 + fsurv(p) = 1._r8 + accfsurv(p) = 1._r8 + countfsurv(p) = 1._r8 vf(p) = 0._r8 croplive(p) = .true. cropplant(p) = .true. @@ -1838,13 +1866,19 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & offset_flag(p) = 0._r8 ! carbon and nitrogen transfers if (croplive(p)) then + cphase(p) = 1._r8 ! call vernalization if winter temperate cereal planted, living, and the ! vernalization factor is not 1; ! vf affects the calculation of gddtsoi & gddplant - if (t_ref2m_min(p) < 1.e30_r8 .and. vf(p) /= 1._r8 .and. (ivt(p) == nwcereal .or. ivt(p) == nwcerealirrig)) then - call vernalization(p,canopystate_vars, cnstate_vars, crop_vars) + ! Modifications based on Yaqiong Lu et al., 2017 in Geosci. Model Dev. + if (vf(p) /= 1._r8 .and. (ivt(p) == nwcereal .or. ivt(p) == nwcerealirrig) .and. hui(p) < 0.7_r8*huigrain(p)) then + call vernalization(p, cnstate_vars, crop_vars) + end if + + if (ivt(p) == nwcereal .or. ivt(p) == nwcerealirrig) then + call coldtolerance(p, cnstate_vars, crop_vars) end if ! days past planting may determine harvest @@ -1864,6 +1898,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & ! transfer seed carbon to leaf emergence if (leafout(p) >= huileaf(p) .and. hui(p) < huigrain(p) .and. idpp < mxmat(ivt(p))) then + cphase(p) = 2._r8 if (abs(onset_counter(p)) > 1.e-6_r8) then onset_flag(p) = 1._r8 onset_counter(p) = dt @@ -1895,6 +1930,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & if (harvdate(p) >= NOT_Harvested) harvdate(p) = jday if (harvday(p) >= NOT_Harvested) harvday(p) = jday croplive(p) = .false. ! no re-entry in greater if-block + cphase(p) = 4._r8 if (tlai(p) > 0._r8) then ! plant had emerged before harvest offset_flag(p) = 1._r8 offset_counter(p) = dt @@ -1914,6 +1950,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 bglfr_leaf(p) = 1._r8/(leaf_long(ivt(p))*dayspyr*secspday) bglfr_froot(p) = 1._r8/(froot_long(ivt(p))*dayspyr*secspday) end if @@ -2244,9 +2281,7 @@ subroutine CropPhenologyInit(bounds) end subroutine CropPhenologyInit !----------------------------------------------------------------------- - subroutine vernalization(p, & - canopystate_vars,cnstate_vars, & - crop_vars) + subroutine vernalization(p, cnstate_vars, crop_vars) ! ! !DESCRIPTION: ! @@ -2260,32 +2295,26 @@ subroutine vernalization(p, & ! ! !ARGUMENTS: !$acc routine seq - integer , intent(in) :: p ! PATCH index running over - type(canopystate_type) , intent(in) :: canopystate_vars + integer , intent(in) :: p ! PATCH index running over type(cnstate_type) , intent(inout) :: cnstate_vars type(crop_type) , intent(inout) :: crop_vars ! ! LOCAL VARAIBLES: - real(r8) tcrown ! ? - real(r8) vd, vd1, vd2 ! vernalization dependence - real(r8) tkil ! Freeze kill threshold - integer c,g ! indices + real(r8), parameter :: vtmin = -1.3_r8 ! vernalization minimum temp based on Lu et al., 2017 + real(r8), parameter :: vtopt = 4.9_r8 ! vernalization optimum temp based on Lu et al., 2017 + real(r8), parameter :: vtmax = 15.7_r8 ! vernalization maximum temp based on Lu et al., 2017 + real(r8) alpha ! parameter in calculating vernalization rate + real(r8) tc ! t_ref2m in degree C + integer c ! indices + real(r8) dt_hr ! convert dt from sec to hour !------------------------------------------------------------------------ associate( & - tlai => canopystate_vars%tlai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index, no burying by snow - - t_ref2m => veg_es%t_ref2m , & ! Input: [real(r8) (:) ] 2 m height surface air temperature (K) - t_ref2m_min => veg_es%t_ref2m_min , & ! Input: [real(r8) (:) ] daily minimum of average 2 m height surface air temperature (K) - t_ref2m_max => veg_es%t_ref2m_max , & ! Input: [real(r8) (:) ] daily maximum of average 2 m height surface air temperature (K) - - snow_depth => col_ws%snow_depth , & ! Input: [real(r8) (:) ] snow height (m) - - hdidx => cnstate_vars%hdidx_patch , & ! Output: [real(r8) (:) ] cold hardening index? + tcrown => crop_vars%tcrown_patch , & ! Output: [integer (:) ] crown temperature + t_ref2m => veg_es%t_ref2m , & ! Input: [real(r8) (:) ] 2 m height surface air temperature (K) + snow_depth => col_ws%snow_depth , & ! Input: [real(r8) (:) ] snow height (m) cumvd => cnstate_vars%cumvd_patch , & ! Output: [real(r8) (:) ] cumulative vernalization d?ependence? - vf => crop_vars%vf_patch , & ! Output: [real(r8) (:) ] vernalization factor for cereal - gddmaturity => cnstate_vars%gddmaturity_patch , & ! Output: [real(r8) (:) ] gdd needed to harvest - huigrain => cnstate_vars%huigrain_patch & ! Output: [real(r8) (:) ] heat unit index needed to reach vegetative maturity + vf => crop_vars%vf_patch & ! Output: [real(r8) (:) ] vernalization factor for cereal ) c = veg_pp%column(p) @@ -2305,85 +2334,161 @@ subroutine vernalization(p, & ! if vf(p) = 1. then plant is fully vernalized - and thermal time ! accumulation in phase 1 will be unaffected ! refers to gddtsoi & gddplant, defined in the accumulation routines (slevis) - ! reset vf, cumvd, and hdidx to 0 at planting of crop (slevis) - - if (t_ref2m_max(p) > tfrz) then - if (t_ref2m_min(p) <= tfrz+15._r8) then - vd1 = 1.4_r8 - 0.0778_r8 * tcrown - vd2 = 0.5_r8 + 13.44_r8 / ((t_ref2m_max(p)-t_ref2m_min(p)+3._r8)**2) * tcrown - vd = max(0._r8, min(1._r8, vd1, vd2)) - cumvd(p) = cumvd(p) + vd - end if + ! reset vf and cumvd to 0 at planting of crop (slevis) + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + ! Modifications based on Yaqiong Lu et al., 2017 in Geosci. Model Dev. + ! A generalized vernalization response function for winter wheat (Streck + ! et al.,2003) was used here + ! Streck, N.A., Weiss, A., Baenziger, P.S., 2003. A generalized + ! vernalization response & + ! function for winter wheat. Agron. J. 95 + + dt_hr = dtime_mod/3600.0_r8 ! dt_hr is the time step in hour + + alpha = log(2._r8)/log((vtmax - vtmin)/(vtopt - vtmin)) + + tc = t_ref2m(p) - tfrz + if(tc >= vtmin .and. tc <= vtmax) then + ! Fixes minor typo in Yaqiong Lu et al., 2017 code + cumvd(p) = cumvd(p) + ((2._r8 * ((tc - vtmin)**alpha)*(vtopt - vtmin)**alpha) & + - (tc - vtmin)**(2._r8 * alpha))/(vtopt - vtmin)**(2._r8 * alpha) * (dt_hr/24._r8) + end if - if (cumvd(p) < 10._r8 .and. t_ref2m_max(p) > tfrz+30._r8) then - cumvd(p) = cumvd(p) - 0.5_r8 * (t_ref2m_max(p) - tfrz - 30._r8) - end if - cumvd(p) = max(0._r8, cumvd(p)) ! must be > 0 + vf(p) = (cumvd(p)**5._r8)/(22.5_r8**5._r8 + cumvd(p)**5._r8) + + end associate + + end subroutine vernalization + + !----------------------------------------------------------------------- + ! Added based on Yaqiong Lu et al., 2017 in Geosci. Model Dev. + subroutine coldtolerance(p, cnstate_vars, crop_vars) + ! + ! !DESCRIPTION: + ! + ! * * * only call for winter temperate cereal * * * + ! + !the subroutine calculates the lethal temperature at 50% of crop alive, + !survival rate, winter degree days + ! !USES: + use elm_time_manager, only: get_step_size,get_curr_date + ! + ! !ARGUMENTS: + implicit none + integer , intent(in) :: p ! PATCH index running over + type(cnstate_type) , intent(inout) :: cnstate_vars + type(crop_type) , intent(inout) :: crop_vars + ! + ! LOCAL VARAIBLES: + logical :: end_cd ! temporary for is_end_curr_day() value + real(r8) tc ! t_ref2m in degree C + real(r8) prevleafc ! previous step leafc + real(r8) tempfsurv ! averaged survival rate + integer c,g ! indices + integer kyr ! current year + integer kmo ! month of year (1, ..., 12) + integer kda ! day of month (1, ..., 31) + integer mcsec ! seconds of day (0, ..., seconds/day) + real(r8) :: dt + real(r8) dt_hr ! convert dt from sec to hour + real(r8), parameter :: Hparam=0.0093 ! based on Lu et al., 2017 which was based on Bergjord et al. (2008) + real(r8), parameter :: Dparam=2.7e-5 ! based on Lu et al., 2017 which was based on Bergjord et al. (2008) + real(r8), parameter :: Sparam=1.9 + real(r8), parameter :: Rparam=0.54 + real(r8), parameter :: T_S_max=12.5 + real(r8), parameter :: lt50max=-23 + + !the calculation of frost tolerance is based on Bergjord et + !al.,(2008), Europ. J. Agronomy + !the calculation of survival rate and WDD is based on Vico et al.,(2014),Agri + !and Forest Metero. + + !------------------------------------------------------------------------ + + associate( & + ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type + rateh => crop_vars%rateh_patch , & ! Output: [integer (:) ] increase of tolerance caused by cold hardening index + rated => crop_vars%rated_patch , & ! Output: [integer (:) ] loss of tolerance caused by dehardening + rates => crop_vars%rates_patch , & ! Output: [integer (:) ] loss of tolerance caused by low temperature + rater => crop_vars%rater_patch , & ! Output: [integer (:) ] loss of tolerance caused by respiration under snow + lt50 => crop_vars%lt50_patch , & ! Output: [integer (:) ] the lethal temperature at which 50% of the individuals are damaged + fsurv => crop_vars%fsurv_patch , & ! Output: [integer (:) ] winter wheat survival rate + accfsurv => crop_vars%accfsurv_patch , & ! Output: [integer (:) ] accumulated winter wheat survival rate + countfsurv => crop_vars%countfsurv_patch , & ! Output: [integer (:) ] count of accumulated winter wheat survival rate + wdd => crop_vars%wdd_patch , & ! Output: [integer (:) ] winter wheat weighted cumulated degree days + tcrown => crop_vars%tcrown_patch , & ! Output: [integer (:) ] crown temperature + vf => crop_vars%vf_patch , & ! Output: [real(r8) (:) ] vernalization factor + leafc_to_litter => veg_cf%leafc_to_litter , & ! Input: [real(r8) (:) ] leaf C litterfall (gC/m2/s) + leafn_to_litter => veg_nf%leafn_to_litter , & ! Output: [real(r8) (:) ] leaf N litterfall (gN/m2/s) + leafc => veg_cs%leafc , & ! Input: [real(r8) (:) ] (gC/m2) leaf C + leafn => veg_ns%leafn , & ! Input: [real(r8) (:) ] (gN/m2) leaf N + lflitcn => veg_vp%lflitcn , & ! Input: [real(r8) (:) ] leaf litter C:N (gC/gN) + t_ref2m => veg_es%t_ref2m , & ! Input: [real(r8) (:) ] 2m air temperature (K) + snow_depth => col_ws%snow_depth & ! Input: [real(r8) (:) ] snow height (m) + ) + + dt = dtime_mod + dt_hr = dt/3600.0_r8 ! dt_hr is the time step in hour + + c = veg_pp%column(p) + + ! for all equations - temperatures must be in degrees (C) + ! calculate temperature of crown of crop (e.g., 3 cm soil temperature) + ! snow depth in centimeters + + if (t_ref2m(p) < tfrz) then !slevis: t_ref2m inst of td=daily avg (K) + tcrown = 2._r8 + (t_ref2m(p) - tfrz) * (0.4_r8 + 0.0018_r8 * & + (min(snow_depth(c)*100._r8, 15._r8) - 15._r8)**2) + else !slevis: snow_depth inst of adsnod=daily average (m) + tcrown = t_ref2m(p) - tfrz + end if + + ! frost tolerance and survival rate calculation - vf(p) = 1._r8 - p1v * (50._r8 - cumvd(p)) - vf(p) = max(0._r8, min(vf(p), 1._r8)) ! must be between 0 - 1 + if(tcrown(p) < 10._r8) then + ! Fixes minor typo in Yaqiong Lu et al., 2017 code + rateh(p) = Hparam * (10._r8 - max(tcrown(p), 0._r8)) * (lt50(p) - lt50max) end if - ! calculate cold hardening of plant - ! determines for winter cereal varieties whether the plant has completed - ! a period of cold hardening to protect it from freezing temperatures. If - ! not, then exposure could result in death or killing of plants. + if((tcrown(p) >= -4._r8 .and. vf(p) == 1._r8) .or. (tcrown(p) >= 10._r8 .and. vf(p) < 1._r8)) then + rated(p) = Dparam * (-0.6_r8 + 0.142_r8 * lt50max - lt50(p)) * (tcrown(p) + 4._r8)**3._r8 + end if - ! there are two distinct phases of hardening + ! Fixes minor typo in Yaqiong Lu et al., 2017 code + rater(p) = Rparam * (exp(0.84 + 0.051 * tcrown(p)) - 2._r8)/1.85_r8 * (min(snow_depth(c)*100._r8, 12.5_r8))/12.5 + rates(p) = (lt50(p) - tcrown(p))/exp(-Sparam * (lt50(p) - tcrown(p)) - 3.74_r8) + lt50(p) = lt50(p) + (rated(p) + rates(p) + rater(p) - rateh(p)) * (dt_hr/24._r8) - if (t_ref2m_min(p) <= tfrz-3._r8 .or. hdidx(p) /= 0._r8) then - if (hdidx(p) >= hti) then ! done with phase 1 - hdidx(p) = hdidx(p) + 0.083_r8 - hdidx(p) = min(hdidx(p), hti*2._r8) - end if + fsurv(p) = 2._r8**(-(abs(tcrown(p))/abs(lt50(p)))**4._r8) + wdd(p) = wdd(p) + (max(tbase - tcrown(p), 0._r8) * (1._r8 - fsurv(p))) * (dt_hr/24._r8) - if (t_ref2m_max(p) >= tbase + tfrz + 10._r8) then - hdidx(p) = hdidx(p) - 0.02_r8 * (t_ref2m_max(p)-tbase-tfrz-10._r8) - if (hdidx(p) > hti) hdidx(p) = hdidx(p) - 0.02_r8 * (t_ref2m_max(p)-tbase-tfrz-10._r8) - hdidx(p) = max(0._r8, hdidx(p)) - end if + if(wdd(p) > 0._r8) then + accfsurv(p) = accfsurv(p) + fsurv(p) + countfsurv(p) = countfsurv(p) + 1._r8 + end if - else if (tcrown >= tbase-1._r8) then - if (tcrown <= tbase+8._r8) then - hdidx(p) = hdidx(p) + 0.1_r8 - (tcrown-tbase+3.5_r8)**2 / 506._r8 - if (hdidx(p) >= hti .and. tcrown <= tbase + 0._r8) then - hdidx(p) = hdidx(p) + 0.083_r8 - hdidx(p) = min(hdidx(p), hti*2._r8) - end if - end if + call get_curr_date(kyr, kmo, kda, mcsec) - if (t_ref2m_max(p) >= tbase + tfrz + 10._r8) then - hdidx(p) = hdidx(p) - 0.02_r8 * (t_ref2m_max(p)-tbase-tfrz-10._r8) - if (hdidx(p) > hti) hdidx(p) = hdidx(p) - 0.02_r8 * (t_ref2m_max(p)-tbase-tfrz-10._r8) - hdidx(p) = max(0._r8, hdidx(p)) - end if + end_cd = (mcsec == 0) + if(end_cd .and. wdd(p) > 0._r8 .and. vf(p) < 0.9_r8 .and. leafc(p) > 10._r8 ) then + leafc_to_litter(p) = leafc_to_litter(p) + 5._r8 * (1._r8 - fsurv(p))/dt + leafn_to_litter(p) = leafn_to_litter(p) + (5._r8 * (1._r8 - fsurv(p)))/(dt * lflitcn(ivt(p))) end if - ! calculate what the cereal killing temperature - ! there is a linear inverse relationship between - ! hardening of the plant and the killing temperature or - ! threshold that the plant can withstand - ! when plant is fully-hardened (hdidx = 2), the killing threshold is -18 C - - ! will have to develop some type of relationship that reduces LAI and - ! biomass pools in response to cold damaged crop - - if (t_ref2m_min(p) <= tfrz - 6._r8) then - tkil = (tbase - 6._r8) - 6._r8 * hdidx(p) - if (tkil >= tcrown) then - if ((0.95_r8 - 0.02_r8 * (tcrown - tkil)**2) >= 0.02_r8) then - write (iulog,*) 'crop damaged by cold temperatures at p,c =', p,c - else if (tlai(p) > 0._r8) then ! slevis: kill if past phase1 - gddmaturity(p) = 0._r8 ! by forcing through - huigrain(p) = 0._r8 ! harvest - write (iulog,*) '95% of crop killed by cold temperatures at p,c =', p,c - end if - end if + if(end_cd .and. wdd(p) > 1.0_r8 .and. vf(p) > 0.9_r8 ) then + tempfsurv = accfsurv(p)/countfsurv(p) + prevleafc = leafc(p) + leafc_to_litter(p) = leafc_to_litter(p) + prevleafc * (1._r8 - tempfsurv)/dt + leafn_to_litter(p) = leafn_to_litter(p) + (prevleafc * (1._r8 - tempfsurv))/(dt * lflitcn(ivt(p))) + accfsurv(p) = 1._r8 + countfsurv(p) = 1._r8 + wdd(p) = 0._r8 end if end associate - end subroutine vernalization + end subroutine coldtolerance !----------------------------------------------------------------------- subroutine CropPlantDate (num_soilp, filter_soilp, num_pcropp, filter_pcropp, & diff --git a/components/elm/src/data_types/VegetationDataType.F90 b/components/elm/src/data_types/VegetationDataType.F90 index 3bc7021f1509..a1b90d925b78 100644 --- a/components/elm/src/data_types/VegetationDataType.F90 +++ b/components/elm/src/data_types/VegetationDataType.F90 @@ -1567,6 +1567,7 @@ subroutine update_acc_vars_veg_es (this, bounds) use shr_const_mod , only : SHR_CONST_CDAY, SHR_CONST_TKFRZ use elm_time_manager , only : get_step_size, get_nstep, is_end_curr_day, get_curr_date use accumulMod , only : update_accum_field, extract_accum_field, accumResetVal + use pftvarcon , only: nwcereal, nwcerealirrig ! ! !ARGUMENTS: class(vegetation_energy_state) :: this @@ -1586,6 +1587,10 @@ subroutine update_acc_vars_veg_es (this, bounds) real(r8), pointer :: rbufslp(:) ! temporary single level - pft level !--------------------------------------------------------------------- + associate( & + ivt => veg_pp%itype & ! Input: [integer (:) ] pft vegetation type + ) + begp = bounds%begp; endp = bounds%endp dtime = get_step_size() @@ -1718,13 +1723,29 @@ subroutine update_acc_vars_veg_es (this, bounds) do p = begp,endp g = veg_pp%gridcell(p) - if (month==1 .and. day==1 .and. secs==int(dtime)) then - rbufslp(p) = accumResetVal ! reset gdd - else if (( month > 3 .and. month < 10 .and. grc_pp%latdeg(g) >= 0._r8) .or. & - ((month > 9 .or. month < 4) .and. grc_pp%latdeg(g) < 0._r8) ) then - rbufslp(p) = max(0._r8, min(26._r8, this%t_ref2m(p)-SHR_CONST_TKFRZ)) * dtime/SHR_CONST_CDAY + + ! Added based on Yaqiong Lu et al., 2017 in Geosci. Model Dev. + ! Accumulate GDD0 for winter wheat from September to June in NH + ! April to December in SH, the accumulated period may be less than the actual + ! winter wheat growing season + if(ivt(p) == nwcereal .or. ivt(p) == nwcerealirrig) then + if (month==9 .and. day==1 .and. secs==int(dtime)) then + rbufslp(p) = accumResetVal ! reset gdd + else if (( month > 8 .or. month < 7 .and. grc_pp%latdeg(g) >= 0._r8) .or. & + ((month > 3 .and. month < 10) .and. grc_pp%latdeg(g) < 0._r8)) then + rbufslp(p) = max(0._r8, min(26._r8, this%t_ref2m(p)-SHR_CONST_TKFRZ)) * dtime/SHR_CONST_CDAY + else + rbufslp(p) = 0._r8 ! keeps gdd unchanged at other times (eg,through Dec in NH) + end if else - rbufslp(p) = 0._r8 ! keeps gdd unchanged at other times (eg, through Dec in NH) + if (month==1 .and. day==1 .and. secs==int(dtime)) then + rbufslp(p) = accumResetVal ! reset gdd + else if (( month > 3 .and. month < 10 .and. grc_pp%latdeg(g) >= 0._r8) .or. & + ((month > 9 .or. month < 4) .and. grc_pp%latdeg(g) < 0._r8) ) then + rbufslp(p) = max(0._r8, min(26._r8, this%t_ref2m(p)-SHR_CONST_TKFRZ)) * dtime/SHR_CONST_CDAY + else + rbufslp(p) = 0._r8 ! keeps gdd unchanged at other times (eg, through Dec in NH) + end if end if end do call update_accum_field ('GDD0', rbufslp, nstep) @@ -1766,6 +1787,7 @@ subroutine update_acc_vars_veg_es (this, bounds) end if deallocate(rbufslp) + end associate end subroutine update_acc_vars_veg_es