diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index c12c54645f..65cef8b77e 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -4010,21 +4010,28 @@ sub setup_logic_cropcal_streams { # Option checks my $generate_crop_gdds = $nl->get_value('generate_crop_gdds') ; my $use_mxmat = $nl->get_value('use_mxmat') ; - my $sdate_file = $nl->get_value('stream_fldFileName_sdate') ; + my $swindow_start_file = $nl->get_value('stream_fldFileName_swindow_start') ; + my $swindow_end_file = $nl->get_value('stream_fldFileName_swindow_end') ; my $gdd_file = $nl->get_value('stream_fldFileName_cultivar_gdds') ; my $mesh_file = $nl->get_value('stream_meshfile_cropcal') ; + if ( ($swindow_start_file eq '' and $swindow_start_file ne '') or ($swindow_start_file ne '' and $swindow_start_file eq '') ) { + $log->fatal_error("When specifying sowing window dates, you must provide both swindow_start_file and swindow_end_file. To specify exact sowing dates, use the same file." ); + } if ( $generate_crop_gdds eq '.true.' ) { if ( $use_mxmat eq '.true.' ) { $log->fatal_error("If generate_crop_gdds is true, you must also set use_mxmat to false" ); } - if ( $sdate_file eq '' ) { - $log->fatal_error("If generate_crop_gdds is true, you must specify stream_fldFileName_sdate") + if ( $swindow_start_file eq '' or $swindow_end_file eq '' ) { + $log->fatal_error("If generate_crop_gdds is true, you must specify stream_fldFileName_swindow_start and stream_fldFileName_swindow_end") + } + if ( $swindow_start_file ne $swindow_end_file ) { + $log->fatal_error("If generate_crop_gdds is true, you must specify exact sowing dates by setting stream_fldFileName_swindow_start and stream_fldFileName_swindow_end to the same file") } if ( $gdd_file ne '' ) { $log->fatal_error("If generate_crop_gdds is true, do not specify stream_fldFileName_cultivar_gdds") } } - if ( $mesh_file eq '' and ( $sdate_file ne '' or $gdd_file ne '' ) ) { + if ( $mesh_file eq '' and ( $swindow_start_file ne '' or $gdd_file ne '' ) ) { $log->fatal_error("If prescribing crop sowing dates and/or maturity requirements, you must specify stream_meshfile_cropcal") } } diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index a3c34dc337..bca1beca6a 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1893,9 +1893,19 @@ Last year to loop over for crop calendar data Simulation year that aligns with stream_year_first_cropcal value - +By default, a value in stream_fldFileName_swindow_start or _end outside the range [1, 365] (or 366 in leap years) will cause the run to fail. Set this to .true. to instead fall back on the paramfile sowing windows. + + + +Filename of input stream data for date (day of year) of start of sowing window. Cells with the same sowing window start and end date are always planted on that date, regardless of climatic conditions/history. + + + -Filename of input stream data for sowing dates +Filename of input stream data for date (day of year) of end of sowing window. Cells with the same sowing window start and end date are always planted on that date, regardless of climatic conditions/history. maxval(rx_ends)) then + w = 1 + start_w = rx_starts(w) + end_w = rx_ends(w) + return + end if + + ! Otherwise, use the first prescribed sowing window we find whose end is >= today. This works only if sowing windows that span the new year are located at index w = 1. + do w = 1, mxsowings_in + ! If nothing prescribed at this w, stop looking and exit loop. Will trigger "No sowing window found" error, which we do not move here because it's possible that no start or end date is < 1. + if (min(rx_starts(w), rx_ends(w)) < 1) then + exit + end if + + if (jday <= rx_ends(w)) then + start_w = rx_starts(w) + end_w = rx_ends(w) + exit + end if + end do + + ! Ensure that a window was found. + ! SSR 2023-10-17: This shouldn't currently be reachable, but its being here casts the widest possible net in case code changes in future. + if (start_w < 1 .or. end_w < 1) then + call endrun(msg="get_swindow(): No sowing window found") + end if + + end subroutine get_swindow + + + !----------------------------------------------------------------------- + function was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, sown_in_this_window) + ! !DESCRIPTION: + ! Determine whether the crop was sown in the current sowing window. Although sown_in_this_window is set to false in last timestep of sowing window at the end of CropPhenology(), these extra checks may be necessary if sowing windows change. + ! + ! !USES: + use clm_time_manager , only : is_doy_in_interval + ! !ARGUMENTS: + integer, intent(in) :: sowing_window_startdate, sowing_window_enddate, jday, idop + logical, intent(in) :: sown_in_this_window + ! !LOCAL VARIABLES + logical :: is_in_sowing_window, idop_in_sowing_window + ! !RESULT + logical :: was_sown_in_this_window + + was_sown_in_this_window = sown_in_this_window + + ! If not in a sowing window, sown_in_this_window must be false. + is_in_sowing_window = is_doy_in_interval(sowing_window_startdate, sowing_window_enddate, jday) + if (.not. is_in_sowing_window) then + was_sown_in_this_window = .false. + return + end if + + ! If we're in a sowing window but the day of planting isn't in the active sowing window, we must be in a different sowing window. + idop_in_sowing_window = is_doy_in_interval(sowing_window_startdate, sowing_window_enddate, idop) + if (is_in_sowing_window .and. .not. idop_in_sowing_window) then + was_sown_in_this_window = .false. + return + end if + + ! Sometimes we're in an active sowing window, and the patch was sown between the start and end dates of the window, but not *the currently active* window. Note that windows with start==end are not checked here; we always trust the input value of sown_in_this_window in such cases. + if (sowing_window_startdate < sowing_window_enddate .and. idop > jday) then + was_sown_in_this_window = .false. + else if (sowing_window_startdate > sowing_window_enddate) then + if (jday <= sowing_window_enddate .and. idop <= sowing_window_enddate .and. idop > jday) then + was_sown_in_this_window = .false. + else if (jday >= sowing_window_startdate .and. (idop > jday .or. idop <= sowing_window_enddate)) then + was_sown_in_this_window = .false. + end if + end if + + end function was_sown_in_this_window + + !----------------------------------------------------------------------- subroutine CropPhenology(num_pcropp, filter_pcropp , & waterdiagnosticbulk_inst, temperature_inst, crop_inst, canopystate_inst, cnveg_state_inst , & @@ -1726,7 +1840,8 @@ 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 clm_time_manager , only : get_doy_tomorrow 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 @@ -1738,7 +1853,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & use clm_varctl , only : use_fertilizer use clm_varctl , only : use_c13, use_c14 use clm_varcon , only : c13ratio, c14ratio - use clm_varctl , only : use_cropcal_rx_sdates, use_cropcal_rx_cultivar_gdds, use_cropcal_streams + use clm_varctl , only : use_cropcal_rx_swindows, use_cropcal_rx_cultivar_gdds, use_cropcal_streams ! ! !ARGUMENTS: integer , intent(in) :: num_pcropp ! number of prog crop patches in filter @@ -1763,30 +1878,34 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & integer h ! hemisphere indices integer s ! growing season indices integer k ! grain pool indices + integer w ! sowing window index integer idpp ! number of days past planting integer mxmat ! maximum growing season length 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) + integer sowing_window_startdate ! date (day of year) of first day of sowing window + integer sowing_window_enddate ! date (day of year) of last day of sowing window real(r8) harvest_reason real(r8) dayspyr ! days per year in this year real(r8) avg_dayspyr ! average number of days per year real(r8) crmcorn ! comparitive relative maturity for corn real(r8) ndays_on ! number of days to fertilize + logical has_rx_sowing_date ! does the crop have a single sowing date instead of a window? logical is_in_sowing_window ! is the crop in its sowing window? logical is_end_sowing_window ! is it the last day of the crop's sowing window? logical sowing_gdd_requirement_met ! has the gridcell historically been warm enough to support the crop? logical do_plant_normal ! are the normal planting rules defined and satisfied? logical do_plant_lastchance ! if not the above, what about relaxed rules for the last day of the planting window? logical do_plant_prescribed ! is today the prescribed sowing date? + logical do_plant_prescribed_tomorrow ! is tomorrow the prescribed sowing date? logical do_plant ! are we planting in this time step for any reason? logical did_plant ! did we plant the crop in this time step? logical allow_unprescribed_planting ! should crop be allowed to be planted according to sowing window rules? logical do_harvest ! Are harvest conditions satisfied? logical fake_harvest ! Dealing with incorrect Dec. 31 planting logical did_plant_prescribed_today ! Was the crop sown today? - logical will_plant_prescribed_tomorrow ! Is tomorrow a prescribed sowing day? logical vernalization_forces_harvest ! Was the crop killed by freezing during vernalization? !------------------------------------------------------------------------ @@ -1811,11 +1930,10 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & fertnitro => crop_inst%fertnitro_patch , & ! Input: [real(r8) (:) ] fertilizer nitrogen hui => crop_inst%hui_patch , & ! Input: [real(r8) (:) ] crop patch heat unit index (growing degree-days); set to 0 at sowing and accumulated until harvest - leafout => crop_inst%gddtsoi_patch , & ! Input: [real(r8) (:) ] gdd from top soil layer temperature + leafout => crop_inst%gddtsoi_patch , & ! Input: [real(r8) (:) ] gdd from top soil layer temperature harvdate => crop_inst%harvdate_patch , & ! Output: [integer (:) ] harvest date - croplive => crop_inst%croplive_patch , & ! Output: [logical (:) ] Flag, true if planted, not harvested + croplive => crop_inst%croplive_patch , & ! Output: [logical (:) ] Flag, true if planted, not harvested vf => crop_inst%vf_patch , & ! Output: [real(r8) (:) ] vernalization factor - next_rx_sdate => crop_inst%next_rx_sdate_patch , & ! Inout: [integer (:) ] prescribed sowing date of next growing season this year sowing_count => crop_inst%sowing_count , & ! Inout: [integer (:) ] number of sowing events this year for this patch harvest_count => crop_inst%harvest_count , & ! Inout: [integer (:) ] number of harvest events this year for this patch peaklai => cnveg_state_inst%peaklai_patch , & ! Output: [integer (:) ] 1: max allowed lai; 0: not at max @@ -1889,6 +2007,8 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & harvest_count(p) = 0 do s = 1, mxsowings crop_inst%sdates_thisyr_patch(p,s) = -1._r8 + crop_inst%swindow_starts_thisyr_patch(p,s) = -1._r8 + crop_inst%swindow_ends_thisyr_patch (p,s) = -1._r8 crop_inst%sowing_reason_thisyr_patch(p,s) = -1._r8 end do do s = 1, mxharvests @@ -1913,17 +2033,26 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & cnveg_nitrogenflux_inst%repr_grainn_to_food_thisyr_patch(p,k) = 0._r8 cnveg_nitrogenflux_inst%repr_grainn_to_seed_thisyr_patch(p,k) = 0._r8 end do - next_rx_sdate(p) = crop_inst%rx_sdates_thisyr_patch(p,1) end if - ! Get next sowing date - if (sowing_count(p) < mxsowings) then - next_rx_sdate(p) = crop_inst%rx_sdates_thisyr_patch(p,sowing_count(p)+1) - end if + ! Get dates of current or next sowing window. + call get_swindow(jday, crop_inst%rx_swindow_starts_thisyr_patch(p,:), crop_inst%rx_swindow_ends_thisyr_patch(p,:), minplantjday(ivt(p),h), maxplantjday(ivt(p),h), w, sowing_window_startdate, sowing_window_enddate) + + ! Are we currently in a sowing window? + ! This is outside the croplive check so that the "harvest if planting conditions were met today" conditional works. + is_in_sowing_window = is_doy_in_interval(sowing_window_startdate, sowing_window_enddate, jday) + crop_inst%sown_in_this_window(p) = was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop(p), crop_inst%sown_in_this_window(p)) + is_end_sowing_window = jday == sowing_window_enddate + + ! We only want to plant on a specific day if the prescribed sowing window starts AND ends on the same day. Also make sure we haven't planted yet today. + has_rx_sowing_date = sowing_window_startdate == sowing_window_enddate + do_plant_prescribed = has_rx_sowing_date .and. & + sowing_window_startdate == jday .and. & + .not. crop_inst%sown_in_this_window(p) + do_plant_prescribed_tomorrow = & + has_rx_sowing_date .and. & + sowing_window_startdate == get_doy_tomorrow(jday) - ! CropType->InitAllocate() initializes next_rx_sdate to -1. It's only changed from that, by cropCalStreamMod->cropcal_interp(), when use_cropcal_rx_sdates is true. So if not using prescribed sowing dates, this boolean will always be false, because jday can never be -1. - do_plant_prescribed = next_rx_sdate(p) == jday .and. sowing_count(p) < mxsowings - ! BACKWARDS_COMPATIBILITY(wjs/ssr, 2022-02-18) ! When resuming from a run with old code, may need to manually set these. ! Will be needed until we can rely on all restart files have been generated @@ -1936,14 +2065,15 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & crop_inst%sdates_thisyr_patch(p,1) = real(idop(p), r8) end if - ! This is outside the croplive check so that the "harvest if planting conditions were met today" conditional works. - is_in_sowing_window = is_doy_in_interval(minplantjday(ivt(p),h), maxplantjday(ivt(p),h), jday) - is_end_sowing_window = jday == maxplantjday(ivt(p),h) + ! Save these diagnostic variables only on the last day of the window to ensure that windows spanning the new year aren't double-counted. Doing this on the last day ensures that outputs are ordered as inputs should be. + if (jday == sowing_window_enddate) then + crop_inst%swindow_starts_thisyr_patch(p,w) = sowing_window_startdate + crop_inst%swindow_ends_thisyr_patch (p,w) = sowing_window_enddate + end if ! ! Only allow sowing according to normal "window" rules if not using prescribed - ! sowing dates at all, or if this cell had no values in the prescribed sowing - ! date file. - allow_unprescribed_planting = (.not. use_cropcal_rx_sdates) .or. crop_inst%rx_sdates_thisyr_patch(p,1)<0 + ! sowing dates. + allow_unprescribed_planting = .not. has_rx_sowing_date if (sowing_count(p) == mxsowings) then do_plant_normal = .false. do_plant_lastchance = .false. @@ -1980,6 +2110,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & gdd820(p) /= spval end if do_plant = do_plant_prescribed .or. do_plant_normal .or. do_plant_lastchance + do_plant = do_plant .and. .not. crop_inst%sown_in_this_window(p) did_plant = .false. ! Once outputs can handle >1 planting per year, remove 2nd condition. @@ -2146,7 +2277,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & do_harvest = .false. fake_harvest = .false. did_plant_prescribed_today = .false. - if (use_cropcal_rx_sdates .and. sowing_count(p) > 0) then + if (use_cropcal_rx_swindows .and. sowing_count(p) > 0) then did_plant_prescribed_today = crop_inst%sdates_thisyr_patch(p,sowing_count(p)) == real(jday, r8) end if @@ -2167,45 +2298,23 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & harvest_reason = HARVEST_REASON_SOWTODAY ! If generate_crop_gdds and this patch has prescribed sowing inputs - else if (generate_crop_gdds .and. crop_inst%rx_sdates_thisyr_patch(p,1) .gt. 0) then - if (next_rx_sdate(p) >= 0) then - ! Harvest the day before the next sowing date this year. - do_harvest = jday == next_rx_sdate(p) - 1 - - ! ... unless that will lead to growing season length 365 (or 366, - ! if last year was a leap year). This would result in idop==jday, - ! which would invoke the "manually setting sowing_count and - ! sdates_thisyr" code. This would lead to crops never getting - ! harvested. Instead, always harvest the day before idop. - if ((.not. do_harvest) .and. & - (idop(p) > 1 .and. jday == idop(p) - 1) .or. & - (idop(p) == 1 .and. jday == dayspyr)) then - do_harvest = .true. - if (do_harvest) then - harvest_reason = HARVEST_REASON_IDOPTOMORROW - end if - else if (do_harvest) then - harvest_reason = HARVEST_REASON_SOWTOMORROW - end if - - else - ! If this patch has already had all its plantings for the year, don't harvest - ! until some time next year. - do_harvest = .false. - - ! ... unless first sowing next year happens Jan. 1. - ! WARNING: This implementation assumes that sowing dates don't change over time! - ! In order to avoid this, you'd have to read this year's AND next year's prescribed - ! sowing dates. - if (crop_inst%rx_sdates_thisyr_patch(p,1) == 1) then - do_harvest = jday == dayspyr - end if - - if (do_harvest) then - harvest_reason = HARVEST_REASON_SOWTOMORROW - end if - - endif + else if (generate_crop_gdds .and. crop_inst%rx_swindow_starts_thisyr_patch(p,1) .gt. 0) then + ! Harvest the day before the next prescribed sowing. + do_harvest = do_plant_prescribed_tomorrow + + ! ... unless that will lead to growing season length 365 (or 366, + ! if last year was a leap year). This would result in idop==jday, + ! which would invoke the "manually setting sowing_count and + ! sdates_thisyr" code. This would lead to crops never getting + ! harvested. Instead, always harvest the day before idop. + if ((.not. do_harvest) .and. & + (idop(p) > 1 .and. jday == idop(p) - 1) .or. & + (idop(p) == 1 .and. jday == dayspyr)) then + do_harvest = .true. + harvest_reason = HARVEST_REASON_IDOPTOMORROW + else if (do_harvest) then + harvest_reason = HARVEST_REASON_SOWTOMORROW + end if else if (did_plant_prescribed_today) then ! Do not harvest on the day this growing season began; @@ -2218,25 +2327,17 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & ! Original harvest rule do_harvest = hui(p) >= gddmaturity(p) .or. idpp >= mxmat - ! Always harvest the day before the next prescribed sowing, if still alive. - ! WARNING: This implementation assumes that sowing dates don't change over time! + ! Always harvest the day before the next prescribed sowing date, if still alive. + ! WARNING: This implementation assumes that prescribed sowing dates don't change over time! ! In order to avoid this, you'd have to read this year's AND next year's prescribed ! sowing dates. - ! WARNING: This implementation assumes that all patches use prescribed sowing dates. - if (use_cropcal_rx_sdates) then - will_plant_prescribed_tomorrow = (jday == next_rx_sdate(p) - 1) .or. & - (crop_inst%sdates_thisyr_patch(p,1) == 1 .and. & - jday == dayspyr) - else - will_plant_prescribed_tomorrow = .false. - end if - do_harvest = do_harvest .or. will_plant_prescribed_tomorrow + do_harvest = do_harvest .or. do_plant_prescribed_tomorrow if (hui(p) >= gddmaturity(p)) then harvest_reason = HARVEST_REASON_MATURE else if (idpp >= mxmat) then harvest_reason = HARVEST_REASON_MAXSEASLENGTH - else if (will_plant_prescribed_tomorrow) then + else if (do_plant_prescribed_tomorrow) then harvest_reason = HARVEST_REASON_SOWTOMORROW end if endif @@ -2286,7 +2387,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & crop_inst%gddaccum_thisyr_patch(p, harvest_count(p)) = crop_inst%gddaccum_patch(p) crop_inst%hui_thisyr_patch(p, harvest_count(p)) = hui(p) crop_inst%sowing_reason_perharv_patch(p, harvest_count(p)) = real(crop_inst%sowing_reason_patch(p), r8) - crop_inst%sowing_reason_patch(p) = -1 + crop_inst%sowing_reason_patch(p) = -1 ! "Reason for most recent sowing of this patch." So in the line above we save, and here we reset. crop_inst%harvest_reason_thisyr_patch(p, harvest_count(p)) = harvest_reason endif @@ -2354,6 +2455,11 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & endif end if ! croplive + ! At the end of the sowing window, AFTER we've done everything crop-related, set this to false + if (is_end_sowing_window .and. is_end_curr_day()) then + crop_inst%sown_in_this_window(p) = .false. + end if + end do ! prognostic crops loop end associate @@ -2503,7 +2609,7 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & ! !USES: use clm_varctl , only : use_c13, use_c14 - use clm_varctl , only : use_cropcal_rx_sdates, use_cropcal_rx_cultivar_gdds, use_cropcal_streams + use clm_varctl , only : use_cropcal_rx_cultivar_gdds, use_cropcal_streams use clm_varcon , only : c13ratio, c14ratio use clm_varpar , only : mxsowings use pftconMod , only : ntmp_corn, nswheat, nwwheat, ntmp_soybean @@ -2543,7 +2649,6 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & ivt => patch%itype , & ! Input: [integer (:) ] patch vegetation type croplive => crop_inst%croplive_patch , & ! Output: [logical (:) ] Flag, true if planted, not harvested harvdate => crop_inst%harvdate_patch , & ! Output: [integer (:) ] harvest date - next_rx_sdate => crop_inst%next_rx_sdate_patch , & ! Inout: [integer (:) ] prescribed sowing date of next growing season this year sowing_count => crop_inst%sowing_count , & ! Inout: [integer (:) ] number of sowing events this year for this patch sowing_reason => crop_inst%sowing_reason_thisyr_patch , & ! Output: [real(r8) (:) ] reason for each sowing this year for this patch gddmaturity => cnveg_state_inst%gddmaturity_patch , & ! Output: [real(r8) (:) ] gdd needed to harvest @@ -2569,16 +2674,12 @@ 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(p) = .true. idop(p) = jday iyop(p) = kyr harvdate(p) = NOT_Harvested sowing_count(p) = sowing_count(p) + 1 - if (use_cropcal_rx_sdates .and. sowing_count(p) < mxsowings) then - next_rx_sdate(p) = crop_inst%rx_sdates_thisyr_patch(p, sowing_count(p)+1) - else - next_rx_sdate(p) = -1 - endif crop_inst%sdates_thisyr_patch(p,sowing_count(p)) = real(jday, r8) this_sowing_reason = 0._r8 diff --git a/src/biogeochem/CropType.F90 b/src/biogeochem/CropType.F90 index f84c4c1821..29c4717ab3 100644 --- a/src/biogeochem/CropType.F90 +++ b/src/biogeochem/CropType.F90 @@ -47,10 +47,13 @@ module CropType character(len=20) :: baset_mapping real(r8) :: baset_latvary_intercept real(r8) :: baset_latvary_slope - integer , pointer :: next_rx_sdate_patch (:) ! prescribed sowing date for the next growing season this year - integer , pointer :: rx_sdates_thisyr_patch (:,:) ! all prescribed sowing dates for this patch this year (day of year) [patch, mxsowings] + 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 :: rx_swindow_starts_thisyr_patch(:,:) ! all prescribed sowing window start dates for this patch this year (day of year) [patch, mxsowings] + integer , pointer :: rx_swindow_ends_thisyr_patch (:,:) ! all prescribed sowing window end dates for this patch this year (day of year) [patch, mxsowings] real(r8), pointer :: rx_cultivar_gdds_thisyr_patch (:,:) ! all cultivar GDD targets for this patch this year (ddays) [patch, mxsowings] real(r8), pointer :: sdates_thisyr_patch (:,:) ! all actual sowing dates for this patch this year (day of year) [patch, mxsowings] + real(r8), pointer :: swindow_starts_thisyr_patch(:,:) ! all sowing window start dates for this patch this year (day of year) [patch, mxsowings] + real(r8), pointer :: swindow_ends_thisyr_patch (:,:) ! all sowing window end dates for this patch this year (day of year) [patch, mxsowings] real(r8), pointer :: sdates_perharv_patch (:,:) ! all actual sowing dates for crops *harvested* this year (day of year) [patch, mxharvests] real(r8), pointer :: syears_perharv_patch (:,:) ! all actual sowing years for crops *harvested* this year (day of year) [patch, mxharvests] real(r8), pointer :: hdates_thisyr_patch (:,:) ! all actual harvest dates for this patch this year (day of year) [patch, mxharvests] @@ -227,10 +230,13 @@ 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%next_rx_sdate_patch(begp:endp)) ; this%next_rx_sdate_patch(:) = -1 - allocate(this%rx_sdates_thisyr_patch(begp:endp,1:mxsowings)) ; this%rx_sdates_thisyr_patch(:,:) = -1 + allocate(this%sown_in_this_window(begp:endp)) ; this%sown_in_this_window(:) = .false. + allocate(this%rx_swindow_starts_thisyr_patch(begp:endp,1:mxsowings)); this%rx_swindow_starts_thisyr_patch(:,:) = -1 + allocate(this%rx_swindow_ends_thisyr_patch(begp:endp,1:mxsowings)) ; this%rx_swindow_ends_thisyr_patch (:,:) = -1 allocate(this%rx_cultivar_gdds_thisyr_patch(begp:endp,1:mxsowings)) ; this%rx_cultivar_gdds_thisyr_patch(:,:) = spval allocate(this%sdates_thisyr_patch(begp:endp,1:mxsowings)) ; this%sdates_thisyr_patch(:,:) = spval + allocate(this%swindow_starts_thisyr_patch(begp:endp,1:mxsowings)) ; this%swindow_starts_thisyr_patch(:,:) = spval + allocate(this%swindow_ends_thisyr_patch (begp:endp,1:mxsowings)) ; this%swindow_ends_thisyr_patch (:,:) = spval allocate(this%sdates_perharv_patch(begp:endp,1:mxharvests)) ; this%sdates_perharv_patch(:,:) = spval allocate(this%syears_perharv_patch(begp:endp,1:mxharvests)) ; this%syears_perharv_patch(:,:) = spval allocate(this%hdates_thisyr_patch(begp:endp,1:mxharvests)) ; this%hdates_thisyr_patch(:,:) = spval @@ -242,9 +248,6 @@ subroutine InitAllocate(this, bounds) allocate(this%sowing_count(begp:endp)) ; this%sowing_count(:) = 0 allocate(this%harvest_count(begp:endp)) ; this%harvest_count(:) = 0 - this%rx_sdates_thisyr_patch(:,:) = -1 - this%rx_cultivar_gdds_thisyr_patch(:,:) = spval - end subroutine InitAllocate !----------------------------------------------------------------------- @@ -302,6 +305,16 @@ subroutine InitHistory(this, bounds) avgflag='I', long_name='actual crop sowing dates; should only be output annually', & ptr_patch=this%sdates_thisyr_patch, default='inactive') + this%swindow_starts_thisyr_patch(begp:endp,:) = spval + call hist_addfld2d (fname='SWINDOW_STARTS', units='day of year', type2d='mxsowings', & + avgflag='I', long_name='crop sowing window start dates; should only be output annually', & + ptr_patch=this%swindow_starts_thisyr_patch, default='inactive') + + this%swindow_ends_thisyr_patch(begp:endp,:) = spval + call hist_addfld2d (fname='SWINDOW_ENDS', units='day of year', type2d='mxsowings', & + avgflag='I', long_name='crop sowing window end dates; should only be output annually', & + ptr_patch=this%swindow_ends_thisyr_patch, default='inactive') + this%sdates_perharv_patch(begp:endp,:) = spval call hist_addfld2d (fname='SDATES_PERHARV', units='day of year', type2d='mxharvests', & avgflag='I', long_name='actual sowing dates for crops harvested this year; should only be output annually', & @@ -586,6 +599,31 @@ subroutine Restart(this, bounds, ncid, cnveg_state_inst, flag) end if deallocate(temp1d) + allocate(temp1d(bounds%begp:bounds%endp)) + if (flag == 'write') then + do p= bounds%begp,bounds%endp + if (this%sown_in_this_window(p)) then + temp1d(p) = 1 + else + temp1d(p) = 0 + end if + end do + end if + call restartvar(ncid=ncid, flag=flag, varname='sown_in_this_window', xtype=ncd_log, & + dim1name='pft', & + long_name='Flag that patch was sown already during the current sowing window', & + interpinic_flag='interp', readvar=readvar, data=temp1d) + if (flag == 'read') then + do p= bounds%begp,bounds%endp + if (temp1d(p) == 1) then + this%sown_in_this_window(p) = .true. + else + this%sown_in_this_window(p) = .false. + end if + end do + end if + deallocate(temp1d) + call restartvar(ncid=ncid, flag=flag, varname='harvdate', xtype=ncd_int, & dim1name='pft', long_name='harvest date', units='jday', nvalid_range=(/1,366/), & interpinic_flag='interp', readvar=readvar, data=this%harvdate_patch) diff --git a/src/biogeochem/test/CNPhenology_test/test_CNPhenology.pf b/src/biogeochem/test/CNPhenology_test/test_CNPhenology.pf index 135364e474..9e06bc74e2 100644 --- a/src/biogeochem/test/CNPhenology_test/test_CNPhenology.pf +++ b/src/biogeochem/test/CNPhenology_test/test_CNPhenology.pf @@ -293,6 +293,212 @@ contains end do end subroutine check_crit_dayl_dependsonlatnveg + @Test + subroutine test_get_swindow_startend(this) + use clm_time_manager, only : get_curr_days_per_year + class(TestCNPhenology), intent(inout) :: this + integer, dimension(3), parameter :: rx_starts = (/1, 150, -1/) + integer, dimension(3), parameter :: rx_ends = (/45, 180, -1/) + integer, parameter :: param_start = 200 + integer, parameter :: param_end = 250 + integer :: w, start_w, end_w + + call get_swindow(1, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(1, start_w) + @assertEqual(45, end_w) + + call get_swindow(45, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(1, start_w) + @assertEqual(45, end_w) + + call get_swindow(149, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(2, w) + @assertEqual(150, start_w) + @assertEqual(180, end_w) + + call get_swindow(175, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(2, w) + @assertEqual(150, start_w) + @assertEqual(180, end_w) + + call get_swindow(229, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(1, start_w) + @assertEqual(45, end_w) + + call get_swindow(get_curr_days_per_year(), rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(1, start_w) + @assertEqual(45, end_w) + end subroutine test_get_swindow_startend + + @Test + subroutine test_get_swindow_endstart(this) + use clm_time_manager, only : get_curr_days_per_year + class(TestCNPhenology), intent(inout) :: this + integer, dimension(3), parameter :: rx_starts = (/360, 150, -1/) + integer, dimension(3), parameter :: rx_ends = (/45, 180, -1/) + integer, parameter :: param_start = 200 + integer, parameter :: param_end = 250 + integer :: w, start_w, end_w + + call get_swindow(1, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(360, start_w) + @assertEqual(45, end_w) + + call get_swindow(45, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(start_w, 360) + @assertEqual(end_w, 45) + + call get_swindow(149, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(2, w) + @assertEqual(150, start_w) + @assertEqual(180, end_w) + + call get_swindow(175, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(2, w) + @assertEqual(150, start_w) + @assertEqual(180, end_w) + + call get_swindow(229, rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(360, start_w) + @assertEqual(45, end_w) + + call get_swindow(get_curr_days_per_year(), rx_starts, rx_ends, param_start, param_end, w, start_w, end_w) + @assertEqual(1, w) + @assertEqual(360, start_w) + @assertEqual(45, end_w) + end subroutine test_get_swindow_endstart + + @Test + subroutine test_was_sown_in_this_window_startend(this) + use clm_time_manager , only : is_doy_in_interval + class(TestCNPhenology), intent(inout) :: this + integer, parameter :: sowing_window_startdate = 84 + integer, parameter :: sowing_window_enddate = 205 + integer :: jday, idop + + jday = 100 + + idop = 90 + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + ! With jday 100 and idop 90 (as well as idop 100 below), if the existing value of sown_in_this_window is false, then even though it LOOKS like it was sown in the current window, it must have been in a previous year. If it had been sown in this window this year, sown_in_this_window would be true. + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 100 + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 101 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + idop = 120 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + idop = 360 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + jday = 300 + + idop = 90 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + idop = 360 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + end subroutine test_was_sown_in_this_window_startend + + @Test + subroutine test_was_sown_in_this_window_endstart(this) + use clm_time_manager , only : is_doy_in_interval + class(TestCNPhenology), intent(inout) :: this + integer, parameter :: sowing_window_startdate = 205 + integer, parameter :: sowing_window_enddate = 84 + integer :: jday, idop + + jday = 300 + + idop = 60 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + idop = 205 + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + ! With jday 300 and idop 205 (as well as jday 70 idop 60 below), if the existing value of sown_in_this_window is false, then even though it LOOKS like it was sown in the current window, it must have been in a previous year. If it had been sown in the actual current window, sown_in_this_window would be true. + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 300 + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 301 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + jday = 70 + + idop = 60 + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 75 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + idop = 301 + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + end subroutine test_was_sown_in_this_window_endstart + + @Test + subroutine test_was_sown_in_this_window_sameday(this) + use clm_time_manager , only : is_doy_in_interval + class(TestCNPhenology), intent(inout) :: this + integer, parameter :: sowing_window_startdate = 205 + integer, parameter :: sowing_window_enddate = 205 + integer :: jday, idop + + ! If today == start == end, we trust whatever the current value of sown_in_this_window is. + jday = 205 + idop = 205 + ! If it's false, then even if idop == jday, it that idop value must be left over from planting in a PREVIOUS year. + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + ! The ONLY way was_sown_in_this_window() should return true is if today == start == end == idop AND the current value is true. + @assertTrue(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + ! There is no other situation where was_sown_in_this_window() should return true. + + jday = 300 + idop = 60 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 205 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 300 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 301 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + jday = 70 + + idop = 60 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + + idop = 75 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + + idop = 301 + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .true.)) + @assertFalse(was_sown_in_this_window(sowing_window_startdate, sowing_window_enddate, jday, idop, .false.)) + end subroutine test_was_sown_in_this_window_sameday end module test_CNPhenology diff --git a/src/cpl/share_esmf/cropcalStreamMod.F90 b/src/cpl/share_esmf/cropcalStreamMod.F90 index 8c2168d01d..46696eeba9 100644 --- a/src/cpl/share_esmf/cropcalStreamMod.F90 +++ b/src/cpl/share_esmf/cropcalStreamMod.F90 @@ -13,8 +13,9 @@ module cropcalStreamMod use decompMod , only : bounds_type use abortutils , only : endrun use clm_varctl , only : iulog - use clm_varctl , only : use_cropcal_rx_sdates, use_cropcal_rx_cultivar_gdds, use_cropcal_streams + use clm_varctl , only : use_cropcal_rx_swindows, use_cropcal_rx_cultivar_gdds, use_cropcal_streams use clm_varpar , only : mxpft + use clm_varpar , only : mxsowings use perf_mod , only : t_startf, t_stopf use spmdMod , only : masterproc, mpicom, iam use pftconMod , only : npcropmin @@ -31,12 +32,15 @@ module cropcalStreamMod ! !PRIVATE MEMBER DATA: integer, allocatable :: g_to_ig(:) ! Array matching gridcell index to data index - type(shr_strdata_type) :: sdat_cropcal_sdate ! sdate input data stream + type(shr_strdata_type) :: sdat_cropcal_swindow_start ! sowing window start input data stream + type(shr_strdata_type) :: sdat_cropcal_swindow_end ! sowing window end input data stream type(shr_strdata_type) :: sdat_cropcal_cultivar_gdds ! sdate input data stream - character(len=CS), allocatable :: stream_varnames_sdate(:) + character(len=CS), allocatable :: stream_varnames_sdate(:) ! used for both start and end dates character(len=CS), allocatable :: stream_varnames_cultivar_gdds(:) integer :: ncft ! Number of crop functional types (excl. generic crops) - character(len=CL) :: stream_fldFileName_sdate ! sdate stream filename to read + logical :: allow_invalid_swindow_inputs ! Fall back on paramfile sowing windows in cases of invalid values in stream_fldFileName_swindow_start and _end? + character(len=CL) :: stream_fldFileName_swindow_start ! sowing window start stream filename to read + character(len=CL) :: stream_fldFileName_swindow_end ! sowing window end stream filename to read character(len=CL) :: stream_fldFileName_cultivar_gdds ! cultivar growing degree-days stream filename to read character(len=*), parameter :: sourcefile = & @@ -81,7 +85,9 @@ subroutine cropcal_init(bounds) stream_year_first_cropcal, & stream_year_last_cropcal, & model_year_align_cropcal, & - stream_fldFileName_sdate, & + allow_invalid_swindow_inputs, & + stream_fldFileName_swindow_start, & + stream_fldFileName_swindow_end, & stream_fldFileName_cultivar_gdds, & stream_meshfile_cropcal @@ -89,8 +95,10 @@ subroutine cropcal_init(bounds) stream_year_first_cropcal = 1 ! first year in stream to use stream_year_last_cropcal = 1 ! last year in stream to use model_year_align_cropcal = 1 ! align stream_year_first_cropcal with this model year + allow_invalid_swindow_inputs = .false. stream_meshfile_cropcal = '' - stream_fldFileName_sdate = '' + stream_fldFileName_swindow_start = '' + stream_fldFileName_swindow_end = '' stream_fldFileName_cultivar_gdds = '' ! Will need modification to work with mxsowings > 1 ncft = mxpft - npcropmin + 1 ! Ignores generic crops @@ -119,7 +127,9 @@ subroutine cropcal_init(bounds) call shr_mpi_bcast(stream_year_first_cropcal , mpicom) call shr_mpi_bcast(stream_year_last_cropcal , mpicom) call shr_mpi_bcast(model_year_align_cropcal , mpicom) - call shr_mpi_bcast(stream_fldFileName_sdate , mpicom) + call shr_mpi_bcast(allow_invalid_swindow_inputs, mpicom) + call shr_mpi_bcast(stream_fldFileName_swindow_start, mpicom) + call shr_mpi_bcast(stream_fldFileName_swindow_end , mpicom) call shr_mpi_bcast(stream_fldFileName_cultivar_gdds, mpicom) call shr_mpi_bcast(stream_meshfile_cropcal , mpicom) @@ -129,7 +139,9 @@ subroutine cropcal_init(bounds) write(iulog,'(a,i8)') ' stream_year_first_cropcal = ',stream_year_first_cropcal write(iulog,'(a,i8)') ' stream_year_last_cropcal = ',stream_year_last_cropcal write(iulog,'(a,i8)') ' model_year_align_cropcal = ',model_year_align_cropcal - write(iulog,'(a,a)' ) ' stream_fldFileName_sdate = ',trim(stream_fldFileName_sdate) + write(iulog,'(a,l1)') ' allow_invalid_swindow_inputs = ',allow_invalid_swindow_inputs + write(iulog,'(a,a)' ) ' stream_fldFileName_swindow_start = ',trim(stream_fldFileName_swindow_start) + write(iulog,'(a,a)' ) ' stream_fldFileName_swindow_end = ',trim(stream_fldFileName_swindow_end) write(iulog,'(a,a)' ) ' stream_fldFileName_cultivar_gdds = ',trim(stream_fldFileName_cultivar_gdds) write(iulog,'(a,a)' ) ' stream_meshfile_cropcal = ',trim(stream_meshfile_cropcal) do n = 1,ncft @@ -139,14 +151,15 @@ subroutine cropcal_init(bounds) write(iulog,*) endif - use_cropcal_rx_sdates = stream_fldFileName_sdate /= '' + ! CLMBuildNamelist checks that both start and end files are provided if either is + use_cropcal_rx_swindows = stream_fldFileName_swindow_start /= '' use_cropcal_rx_cultivar_gdds = stream_fldFileName_cultivar_gdds /= '' - use_cropcal_streams = use_cropcal_rx_sdates .or. use_cropcal_rx_cultivar_gdds + use_cropcal_streams = use_cropcal_rx_swindows .or. use_cropcal_rx_cultivar_gdds - ! Initialize the cdeps data type sdat_cropcal_sdate - ! NOTE: stream_dtlimit 1.5 didn't work for some reason - if (use_cropcal_rx_sdates) then - call shr_strdata_init_from_inline(sdat_cropcal_sdate, & + if (use_cropcal_rx_swindows) then + ! Initialize the cdeps data type sdat_cropcal_swindow_start + ! NOTE: stream_dtlimit 1.5 didn't work for some reason + call shr_strdata_init_from_inline(sdat_cropcal_swindow_start, & my_task = iam, & logunit = iulog, & compname = 'LND', & @@ -155,7 +168,7 @@ subroutine cropcal_init(bounds) stream_meshfile = trim(stream_meshfile_cropcal), & stream_lev_dimname = 'null', & stream_mapalgo = trim(cropcal_mapalgo), & - stream_filenames = (/trim(stream_fldFileName_sdate)/), & + stream_filenames = (/trim(stream_fldFileName_swindow_start)/), & stream_fldlistFile = stream_varnames_sdate, & stream_fldListModel = stream_varnames_sdate, & stream_yearFirst = stream_year_first_cropcal, & @@ -165,7 +178,34 @@ subroutine cropcal_init(bounds) stream_taxmode = 'extend', & stream_dtlimit = 1.0e30_r8, & stream_tintalgo = cropcal_tintalgo, & - stream_name = 'sowing date data', & + stream_name = 'sowing window start data', & + rc = rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Initialize the cdeps data type sdat_cropcal_swindow_end + ! NOTE: stream_dtlimit 1.5 didn't work for some reason + call shr_strdata_init_from_inline(sdat_cropcal_swindow_end, & + my_task = iam, & + logunit = iulog, & + compname = 'LND', & + model_clock = model_clock, & + model_mesh = mesh, & + stream_meshfile = trim(stream_meshfile_cropcal), & + stream_lev_dimname = 'null', & + stream_mapalgo = trim(cropcal_mapalgo), & + stream_filenames = (/trim(stream_fldFileName_swindow_end)/), & + stream_fldlistFile = stream_varnames_sdate, & + stream_fldListModel = stream_varnames_sdate, & + stream_yearFirst = stream_year_first_cropcal, & + stream_yearLast = stream_year_last_cropcal, & + stream_yearAlign = model_year_align_cropcal, & + stream_offset = cropcal_offset, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = cropcal_tintalgo, & + stream_name = 'sowing window end data', & rc = rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then call ESMF_Finalize(endflag=ESMF_END_ABORT) @@ -227,8 +267,12 @@ subroutine cropcal_advance( bounds ) call get_curr_date(year, mon, day, sec) mcdate = year*10000 + mon*100 + day - if (use_cropcal_rx_sdates) then - call shr_strdata_advance(sdat_cropcal_sdate, ymd=mcdate, tod=sec, logunit=iulog, istr='cropcaldyn', rc=rc) + if (use_cropcal_rx_swindows) then + call shr_strdata_advance(sdat_cropcal_swindow_start, ymd=mcdate, tod=sec, logunit=iulog, istr='cropcaldyn', rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + call shr_strdata_advance(sdat_cropcal_swindow_end, ymd=mcdate, tod=sec, logunit=iulog, istr='cropcaldyn', rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then call ESMF_Finalize(endflag=ESMF_END_ABORT) end if @@ -260,6 +304,8 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! !USES: use CropType , only : crop_type use PatchType , only : patch + use clm_time_manager, only : get_curr_days_per_year + use pftconMod , only : pftname use dshr_methods_mod , only : dshr_fldbun_getfldptr ! ! !ARGUMENTS: @@ -272,15 +318,23 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! !LOCAL VARIABLES: integer :: ivt, p, ip, ig integer :: nc, fp + integer :: dayspyr integer :: n, g integer :: lsize integer :: rc - real(r8), pointer :: dataptr1d_sdate(:) - real(r8), pointer :: dataptr2d_sdate(:,:) + real(r8), pointer :: dataptr1d_swindow_start(:) + real(r8), pointer :: dataptr1d_swindow_end (:) real(r8), pointer :: dataptr1d_cultivar_gdds(:) + real(r8), pointer :: dataptr2d_swindow_start(:,:) + real(r8), pointer :: dataptr2d_swindow_end (:,:) real(r8), pointer :: dataptr2d_cultivar_gdds(:,:) !----------------------------------------------------------------------- + associate( & + starts => crop_inst%rx_swindow_starts_thisyr_patch, & + ends => crop_inst%rx_swindow_ends_thisyr_patch & + ) + SHR_ASSERT_FL( (lbound(g_to_ig,1) <= bounds%begg ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(g_to_ig,1) >= bounds%endg ), sourcefile, __LINE__) @@ -288,30 +342,43 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! Place all data from each type into a temporary 2d array lsize = bounds%endg - bounds%begg + 1 - ! Read prescribed sowing dates from input files - allocate(dataptr2d_sdate(lsize, ncft)) - if (use_cropcal_rx_sdates) then + dayspyr = get_curr_days_per_year() + + ! Read prescribed sowing window start dates from input files + allocate(dataptr2d_swindow_start(lsize, ncft)) + dataptr2d_swindow_start(:,:) = -1._r8 + allocate(dataptr2d_swindow_end (lsize, ncft)) + dataptr2d_swindow_end(:,:) = -1._r8 + if (use_cropcal_rx_swindows) then ! Starting with npcropmin will skip generic crops do n = 1, ncft - call dshr_fldbun_getFldPtr(sdat_cropcal_sdate%pstrm(1)%fldbun_model, trim(stream_varnames_sdate(n)), & - fldptr1=dataptr1d_sdate, rc=rc) + call dshr_fldbun_getFldPtr(sdat_cropcal_swindow_start%pstrm(1)%fldbun_model, trim(stream_varnames_sdate(n)), & + fldptr1=dataptr1d_swindow_start, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + call dshr_fldbun_getFldPtr(sdat_cropcal_swindow_end%pstrm(1)%fldbun_model, trim(stream_varnames_sdate(n)), & + fldptr1=dataptr1d_swindow_end, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then call ESMF_Finalize(endflag=ESMF_END_ABORT) end if ! Note that the size of dataptr1d includes ocean points so it will be around 3x larger than lsize ! So an explicit loop is required here do g = 1,lsize - - ! If read-in value is invalid, allow_unprescribed_planting in CropPhenology() - if (dataptr1d_sdate(g) <= 0 .or. dataptr1d_sdate(g) > 365) then - dataptr1d_sdate(g) = -1 + + ! If read-in value is invalid, set to -1. Will be handled later in this subroutine. + if (dataptr1d_swindow_start(g) <= 0 .or. dataptr1d_swindow_start(g) > dayspyr & + .or. dataptr1d_swindow_end(g) <= 0 .or. dataptr1d_swindow_end(g) > dayspyr) then + dataptr1d_swindow_start(g) = -1 + dataptr1d_swindow_end (g) = -1 end if - - dataptr2d_sdate(g,n) = dataptr1d_sdate(g) + + dataptr2d_swindow_start(g,n) = dataptr1d_swindow_start(g) + dataptr2d_swindow_end (g,n) = dataptr1d_swindow_end (g) end do end do - - ! Set rx_sdate for each gridcell/patch combination + + ! Set sowing window for each gridcell/patch combination do fp = 1, num_pcropp p = filter_pcropp(fp) ivt = patch%itype(p) @@ -320,27 +387,50 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) n = ivt - npcropmin + 1 ! vegetated pft ig = g_to_ig(patch%gridcell(p)) - crop_inst%rx_sdates_thisyr_patch(p,1) = dataptr2d_sdate(ig,n) - - ! Sanity check: Should only read in valid values - if (crop_inst%rx_sdates_thisyr_patch(p,1) > 365) then - write(iulog,'(a,i0,a,i0)') 'cropcal_interp(): Crop patch (ivt ',ivt,') has dataptr2d prescribed sowing date ',& - crop_inst%rx_sdates_thisyr_patch(p,1) - call ESMF_Finalize(endflag=ESMF_END_ABORT) - end if - - ! Only for first sowing date of the year - ! The conditional here is to ensure nothing weird happens if it's called incorrectly on day 365 - if (crop_inst%sdates_thisyr_patch(p,1) <= 0) then - crop_inst%next_rx_sdate_patch(p) = crop_inst%rx_sdates_thisyr_patch(p,1) - end if + starts(p,1) = dataptr2d_swindow_start(ig,n) + ends(p,1) = dataptr2d_swindow_end (ig,n) else - write(iulog,'(a,i0)') 'cropcal_interp(), rx_sdates: Crop patch has ivt ',ivt + write(iulog,'(a,i0)') 'cropcal_interp(), prescribed sowing windows: Crop patch has ivt ',ivt call ESMF_Finalize(endflag=ESMF_END_ABORT) endif end do - end if ! use_cropcal_rx_sdates - deallocate(dataptr2d_sdate) + + ! Ensure that, if mxsowings > 1, sowing windows are ordered such that ENDS are monotonically increasing. This is necessary because of how get_swindow() works. + if (mxsowings > 1) then + if (any(ends(:,2:mxsowings) <= ends(:,1:mxsowings-1) .and. & + ends(:,2:mxsowings) >= 1)) then + write(iulog, *) 'Sowing window inputs must be ordered such that end dates are monotonically increasing.' + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + end if + + ! Handle invalid sowing window values + if (any(starts < 1 .or. ends < 1)) then + ! Fail if not allowing fallback to paramfile sowing windows + if ((.not. allow_invalid_swindow_inputs) .and. any(all(starts < 1, dim=2) .and. patch%wtgcell > 0._r8 .and. patch%itype >= npcropmin)) then + write(iulog, *) 'At least one crop in one gridcell has invalid prescribed sowing window start date(s). To ignore and fall back to paramfile sowing windows, set allow_invalid_swindow_inputs to .true.' + write(iulog, *) 'Affected crops:' + do ivt = npcropmin, mxpft + do fp = 1, num_pcropp + p = filter_pcropp(fp) + if (ivt == patch%itype(p) .and. patch%wtgcell(p) > 0._r8 .and. all(starts(p,:) < 1)) then + write(iulog, *) ' ',pftname(ivt),' (',ivt,')' + exit ! Stop looking for patches of this type + end if + end do + end do + call ESMF_Finalize(endflag=ESMF_END_ABORT) + + ! Fail if a sowing window start date is prescribed without an end date (or vice versa) + else if (any((starts >= 1 .and. ends < 1) .or. (starts < 1 .and. ends >= 1))) then + write(iulog, *) 'Every prescribed sowing window start date must have a corresponding end date.' + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + end if + + end if ! use_cropcal_rx_swindows + deallocate(dataptr2d_swindow_start) + deallocate(dataptr2d_swindow_end) allocate(dataptr2d_cultivar_gdds(lsize, ncft)) if (use_cropcal_rx_cultivar_gdds) then @@ -370,10 +460,6 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) do fp = 1, num_pcropp p = filter_pcropp(fp) -! if (.not. patch%active(p)) then -! continue -! end if - ivt = patch%itype(p) ! Will skip generic crops if (ivt >= npcropmin) then @@ -404,6 +490,8 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) deallocate(dataptr2d_cultivar_gdds) + end associate + end subroutine cropcal_interp end module cropcalStreamMod diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 95c800504b..acefe4acda 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -352,7 +352,7 @@ module clm_varctl !---------------------------------------------------------- logical, public :: use_cropcal_streams = .false. - logical, public :: use_cropcal_rx_sdates = .false. + logical, public :: use_cropcal_rx_swindows = .false. logical, public :: use_cropcal_rx_cultivar_gdds = .false. !---------------------------------------------------------- diff --git a/src/utils/clm_time_manager.F90 b/src/utils/clm_time_manager.F90 index 5c65f5decd..955d98057a 100644 --- a/src/utils/clm_time_manager.F90 +++ b/src/utils/clm_time_manager.F90 @@ -42,6 +42,7 @@ module clm_time_manager get_prev_calday, &! return calendar day at beginning of current timestep get_calday, &! return calendar day from input date get_calendar, &! return calendar + get_doy_tomorrow, &! return next day of year get_average_days_per_year,&! return the average number of days per year for the given calendar get_curr_days_per_year, &! return the days per year for year as of the end of the current time step get_prev_days_per_year, &! return the days per year for year as of the beginning of the current time step @@ -1274,6 +1275,34 @@ end function get_calendar !========================================================================================= + function get_doy_tomorrow(doy_today) result(doy_tomorrow) + + !--------------------------------------------------------------------------------- + ! Given a day of the year (doy_today), return the next day of the year + + integer, intent(in) :: doy_today + integer :: doy_tomorrow + integer :: days_in_year + character(len=*), parameter :: sub = 'clm::get_doy_tomorrow' + + ! Use get_prev_days_per_year() instead of get_curr_days_per_year() because the latter, in the last timestep of a year, actually returns the number of days in the NEXT year. + days_in_year = get_prev_days_per_year() + + if ( doy_today < 1 .or. doy_today > days_in_year )then + write(iulog,*) 'doy_today = ', doy_today + write(iulog,*) 'days_in_year = ', days_in_year + call shr_sys_abort( sub//': error doy_today out of range' ) + end if + + if (doy_today == days_in_year) then + doy_tomorrow = 1 + else + doy_tomorrow = doy_today + 1 + end if + end function get_doy_tomorrow + + !========================================================================================= + real(r8) function get_average_days_per_year() !--------------------------------------------------------------------------------- diff --git a/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf b/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf index 78565fd54d..df8a59de4b 100644 --- a/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf +++ b/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf @@ -15,7 +15,10 @@ module test_clm_time_manager save real(r8), parameter :: tol = 1.e-13_r8 - integer, parameter :: dtime = 1800 + integer, parameter :: secs_in_day_int = 86400 + real(r8), parameter :: secs_in_day_r8 = 86400._r8 + integer, parameter :: dtime_int = 1800 + real(r8), parameter :: dtime_r8 = 1800._r8 @TestCase type, extends(TestCase) :: TestTimeManager @@ -42,11 +45,11 @@ contains class(TestTimeManager), intent(inout) :: this integer :: step_size - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) step_size = get_step_size() - @assertEqual(dtime, step_size) + @assertEqual(dtime_int, step_size) end subroutine getStepSize_returnsCorrectValue @Test @@ -54,7 +57,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: calday - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) calday = get_calday(101, 0) @@ -66,7 +69,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: calday - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar=.true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) calday = get_calday(41231, 43200) @@ -78,7 +81,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: calday - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar=.true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) calday = get_calday(41231, 43200, reuse_day_365_for_day_366=.true.) @@ -90,7 +93,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: calday - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2000, mon=1, day=1, tod=0) @@ -104,7 +107,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: calday - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar=.true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) call set_date(yr=2000, mon=12, day=31, tod=43200) @@ -118,7 +121,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: calday - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar=.true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) call set_date(yr=2000, mon=12, day=31, tod=43200) @@ -133,12 +136,12 @@ contains real(r8) :: calday real(r8) :: calday_expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2000, mon=1, day=1, tod=0) calday = get_prev_calday() - calday_expected = 366._r8 - dtime/86400._r8 + calday_expected = 366._r8 - dtime_int/secs_in_day_r8 @assertEqual(calday_expected, calday, tolerance=tol) end subroutine getPrevCalday_jan1Time0_returnsCorrectValue @@ -149,12 +152,12 @@ contains real(r8) :: calday real(r8) :: calday_expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2000, mon=1, day=2, tod=0) calday = get_prev_calday() - calday_expected = 2._r8 - dtime/86400._r8 + calday_expected = 2._r8 - dtime_int/secs_in_day_r8 @assertEqual(calday_expected, calday, tolerance=tol) end subroutine getPrevCalday_jan2Time0_returnsCorrectValue @@ -164,7 +167,7 @@ contains class(TestTimeManager), intent(inout) :: this real(r8) :: yearfrac - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2, mon=1, day=1, tod=0) @@ -179,7 +182,7 @@ contains real(r8) :: yearfrac real(r8) :: yearfrac_expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2, mon=3, day=1, tod=43200) @@ -195,7 +198,7 @@ contains real(r8) :: yearfrac real(r8) :: yearfrac_expected - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar=.true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) call set_date(yr=2000, mon=12, day=31, tod=43200) @@ -209,16 +212,15 @@ contains subroutine getPrevYearfrac_atYearBoundary_returnsLargeValue(this) class(TestTimeManager), intent(inout) :: this real(r8) :: yearfrac - integer, parameter :: secs_in_day = 86400 real(r8) :: yearfrac_expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2, mon=1, day=1, tod=0) yearfrac = get_prev_yearfrac() - yearfrac_expected = (365._r8 - real(dtime, r8) / real(secs_in_day, r8)) / 365._r8 + yearfrac_expected = (365._r8 - dtime_r8 / secs_in_day_r8) / 365._r8 @assertEqual(yearfrac_expected, yearfrac) end subroutine getPrevYearfrac_atYearBoundary_returnsLargeValue @@ -226,16 +228,15 @@ contains subroutine getPrevYearfrac_inMiddleOfYear_returnsCorrectValue(this) class(TestTimeManager), intent(inout) :: this real(r8) :: yearfrac - integer, parameter :: secs_in_day = 86400 real(r8) :: yearfrac_expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2, mon=3, day=1, tod=43200) yearfrac = get_prev_yearfrac() - yearfrac_expected = (59.5_r8 - real(dtime, r8) / real(secs_in_day, r8)) / 365._r8 + yearfrac_expected = (59.5_r8 - dtime_r8 / secs_in_day_r8) / 365._r8 @assertEqual(yearfrac_expected, yearfrac) end subroutine getPrevYearfrac_inMiddleOfYear_returnsCorrectValue @@ -243,16 +244,15 @@ contains subroutine getPrevYearfrac_leapYearInMiddleOfYear_returnsCorrectValue(this) class(TestTimeManager), intent(inout) :: this real(r8) :: yearfrac - integer, parameter :: secs_in_day = 86400 real(r8) :: yearfrac_expected - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar = .true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar = .true.) call set_date(yr=2000, mon=3, day=1, tod=43200) yearfrac = get_prev_yearfrac() - yearfrac_expected = (60.5_r8 - real(dtime, r8) / real(secs_in_day, r8)) / 366._r8 + yearfrac_expected = (60.5_r8 - dtime_r8 / secs_in_day_r8) / 366._r8 @assertEqual(yearfrac_expected, yearfrac) end subroutine getPrevYearfrac_leapYearInMiddleOfYear_returnsCorrectValue @@ -261,10 +261,9 @@ contains ! This ensures that the correct year is used in determining the number of days in the year class(TestTimeManager), intent(inout) :: this real(r8) :: yearfrac - integer, parameter :: secs_in_day = 86400 real(r8) :: yearfrac_expected - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar = .true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar = .true.) call set_date(yr=2000, mon=1, day=1, tod=0) @@ -272,7 +271,7 @@ contains ! In the following, note that we have 365 and not 366, because the prev_yearfrac uses ! year 1999, which is not a leap year: - yearfrac_expected = (365._r8 - real(dtime, r8) / real(secs_in_day, r8)) / 365._r8 + yearfrac_expected = (365._r8 - dtime_r8 / secs_in_day_r8) / 365._r8 @assertEqual(yearfrac_expected, yearfrac) end subroutine getPrevYearfrac_leapYearAtYearBoundary_returnsCorrectValue @@ -281,7 +280,7 @@ contains class(TestTimeManager), intent(inout) :: this integer :: nstep - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) nstep = get_nstep() @@ -294,7 +293,7 @@ contains integer, parameter :: expected_nstep = 3 integer :: nstep - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_nstep(expected_nstep) @@ -308,9 +307,9 @@ contains class(TestTimeManager), intent(inout) :: this logical :: is_beg - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) - call set_date(yr=2, mon=1, day=1, tod=dtime) + call set_date(yr=2, mon=1, day=1, tod=dtime_int) is_beg = is_beg_curr_year() @@ -322,9 +321,9 @@ contains class(TestTimeManager), intent(inout) :: this logical :: is_beg - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) - call set_date(yr=2, mon=1, day=2, tod=dtime) + call set_date(yr=2, mon=1, day=2, tod=dtime_int) is_beg = is_beg_curr_year() @@ -336,7 +335,7 @@ contains class(TestTimeManager), intent(inout) :: this logical :: is_end - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2, mon=1, day=1, tod=0) @@ -350,7 +349,7 @@ contains class(TestTimeManager), intent(inout) :: this logical :: is_end - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_date(yr=2, mon=1, day=2, tod=0) @@ -364,7 +363,7 @@ contains class(TestTimeManager), intent(inout) :: this logical :: is_first - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) is_first = is_first_step() @@ -376,7 +375,7 @@ contains class(TestTimeManager), intent(inout) :: this logical :: is_first - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) call set_nstep(1) @@ -390,7 +389,7 @@ contains class(TestTimeManager), intent(inout) :: this integer :: nstep - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) nstep = 100 call set_nstep(nstep) @@ -419,7 +418,7 @@ contains real(r8) :: londeg integer :: expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) ! Check for local noon at Greenich londeg = 0.0_r8 @@ -446,7 +445,7 @@ contains integer :: secs real(r8) :: londeg - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) londeg = 0.0_r8 do while ( londeg <= 360.0_r8 ) @@ -470,16 +469,16 @@ contains real(r8) :: londeg integer :: expected - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) ! Check for local noon at Greenich for 1 time step ahead londeg = 0.0_r8 - secs = 3600*12 + dtime + secs = 3600*12 + dtime_int call set_date(yr=2018, mon=9, day=3, tod=secs) - expected = secs - dtime - @assertEqual( expected, get_local_time( londeg, offset=-dtime ) ) + expected = secs - dtime_int + @assertEqual( expected, get_local_time( londeg, offset=-dtime_int ) ) londeg = 360.0_r8 - @assertEqual( expected, get_local_time( londeg, offset=-dtime ) ) + @assertEqual( expected, get_local_time( londeg, offset=-dtime_int ) ) end subroutine check_local_time_woffset @@ -490,7 +489,7 @@ contains integer :: secs real(r8) :: londeg - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) ! Check for local noon at Greenich will be true from 11 to 1pm londeg = 0.0_r8 @@ -528,7 +527,7 @@ contains integer :: secs logical :: check - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) londeg = 0.0_r8 secs = get_local_time( londeg ) @@ -549,7 +548,7 @@ contains real(r8) :: londeg integer :: secs - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) londeg = -200.0_r8 secs = get_local_time( londeg ) @@ -567,7 +566,7 @@ contains real(r8) :: londeg integer :: secs - call unittest_timemgr_setup(dtime=dtime) + call unittest_timemgr_setup(dtime=dtime_int) londeg = 400.0_r8 secs = get_local_time( londeg ) @@ -627,7 +626,7 @@ contains integer :: start, end - call unittest_timemgr_setup(dtime=dtime, use_gregorian_calendar=.true.) + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) start = 100 ! April 10 end = 300 ! October 27 @@ -641,7 +640,7 @@ contains @assertFalse(is_today_in_doy_interval(start, end)) ! Test first timestep of interval - call set_date(yr=2009, mon=4, day=10, tod=dtime) + call set_date(yr=2009, mon=4, day=10, tod=dtime_int) @assertTrue(is_today_in_doy_interval(start, end)) ! Test well within interval @@ -653,9 +652,58 @@ contains @assertTrue(is_today_in_doy_interval(start, end)) ! Test first timestep after interval - call set_date(yr=2009, mon=10, day=28, tod=dtime) + call set_date(yr=2009, mon=10, day=28, tod=dtime_int) @assertFalse(is_today_in_doy_interval(start, end)) end subroutine check_is_today_in_doy_interval + @Test + subroutine check_get_doy_tomorrow_noleap(this) + class(TestTimeManager), intent(inout) :: this + character(len=256) :: expected_msg + integer :: doy_tomorrow ! Dummy needed for exception tests + integer :: days_in_year + + ! We don't care about the actual date here; we just want to enable get_prev_days_per_year() + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) + call set_date(yr=2009, mon=10, day=28, tod=dtime_int) + + ! Use get_prev_days_per_year() instead of get_curr_days_per_year() because the latter, in the last timestep of a year, actually returns the number of days in the NEXT year. + days_in_year = get_prev_days_per_year() + + @assertEqual(get_doy_tomorrow(1), 2) + @assertEqual(get_doy_tomorrow(150), 151) + @assertEqual(get_doy_tomorrow(days_in_year), 1) + + doy_tomorrow = get_doy_tomorrow(days_in_year + 1) + expected_msg = endrun_msg("clm::get_doy_tomorrow: error doy_today out of range" ) + @assertExceptionRaised(expected_msg) + + end subroutine check_get_doy_tomorrow_noleap + + @Test + subroutine check_get_doy_tomorrow_leap(this) + class(TestTimeManager), intent(inout) :: this + character(len=256) :: expected_msg + integer :: doy_tomorrow ! Dummy needed for exception tests + integer :: days_in_year + + call unittest_timemgr_setup(dtime=dtime_int, use_gregorian_calendar=.true.) + + ! Ensure that things work even in the last timestep of a leap year... + call set_date(yr=2008, mon=12, day=31, tod=secs_in_day_int - dtime_int) + ! Use get_prev_days_per_year() instead of get_curr_days_per_year() because the latter, in the last timestep of a year, actually returns the number of days in the NEXT year. + days_in_year = get_prev_days_per_year() + @assertEqual(get_doy_tomorrow(days_in_year), 1) + doy_tomorrow = get_doy_tomorrow(days_in_year + 1) + expected_msg = endrun_msg("clm::get_doy_tomorrow: error doy_today out of range") + @assertExceptionRaised(expected_msg) + + ! ... as well as in the first timestep after a leap year + call set_date(yr=2009, mon=1, day=1, tod=0) + @assertEqual(get_doy_tomorrow(1), 2) + + end subroutine check_get_doy_tomorrow_leap + + end module test_clm_time_manager