From e22260b87f91996f2a210d0cffd510e1830c1cd1 Mon Sep 17 00:00:00 2001 From: Sam Rabin Date: Tue, 19 Sep 2023 13:53:05 -0600 Subject: [PATCH] Prevent crops from being sown twice in one sowing window. --- src/biogeochem/CNPhenologyMod.F90 | 21 ++++++++++++++++----- src/biogeochem/CropType.F90 | 6 ++++++ 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/biogeochem/CNPhenologyMod.F90 b/src/biogeochem/CNPhenologyMod.F90 index 608b149dee..5ada26f377 100644 --- a/src/biogeochem/CNPhenologyMod.F90 +++ b/src/biogeochem/CNPhenologyMod.F90 @@ -1726,7 +1726,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & use clm_time_manager , only : get_prev_calday, get_curr_days_per_year, is_beg_curr_year use clm_time_manager , only : get_average_days_per_year use clm_time_manager , only : get_prev_date - use clm_time_manager , only : is_doy_in_interval + use clm_time_manager , only : is_doy_in_interval, is_end_curr_day use pftconMod , only : ntmp_corn, nswheat, nwwheat, ntmp_soybean use pftconMod , only : nirrig_tmp_corn, nirrig_swheat, nirrig_wwheat, nirrig_tmp_soybean use pftconMod , only : ntrp_corn, nsugarcane, ntrp_soybean, ncotton, nrice @@ -1962,6 +1962,9 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & sowing_window_enddate = maxplantjday(ivt(p),h) end if is_in_sowing_window = is_doy_in_interval(sowing_window_startdate, sowing_window_enddate, jday) + if (crop_inst%sown_in_this_window .and. .not. is_in_sowing_window) then + crop_inst%sown_in_this_window = .false. + end if is_end_sowing_window = jday == sowing_window_enddate ! ! Save these diagnostic variables only on the first day of the window to ensure that windows spanning the new year aren't double-counted. @@ -2199,8 +2202,8 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & ! If generate_crop_gdds and this patch has prescribed sowing inputs else if (generate_crop_gdds .and. crop_inst%rx_swindow_starts_thisyr_patch(p,1) .gt. 0) then - if (next_rx_swindow_start(p) >= 0) then - ! Harvest the day before the start of the next sowing window this year. + if (next_rx_swindow_start(p) >= 0 .and. is_end_curr_day()) then + ! Harvest the timestep before the start of the next sowing window this year. do_harvest = jday == next_rx_swindow_start(p) - 1 ! ... unless that will lead to growing season length 365 (or 366, @@ -2229,7 +2232,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & ! In order to avoid this, you'd have to read this year's AND next year's prescribed ! sowing window start dates. if (crop_inst%rx_swindow_starts_thisyr_patch(p,1) == 1) then - do_harvest = jday == dayspyr + do_harvest = jday == dayspyr .and. is_end_curr_day() end if if (do_harvest) then @@ -2257,7 +2260,8 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & if (use_cropcal_rx_swindows) then will_plant_prescribed_tomorrow = (jday == next_rx_swindow_start(p) - 1) .or. & (crop_inst%sdates_thisyr_patch(p,1) == 1 .and. & - jday == dayspyr) + jday == dayspyr .and. + is_end_curr_day()) else will_plant_prescribed_tomorrow = .false. end if @@ -2323,6 +2327,12 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & croplive(p) = .false. ! no re-entry in greater if-block cphase(p) = cphase_harvest + + ! Avoid situation where a crop could be stuck with sown_in_this_window true when moving from one sowing window to another with no days in between. + if (harvest_reason == HARVEST_REASON_SOWTOMORROW .or. harvest_reason == HARVEST_REASON_IDOPTOMORROW) then + crop_inst%sown_in_this_window = .false. + end if + if (tlai(p) > 0._r8) then ! plant had emerged before harvest offset_flag(p) = 1._r8 offset_counter(p) = dt @@ -2600,6 +2610,7 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & ! impose limit on growing season length needed ! for crop maturity - for cold weather constraints croplive(p) = .true. + crop_inst%sown_in_this_window = .true. idop(p) = jday iyop(p) = kyr harvdate(p) = NOT_Harvested diff --git a/src/biogeochem/CropType.F90 b/src/biogeochem/CropType.F90 index 2239d3e9bd..ae56bd8ad6 100644 --- a/src/biogeochem/CropType.F90 +++ b/src/biogeochem/CropType.F90 @@ -47,6 +47,7 @@ module CropType character(len=20) :: baset_mapping real(r8) :: baset_latvary_intercept real(r8) :: baset_latvary_slope + logical , pointer :: sown_in_this_window (:) ! patch flag. True if the crop has already been sown during the current sowing window. False otherwise or if not in a sowing window. integer , pointer :: next_rx_swindow_start_patch (:) ! start of prescribed sowing window for the next growing season this year integer , pointer :: next_rx_swindow_end_patch (:) ! end of prescribed sowing window for the next growing season this year integer , pointer :: rx_swindow_starts_thisyr_patch(:,:) ! all prescribed sowing window start dates for this patch this year (day of year) [patch, mxsowings] @@ -231,6 +232,7 @@ subroutine InitAllocate(this, bounds) allocate(this%cphase_patch (begp:endp)) ; this%cphase_patch (:) = cphase_not_planted allocate(this%sowing_reason_patch (begp:endp)) ; this%sowing_reason_patch (:) = -1 allocate(this%latbaset_patch (begp:endp)) ; this%latbaset_patch (:) = spval + allocate(this%sown_in_this_window(begp:endp)) ; this%sown_in_this_window(:) = .false. allocate(this%next_rx_swindow_start_patch(begp:endp)) ; this%next_rx_swindow_start_patch(:) = -1 allocate(this%next_rx_swindow_end_patch (begp:endp)) ; this%next_rx_swindow_end_patch (:) = -1 allocate(this%rx_swindow_starts_thisyr_patch(begp:endp,1:mxsowings)); this%rx_swindow_starts_thisyr_patch(:,:) = -1 @@ -627,6 +629,10 @@ subroutine Restart(this, bounds, ncid, cnveg_state_inst, flag) dim1name='pft', long_name='sowing reason for this patch', & units='none', & interpinic_flag='interp', readvar=readvar, data=this%sowing_reason_patch) + call restartvar(ncid=ncid, flag=flag, varname='sown_in_this_window', xtype=ncd_log, & + dim1name='pft', & + long_name='whether this patch was sown already during the current sowing window', & + interpinic_flag='interp', readvar=readvar, data=this%sown_in_this_window) ! Read or write variable(s) with mxsowings dimension ! BACKWARDS_COMPATIBILITY(wjs/ssr, 2022-02-02) See note in CallRestartvarDimOK()