diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index ccb310a4f9..b0bf2c5a0c 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -26,6 +26,7 @@ d866510188d26d51bcd6d37239283db690af7e82 0dcd0a3c1abcaffe5529f8d79a6bc34734b195c7 e096358c832ab292ddfd22dd5878826c7c788968 475831f0fb0e31e97f630eac4e078c886558b61c +fd5f177131d63d39e79a13918390bdfb642d781e # Ran SystemTests and python/ctsm through black python formatter 5364ad66eaceb55dde2d3d598fe4ce37ac83a93c 8056ae649c1b37f5e10aaaac79005d6e3a8b2380 @@ -49,3 +50,4 @@ aa04d1f7d86cc2503b98b7e2b2d84dbfff6c316b 665cf86102e09b4c4c5a140700676dca23bc55a9 1a49e547ba3c48fa483f9ae81a8f05adcd6b888c 045d90f1d80f713eb3ae0ac58f6c2352937f1eb0 +753fda3ff0147837231a73c9c728dd9ce47b5997 diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index f92d0811d1..c6a613d449 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1904,9 +1904,10 @@ sub setup_logic_site_specific { $nl->set_variable_value($group, $var, $val); } - if ($nl_flags->{'res'} eq "1x1_smallvilleIA") { + my $res = $nl_flags->{'res'}; + if ($res eq "1x1_smallvilleIA" or $res eq "1x1_cidadinhoBR") { if (! &value_is_true($nl_flags->{'use_cn'}) || ! &value_is_true($nl_flags->{'use_crop'})) { - $log->fatal_error("1x1_smallvilleIA grids must use a compset with CN and CROP turned on."); + $log->fatal_error("${res} grids must use a compset with CN and CROP turned on."); } } @@ -4223,52 +4224,137 @@ sub setup_logic_lai_streams { sub setup_logic_cropcal_streams { my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; - # Set first and last stream years - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_year_first_cropcal', - 'sim_year'=>$nl_flags->{'sim_year'}, - 'sim_year_range'=>$nl_flags->{'sim_year_range'}); - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_year_last_cropcal', - 'sim_year'=>$nl_flags->{'sim_year'}, - 'sim_year_range'=>$nl_flags->{'sim_year_range'}); - - # Set align year, if first and last years are different - if ( $nl->get_value('stream_year_first_cropcal') != - $nl->get_value('stream_year_last_cropcal') ) { - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, - 'model_year_align_cropcal', 'sim_year'=>$nl_flags->{'sim_year'}, - 'sim_year_range'=>$nl_flags->{'sim_year_range'}); - } - # Set up other crop calendar parameters + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'cropcals_rx'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'cropcals_rx_adapt'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_gdd20_seasons'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'flush_gdd20'); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'generate_crop_gdds'); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_mxmat'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_cropcal'); + + # These can't both be true + my $cropcals_rx = $nl->get_value('cropcals_rx') ; + my $cropcals_rx_adapt = $nl->get_value('cropcals_rx_adapt') ; + if (&value_is_true($cropcals_rx) and &value_is_true($cropcals_rx_adapt)) { + $log->fatal_error("cropcals_rx and cropcals_rx_adapt may not both be true" ); + } + + # Add defaults if reading gdd20 seasons from stream files + my $stream_gdd20_seasons = $nl->get_value('stream_gdd20_seasons') ; + if ( &value_is_true($stream_gdd20_seasons)) { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldFileName_gdd20_season_start'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldFileName_gdd20_season_end'); + + # Check + my $gdd20_season_start_file = $nl->get_value('stream_fldFileName_gdd20_season_start') ; + my $gdd20_season_end_file = $nl->get_value('stream_fldFileName_gdd20_season_end') ; + if ( &string_is_undef_or_empty($gdd20_season_start_file) or &string_is_undef_or_empty($gdd20_season_end_file) ) { + $log->message($gdd20_season_start_file); + $log->message($gdd20_season_end_file); + $log->fatal_error("If stream_gdd20_seasons is true, gdd20 season start and end files must be provided." ); + } + } - # Option checks - my $generate_crop_gdds = $nl->get_value('generate_crop_gdds') ; - my $use_mxmat = $nl->get_value('use_mxmat') ; + # Add defaults if using prescribed crop calendars + if ( &value_is_true($cropcals_rx) or &value_is_true($cropcals_rx_adapt) ) { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldFileName_swindow_start'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldFileName_swindow_end'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_cultivar_gdds'); + if ( &value_is_true($cropcals_rx_adapt) ) { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldFileName_gdd20_baseline', 'stream_gdd20_seasons'=>$stream_gdd20_seasons); + } + } + + # Add defaults if using any crop calendar input files 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 $gdd20_baseline_file = $nl->get_value('stream_fldFileName_gdd20_baseline') ; 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 ( !&string_is_undef_or_empty($swindow_start_file) or !&string_is_undef_or_empty($swindow_end_file) or !&string_is_undef_or_empty($gdd_file) or !&string_is_undef_or_empty($gdd20_baseline_file)) { + + # User gave an input file without specifying cropcals_rx or cropcals_rx_adapt = .true. + # Requiring this means nothing to the code, but helps namelist make more sense + if ( !&value_is_true($cropcals_rx) and !&value_is_true($cropcals_rx_adapt) ){ + $log->fatal_error("If providing any crop calendar input file(s), cropcals_rx or cropcals_rx_adapt must be true" ); + } + + # User set cropcals_rx_adapt to true but set stream_fldFileName_gdd20_baseline to empty + if ( &value_is_true($cropcals_rx_adapt) and &string_is_undef_or_empty($gdd20_baseline_file) ) { + $log->fatal_error("If cropcals_rx_adapt is true, stream_fldFileName_gdd20_baseline must be provided" ); + } + + # cropcals_rx_adapt is false but user provided stream_fldFileName_gdd20_baseline + if ( !&value_is_true($cropcals_rx_adapt) and !&string_is_undef_or_empty($gdd20_baseline_file) ) { + $log->fatal_error("If stream_fldFileName_gdd20_baseline provided, cropcals_rx_adapt must be true" ); + } + + # User provided an input file but set mesh file to empty + if ( &string_is_undef_or_empty($mesh_file) ) { + $log->fatal_error("If providing any crop calendar input file(s), you must provide stream_meshfile_cropcal" ); + } + + # Set stream years + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_year_first_cropcal_swindows', + 'sim_year'=>$nl_flags->{'sim_year'}, + 'sim_year_range'=>$nl_flags->{'sim_year_range'}); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_year_last_cropcal_swindows', + 'sim_year'=>$nl_flags->{'sim_year'}, + 'sim_year_range'=>$nl_flags->{'sim_year_range'}); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'model_year_align_cropcal_swindows', + 'sim_year'=>$nl_flags->{'sim_year'}, + 'sim_year_range'=>$nl_flags->{'sim_year_range'}); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_year_first_cropcal_cultivar_gdds', + 'sim_year'=>$nl_flags->{'sim_year'}, + 'sim_year_range'=>$nl_flags->{'sim_year_range'}); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_year_last_cropcal_cultivar_gdds', + 'sim_year'=>$nl_flags->{'sim_year'}, + 'sim_year_range'=>$nl_flags->{'sim_year_range'}); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'model_year_align_cropcal_cultivar_gdds', + 'sim_year'=>$nl_flags->{'sim_year'}, + 'sim_year_range'=>$nl_flags->{'sim_year_range'}); + + # Do not allow maturity requirements to change over time if stream_fldFileName_gdd20_baseline is provided. That would be nonsensical. + if ( $nl->get_value('stream_year_first_cropcal_cultivar_gdds') != + $nl->get_value('stream_year_last_cropcal_cultivar_gdds') + and !&string_is_undef_or_empty($gdd20_baseline_file) ) { + $log->fatal_error("If cropcals_rx_adapt is true (i.e., stream_fldFileName_gdd20_baseline is provided), baseline maturity requirements are allowed to vary over time (i.e., stream_year_first_cropcal_cultivar_gdds and stream_year_last_cropcal_cultivar_gdds must be the same)." ); + } } - if ( $generate_crop_gdds eq '.true.' ) { - if ( $use_mxmat eq '.true.' ) { + + # If running with prescribed crop calendars, certain files must be provided + my $generate_crop_gdds = $nl->get_value('generate_crop_gdds') ; + if ( &value_is_true($cropcals_rx) or &value_is_true($cropcals_rx_adapt) ) { + if ( &string_is_undef_or_empty($swindow_start_file) or &string_is_undef_or_empty($swindow_end_file) ) { + $log->fatal_error("If cropcals_rx or cropcals_rx_adapt is true, sowing window start and end files must be provided. To specify exact sowing dates, use the same file." ); + } + if ( &string_is_undef_or_empty($gdd_file) and (! &value_is_true($generate_crop_gdds)) ){ + $log->fatal_error("If cropcals_rx or cropcals_rx_adapt is true and generate_crop_gdds is false, maturity requirement file stream_fldFileName_cultivar_gdds must be provided" ); + } + } + + # Option checks + if ( &string_is_undef_or_empty($gdd_file) and ! &string_is_undef_or_empty($gdd20_baseline_file) ) { + $log->fatal_error("If not providing stream_fldFileName_cultivar_gdds, don't provide stream_fldFileName_gdd20_baseline"); + } + if ( &value_is_true($generate_crop_gdds) ) { + my $use_mxmat = $nl->get_value('use_mxmat') ; + if ( &value_is_true($use_mxmat) ) { $log->fatal_error("If generate_crop_gdds is true, you must also set use_mxmat to false" ); } - if ( $swindow_start_file eq '' or $swindow_end_file eq '' ) { + if ( &string_is_undef_or_empty($swindow_start_file) or &string_is_undef_or_empty($swindow_end_file) ) { $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 '' ) { + if ( ! &string_is_undef_or_empty($gdd_file) ) { $log->fatal_error("If generate_crop_gdds is true, do not specify stream_fldFileName_cultivar_gdds") } - } - 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") + if ( ! &string_is_undef_or_empty($gdd20_baseline_file) ) { + $log->fatal_error("If generate_crop_gdds is true, do not specify stream_fldFileName_gdd20_baseline") + } } } diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index c94cc47bd3..3594fe5c38 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -1909,6 +1909,8 @@ lnd/clm2/surfdata_esmf/ctsm5.2.0/surfdata_ne3np4.pg3_hist_1850_78pfts_c240216.nc lnd/clm2/surfdata_esmf/ctsm5.2.0/surfdata_C96_hist_1850_78pfts_c240216.nc lnd/clm2/surfdata_esmf/ctsm5.2.0/surfdata_1x1_smallvilleIA_hist_1850_78pfts_c240221.nc + +lnd/clm2/surfdata_esmf/surfdata_1x1_cidadinhoBR_hist_2000_78pfts_c240613.nc lnd/clm2/surfdata_esmf/ctsm5.2.0/surfdata_1x1_brazil_hist_1850_78pfts_c240221.nc @@ -1980,11 +1982,14 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c >lnd/clm2/surfdata_esmf/ctsm5.2.0/landuse.timeseries_C96_SSP2-4.5_1850-2100_78pfts_c240216.nc - lnd/clm2/surfdata_esmf/ctsm5.2.0/landuse.timeseries_1x1_smallvilleIA_SSP2-4.5_1850-1855_78pfts_c240221.nc + + + @@ -2192,10 +2197,34 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c nn nn - -1850 -2100 -1850 + +.false. +.true. +.false. +.false. +.false. +.false. +.false. +2000 +2000 +2000 +2000 +2000 +2000 + + +lnd/clm2/cropdata/calendars/processed/swindow_starts_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/swindow_ends_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/gdds_20230829_161011.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/360x720_120830_ESMFmesh_c20210507_cdf5.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/swindow_starts_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/swindow_ends_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/gdds_20230829_161011.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/20230714_cropcals_pr2_1deg.actually2deg.1980-2009.from_GDDB20.interpd_halfdeg.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/gdd20bl.copied_from.gdds_20230829_161011.v2.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/sdates_ggcmi_crop_calendar_phase3_v1.01_nninterp-hcru_hcru_mt13.2000-2000.20230728_165845.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/hdates_ggcmi_crop_calendar_phase3_v1.01_nninterp-hcru_hcru_mt13.2000-2000.20230728_165845.tweaked_latlons.nc +lnd/clm2/cropdata/calendars/processed/360x720_120830_ESMFmesh_c20210507_cdf5.tweaked_latlons.nc diff --git a/bld/namelist_files/namelist_defaults_overall.xml b/bld/namelist_files/namelist_defaults_overall.xml index 479b2a02b7..5b7ae1bdd9 100644 --- a/bld/namelist_files/namelist_defaults_overall.xml +++ b/bld/namelist_files/namelist_defaults_overall.xml @@ -62,6 +62,7 @@ determine default values for namelists. 1x1_urbanc_alpha 1x1_numaIA 1x1_smallvilleIA +1x1_cidadinhoBR 2000 @@ -110,6 +111,7 @@ determine default values for namelists. test navy test +test gx1v7 diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 8aa324bb16..9cbcc5e692 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1845,19 +1845,44 @@ Mapping method from LAI input file to the model resolution - -First year to loop over for crop calendar data +Flag to enable prescribed crop calendars (sowing window dates and maturity requirement) - -Last year to loop over for crop calendar data +Flag to enable prescribed crop calendars (sowing window dates and maturity requirement), with maturity requirement adaptive based on recent climate - -Simulation year that aligns with stream_year_first_cropcal value +First year to loop over for crop sowing windows + + + +Last year to loop over for crop sowing windows + + + +Simulation year that aligns with stream_year_first_cropcal_swindows value + + + +First year to loop over for crop maturity requirements + + + +Last year to loop over for crop maturity requirements + + + +Simulation year that aligns with stream_year_first_cropcal_cultivar_gdds value + +Filename of input stream data for baseline GDD20 values + + + +Set this to true to read gdd20 accumulation season start and end dates from stream files, rather than using hard-coded hemisphere-specific "warm seasons." + + + +Set this to true to flush the accumulated GDD20 variables as soon as possible. + + + +By default, a value in stream_fldFileName_gdd20_season_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 have such cells fall back to the hard-coded hemisphere-specific "warm seasons." + + + +Filename of input stream data for date (day of year) of start of gdd20 accumulation season. + + + +Filename of input stream data for date (day of year) of end of gdd20 accumulation season. + + Filename of input stream data for crop calendar inputs diff --git a/bld/unit_testers/build-namelist_test.pl b/bld/unit_testers/build-namelist_test.pl index 1533e00d73..c92f428d82 100755 --- a/bld/unit_testers/build-namelist_test.pl +++ b/bld/unit_testers/build-namelist_test.pl @@ -1530,21 +1530,23 @@ sub cat_and_create_namelistinfile { print "==================================================\n"; # Check for crop resolutions -my $crop1850_res = "1x1_smallvilleIA"; -$options = "-bgc bgc -crop -res $crop1850_res -use_case 1850_control -envxml_dir ."; -&make_env_run(); -eval{ system( "$bldnml $options > $tempfile 2>&1 " ); }; -is( $@, '', "$options" ); -$cfiles->checkfilesexist( "$options", $mode ); -$cfiles->shownmldiff( "default", "standard" ); -if ( defined($opts{'compare'}) ) { - $cfiles->doNOTdodiffonfile( "$tempfile", "$options", $mode ); - $cfiles->comparefiles( "$options", $mode, $opts{'compare'} ); -} -if ( defined($opts{'generate'}) ) { - $cfiles->copyfiles( "$options", $mode ); +my @crop1850_res = ( "1x1_smallvilleIA", "1x1_cidadinhoBR" ); +foreach my $res ( @crop1850_res ) { + $options = "-bgc bgc -crop -res $res -use_case 1850_control -envxml_dir ."; + &make_env_run(); + eval{ system( "$bldnml $options > $tempfile 2>&1 " ); }; + is( $@, '', "$options" ); + $cfiles->checkfilesexist( "$options", $mode ); + $cfiles->shownmldiff( "default", "standard" ); + if ( defined($opts{'compare'}) ) { + $cfiles->doNOTdodiffonfile( "$tempfile", "$options", $mode ); + $cfiles->comparefiles( "$options", $mode, $opts{'compare'} ); + } + if ( defined($opts{'generate'}) ) { + $cfiles->copyfiles( "$options", $mode ); + } + &cleanup(); } -&cleanup(); my @crop_res = ( "1x1_numaIA", "4x5", "10x15", "0.9x1.25", "1.9x2.5", "ne3np4.pg3", "ne30np4", "ne30np4.pg3", "C96", "mpasa120" ); foreach my $res ( @crop_res ) { diff --git a/cime_config/SystemTests/rxcropmaturity.py b/cime_config/SystemTests/rxcropmaturity.py index d25bd015ca..fb254c408f 100644 --- a/cime_config/SystemTests/rxcropmaturity.py +++ b/cime_config/SystemTests/rxcropmaturity.py @@ -135,6 +135,7 @@ def _run_phase(self, skip_gen=False): logger.info("RXCROPMATURITY log: modify user_nl files: generate GDDs") self._append_to_user_nl_clm( [ + "stream_fldFileName_cultivar_gdds = ''", "generate_crop_gdds = .true.", "use_mxmat = .false.", " ", @@ -400,12 +401,14 @@ def _run_check_rxboth_run(self, skip_gen): def _modify_user_nl_allruns(self): nl_additions = [ + "cropcals_rx = .true.", + "cropcals_rx_adapt = .false.", "stream_meshfile_cropcal = '{}'".format(self._case.get_value("LND_DOMAIN_MESH")), "stream_fldFileName_swindow_start = '{}'".format(self._sdatefile), "stream_fldFileName_swindow_end = '{}'".format(self._sdatefile), - "stream_year_first_cropcal = 2000", - "stream_year_last_cropcal = 2000", - "model_year_align_cropcal = 2000", + "stream_year_first_cropcal_swindows = 2000", + "stream_year_last_cropcal_swindows = 2000", + "model_year_align_cropcal_swindows = 2000", " ", "! (h1) Annual outputs on sowing or harvest axis", "hist_fincl2 = 'GRAINC_TO_FOOD_PERHARV', 'GRAINC_TO_FOOD_ANN', 'SDATES', 'SDATES_PERHARV', 'SYEARS_PERHARV', 'HDATES', 'GDDHARV_PERHARV', 'GDDACCUM_PERHARV', 'HUI_PERHARV', 'SOWING_REASON_PERHARV', 'HARVEST_REASON_PERHARV'", @@ -442,7 +445,7 @@ def _run_generate_gdds(self, case_gddgen): f"--sdates-file {sdates_file}", f"--hdates-file {hdates_file}", f"--output-dir generate_gdds_out", - f"--skip-crops miscanthus,irrigated_miscanthus", + f"--skip-crops miscanthus,irrigated_miscanthus,switchgrass,irrigated_switchgrass", ] ) stu.run_python_script( diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index d1e0a1a8b8..d5e3e4e378 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -163,6 +163,13 @@ + + + FAIL + #2659 + + + diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 59044feb10..951e0d8d5b 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -863,6 +863,14 @@ + + + + + + + + @@ -1016,6 +1024,7 @@ + @@ -3478,7 +3487,7 @@ - + @@ -3488,9 +3497,10 @@ - + + @@ -3498,10 +3508,11 @@ - + + @@ -3558,5 +3569,82 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/clm/GddGen/README b/cime_config/testdefs/testmods_dirs/clm/GddGen/README new file mode 100644 index 0000000000..3236ca609a --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/GddGen/README @@ -0,0 +1,5 @@ +The GddGen test is set up just like a GDD-Generating run, with two differences: +1) It doesn't include an all-crops-everywhere surface dataset, +2) it doesn't actually run the GDD-generating script, +and +3) it includes some extra outputs. diff --git a/cime_config/testdefs/testmods_dirs/clm/GddGen/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/GddGen/include_user_mods new file mode 100644 index 0000000000..4d75082583 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/GddGen/include_user_mods @@ -0,0 +1,2 @@ +../cropMonthOutput +../oldCropCals \ No newline at end of file diff --git a/cime_config/testdefs/testmods_dirs/clm/GddGen/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/GddGen/user_nl_clm new file mode 100644 index 0000000000..cfde517fd9 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/GddGen/user_nl_clm @@ -0,0 +1,13 @@ +cropcals_rx = .true. +stream_fldFileName_swindow_start = '$DIN_LOC_ROOT/lnd/clm2/cropdata/calendars/processed/sdates_ggcmi_crop_calendar_phase3_v1.01_nninterp-hcru_hcru_mt13.2000-2000.20230728_165845.tweaked_latlons.nc' +stream_fldFileName_swindow_end = '$DIN_LOC_ROOT/lnd/clm2/cropdata/calendars/processed/sdates_ggcmi_crop_calendar_phase3_v1.01_nninterp-hcru_hcru_mt13.2000-2000.20230728_165845.tweaked_latlons.nc' +stream_fldFileName_cultivar_gdds = '' +generate_crop_gdds = .true. +use_mxmat = .false. + +! (h3) Daily outputs for GDD generation and figure-making +hist_fincl4 = 'GDDACCUM', 'GDDHARV' +hist_nhtfrq(4) = -24 +hist_mfilt(4) = 365 +hist_type1d_pertape(4) = 'PFTS' +hist_dov2xy(4) = .false. diff --git a/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/include_user_mods new file mode 100644 index 0000000000..23ea3745e6 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/include_user_mods @@ -0,0 +1 @@ +../crop diff --git a/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/user_nl_clm new file mode 100644 index 0000000000..fdf5a86c26 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/user_nl_clm @@ -0,0 +1,6 @@ +cropcals_rx = .false. +cropcals_rx_adapt = .true. +stream_gdd20_seasons = .true. +flush_gdd20 = .true. +!TODO SSR: Try without this once you have half-degree inputs +allow_invalid_gdd20_season_inputs = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/RxCropCalsNoAdapt/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsNoAdapt/include_user_mods new file mode 100644 index 0000000000..23ea3745e6 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsNoAdapt/include_user_mods @@ -0,0 +1 @@ +../crop diff --git a/cime_config/testdefs/testmods_dirs/clm/RxCropCalsNoAdapt/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsNoAdapt/user_nl_clm new file mode 100644 index 0000000000..625960b389 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsNoAdapt/user_nl_clm @@ -0,0 +1,2 @@ +cropcals_rx = .true. +cropcals_rx_adapt = .false. diff --git a/cime_config/testdefs/testmods_dirs/clm/StreamGdd20Seasons/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/StreamGdd20Seasons/user_nl_clm new file mode 100644 index 0000000000..8b040d9d43 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/StreamGdd20Seasons/user_nl_clm @@ -0,0 +1,2 @@ +stream_gdd20_seasons = .true. +allow_invalid_gdd20_season_inputs = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/crop/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/crop/user_nl_clm index 67042ea01a..aea8e39e6c 100644 --- a/cime_config/testdefs/testmods_dirs/clm/crop/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/crop/user_nl_clm @@ -9,9 +9,12 @@ hist_fincl2 += 'DYN_COL_SOIL_ADJUSTMENTS_C' -! Annual crop variables on per-sowing/per-harvest axes, per PFT. -hist_fincl3 = 'SDATES', 'SDATES_PERHARV', 'SYEARS_PERHARV', 'HDATES', 'GRAINC_TO_FOOD_PERHARV', 'GRAINC_TO_FOOD_ANN', 'GRAINN_TO_FOOD_PERHARV', 'GRAINN_TO_FOOD_ANN', 'GRAINC_TO_SEED_PERHARV', 'GRAINC_TO_SEED_ANN', 'GRAINN_TO_SEED_PERHARV', 'GRAINN_TO_SEED_ANN', 'HDATES', 'GDDHARV_PERHARV', 'GDDACCUM_PERHARV', 'HUI_PERHARV', 'SOWING_REASON_PERHARV', 'HARVEST_REASON_PERHARV', 'SWINDOW_STARTS', 'SWINDOW_ENDS' -hist_nhtfrq(3) = 17520 +! Instantaneous crop variables (including per-sowing/per-harvest axes), per PFT. +! Note that, under normal circumstances, these should only be saved annually. +! That's needed for the mxsowings and mxharvests axes to make sense. +! However, for testing purposes, it makes sense to save more frequently. +hist_fincl3 = 'SDATES', 'SDATES_PERHARV', 'SYEARS_PERHARV', 'HDATES', 'GRAINC_TO_FOOD_PERHARV', 'GRAINC_TO_FOOD_ANN', 'GRAINN_TO_FOOD_PERHARV', 'GRAINN_TO_FOOD_ANN', 'GRAINC_TO_SEED_PERHARV', 'GRAINC_TO_SEED_ANN', 'GRAINN_TO_SEED_PERHARV', 'GRAINN_TO_SEED_ANN', 'HDATES', 'GDDHARV_PERHARV', 'GDDACCUM_PERHARV', 'HUI_PERHARV', 'SOWING_REASON_PERHARV', 'HARVEST_REASON_PERHARV', 'SWINDOW_STARTS', 'SWINDOW_ENDS', 'GDD20_BASELINE', 'GDD20_SEASON_START', 'GDD20_SEASON_END' +hist_nhtfrq(3) = -24 hist_mfilt(3) = 1 hist_type1d_pertape(3) = 'PFTS' hist_dov2xy(3) = .false. diff --git a/cime_config/testdefs/testmods_dirs/clm/cropMonthOutput/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/cropMonthOutput/user_nl_clm index 8f779ed011..18220de5ef 100644 --- a/cime_config/testdefs/testmods_dirs/clm/cropMonthOutput/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/cropMonthOutput/user_nl_clm @@ -1,4 +1,4 @@ - hist_nhtfrq = 0,-240,17520 + hist_nhtfrq = 0,-240,0 hist_mfilt = 1,1,1 ! NOTE slevis (2024/2/23) Adding option for tests to pass. In the long term diff --git a/cime_config/testdefs/testmods_dirs/clm/decStart/README b/cime_config/testdefs/testmods_dirs/clm/decStart/README new file mode 100644 index 0000000000..7cdab6abfd --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/decStart/README @@ -0,0 +1 @@ +Use midDecStart instead of decStart if you want ERP/ERS/etc. tests longer than 2 days to be able to have the split in December instead of January (i.e., before rather than after new year). \ No newline at end of file diff --git a/cime_config/testdefs/testmods_dirs/clm/midDecStart/README b/cime_config/testdefs/testmods_dirs/clm/midDecStart/README new file mode 100644 index 0000000000..7cdab6abfd --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/midDecStart/README @@ -0,0 +1 @@ +Use midDecStart instead of decStart if you want ERP/ERS/etc. tests longer than 2 days to be able to have the split in December instead of January (i.e., before rather than after new year). \ No newline at end of file diff --git a/cime_config/testdefs/testmods_dirs/clm/sowingWindows/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/midDecStart/include_user_mods similarity index 100% rename from cime_config/testdefs/testmods_dirs/clm/sowingWindows/include_user_mods rename to cime_config/testdefs/testmods_dirs/clm/midDecStart/include_user_mods diff --git a/cime_config/testdefs/testmods_dirs/clm/midDecStart/shell_commands b/cime_config/testdefs/testmods_dirs/clm/midDecStart/shell_commands new file mode 100755 index 0000000000..d044ab8c3b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/midDecStart/shell_commands @@ -0,0 +1,2 @@ +./xmlchange RUN_STARTDATE=2001-12-15 +./xmlchange CLM_BLDNML_OPTS=-ignore_warnings --append diff --git a/cime_config/testdefs/testmods_dirs/clm/oldCropCals/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/oldCropCals/user_nl_clm new file mode 100644 index 0000000000..5885f6f2a0 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/oldCropCals/user_nl_clm @@ -0,0 +1,3 @@ +cropcals_rx = .false. +cropcals_rx_adapt = .false. +stream_gdd20_seasons = .false. \ No newline at end of file diff --git a/cime_config/testdefs/testmods_dirs/clm/sowingWindows/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/sowingWindows/user_nl_clm deleted file mode 100644 index 03165bb306..0000000000 --- a/cime_config/testdefs/testmods_dirs/clm/sowingWindows/user_nl_clm +++ /dev/null @@ -1,5 +0,0 @@ -stream_fldFileName_swindow_start = '$DIN_LOC_ROOT/lnd/clm2/cropdata/calendars/processed/swindow_starts_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.nc' -stream_fldFileName_swindow_end = '$DIN_LOC_ROOT/lnd/clm2/cropdata/calendars/processed/swindow_ends_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.nc' -stream_meshfile_cropcal = '$DIN_LOC_ROOT/share/meshes/360x720_120830_ESMFmesh_c20210507_cdf5.nc' -stream_year_first_cropcal = 2000 -stream_year_last_cropcal = 2000 diff --git a/doc/ChangeLog b/doc/ChangeLog index 506b3d5ad0..79d640242d 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,141 @@ =============================================================== +Tag name: ctsm5.2.016 +Originator(s): samrabin (Sam Rabin, UCAR/TSS, samrabin@ucar.edu) +Date: Sat 27 Jul 2024 05:13:08 PM MDT +One-line Summary: Enable new crop calendars for clm60 compsets + +Purpose and description of changes +---------------------------------- + +This commit switches clm60 compsets (really, any compset other than clm45, clm50, and clm51) to use: +- Per-gridcell and -crop sowing windows derived from the GGCMI phase 3b group II static growing seasons +- Per-gridcell and -crop maturity requirements derived from those same growing seasons, over the 1980-2009 growing seasons +- Code to adjust those prescribed maturity requirements based on recent climate + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + +[X] clm6_0 + +[ ] clm5_1 + +[ ] clm5_0 + +[?] ctsm5_0-nwp: Probably? Haven't tested, but didn't exclude that setup. + +[ ] clm4_5 + + +Bugs fixed +---------- + +List of CTSM issues fixed (include CTSM Issue # and description): +- Fixes ESCOMP/CTSM#2584: Reset accumulators to initval instead of 0 (https://github.com/ESCOMP/CTSM/issues/2584) +- Fixes ESCOMP/CTSM#2667: Change out a f19 resolution test on Izumi for f10 (https://github.com/ESCOMP/CTSM/issues/2667) + +Notes of particular relevance for users +--------------------------------------- + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): +- Add cropcals_rx, default false. +- Add cropcals_rx_adapt, default true. This is what changes default behavior. +- Rename stream_year_first_cropcal to stream_year_first_cropcal_swindows +- Rename stream_year_last_cropcal to stream_year_last_cropcal_swindows +- Rename model_year_align_cropcal to model_year_align_cropcal_swindows +- Add stream_year_first_cropcal_cultivar_gdds +- Add stream_year_last_cropcal_cultivar_gdds +- Add model_year_align_cropcal_cultivar_gdds +- Add stream_fldFileName_gdd20_baseline +- Add stream_gdd20_seasons (experimental) +- Add flush_gdd20 (experimental) +- Add stream_fldFileName_gdd20_season_start (experimental) +- Add stream_fldFileName_gdd20_season_end (experimental) +- Add allow_invalid_gdd20_season_inputs + +Changes to the datasets (e.g., parameter, surface or initial files): +- Many new default and optional stream files added. + +Changes to documentation: +- None yet. + +Notes of particular relevance for developers: +--------------------------------------------- + +Caveats for developers: +- tools/contrib/tweak_latlons.py: This is a script I made to avoid the "ambiguous nearest neighbor" issue described in the discussion at https://github.com/orgs/esmf-org/discussions/261. It's hopefully just a temporary thing as we wait for ESMF to fix that bug. If they end up deciding not to, then this should be moved into python/ctsm/ and fully tested. +- Unit and system testing needs to be added for generate_gdd20_baseline.py. +- RXCROPMATURITY SystemTest should be updated to also call generate_gdd20_baseline.py. + +Changes to tests or testing: +- Adds testmods: + - midDecStart: Like decStart, except starts Dec. 15 instead of 30. Useful for ERP/ERS tests where you want the split in December. + - GddGen + - oldCropCals: To recreate previous crop calendars + - RxCropCalsAdaptGGCMI (experimental) + - RxCropCalsNoAdapt + - StreamGDD20Seasons (experimental) +- Removes sowingWindows testmod +- Adds and changes various tests related to crop calendars +- Changes one Izumi test from f19 to 10x15 resolution: + - Was: ERS_D.f19_g17.I1850Clm50BgcCrop.izumi_nag.clm-ciso_monthly_matrixcn_spinup + - Now: ERS_D.f10_f10_mg37.I1850Clm50BgcCrop.izumi_nag.clm-ciso_monthly_matrixcn_spinup + +Note that testmods marked "experimental" are fine to be marked as expected failures, if needed. + + +Testing summary: +---------------- + + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + derecho - FAIL + 1x1_cidadinhoBR tests only. That's fine. Filed Issue #2666: 1x1_cidadinhoBR tests fail in build-namelist_test.pl (https://github.com/ESCOMP/CTSM/issues/2666). + + python testing (if python code has changed; see instructions in python/README.md; document testing done): + + derecho - PASS + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + derecho ----- OK + izumi ------- OK + +Answer changes +-------------- + +Changes answers relative to baseline: + + Summarize any changes to answers, i.e., + - what code configurations: crop-enabled compsets with clm60 (or more specifically, other than clm45, clm50, clm51) + - what platforms/compilers: all + - nature of change (roundoff; larger than roundoff/same climate; new climate): new climate + + Probably? At least, there can be large regional shifts in crop growing season. Results are similar to those for the Prescribed Calendars setup in Rabin et al. (2023; https://gmd.copernicus.org/articles/16/7253/2023/gmd-16-7253-2023.html). + + If this tag changes climate describe the run(s) done to evaluate the new + climate (put details of the simulations in the experiment database) + - None + + +Other details +------------- + +Pull Requests that document the changes (include PR ids): +- ESCOMP/CTSM#2560: Allow prescribed gddmaturity to vary based on recent climate (initial work, v2) (https://github.com/ESCOMP/CTSM/pull/2560) +- ESCOMP/CTSM#2585: Allow "manual" resets of all accumulator types (https://github.com/ESCOMP/CTSM/pull/2585) +- ESCOMP/CTSM#2593: Optionally base GDD20 on per-gridcell windows, not per-hemisphere (https://github.com/ESCOMP/CTSM/pull/2593) +- ESCOMP/CTSM#2661: Merge ctsm5.2.015 into scale-maturity-reqs (https://github.com/ESCOMP/CTSM/pull/2661) +- ESCOMP/CTSM#2664: Enable new crop calendars for clm60 compsets (https://github.com/ESCOMP/CTSM/pull/2664) +- ESCOMP/CTSM#2665: Merge new crop calendars into master (https://github.com/ESCOMP/CTSM/pull/2665) + +=============================================================== +=============================================================== Tag name: ctsm5.2.015 Originator(s): multiple (Samuel Levis,UCAR/TSS,303-665-1310, @mvertens, @jedwards4b, @billsacks, @Katetc) Date: Mon 22 Jul 2024 12:46:17 PM MDT diff --git a/doc/ChangeSum b/doc/ChangeSum index 9c6524b33d..382e5cf7cb 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.2.016 samrabin 07/27/2024 Enable new crop calendars for clm60 compsets ctsm5.2.015 multiple 07/22/2024 Update submodule tags to pass runoff from cism to rof ctsm5.2.014 multiple 07/19/2024 use_matrixcn, use_soil_matrixcn come in as default .false. ctsm5.2.013 glemieux 07/18/2024 FATES Land Use V2 diff --git a/doc/source/users_guide/running-special-cases/Running-with-custom-crop-calendars.rst b/doc/source/users_guide/running-special-cases/Running-with-custom-crop-calendars.rst index b19ebcdb09..878cc0d353 100644 --- a/doc/source/users_guide/running-special-cases/Running-with-custom-crop-calendars.rst +++ b/doc/source/users_guide/running-special-cases/Running-with-custom-crop-calendars.rst @@ -28,8 +28,8 @@ However, geographically- and temporally-varying maps can also be used to prescri stream_meshfile_cropcal = '/glade/p/cesmdata/cseg/inputdata/share/meshes/360x720_120830_ESMFmesh_c20210507_cdf5.nc' ! First and last years on the sowing window datasets - stream_year_first_cropcal = 2000 - stream_year_last_cropcal = 2000 + stream_year_first_cropcal_swindows = 2000 + stream_year_last_cropcal_swindows = 2000 Sowing date ----------- @@ -39,7 +39,7 @@ Specific sowing *dates* can be prescribed for any crop in any gridcell by settin Maturity requirements --------------------- -The heat unit accumulation required for a crop to reach maturity (and thus be harvested) is typically determined by a formula with crop-specific parameters that are specified on the parameter file. However, geographically- and temporally-varying maps of maturity requirement (in units of degree-days) can also be specified using the ``user_nl_clm`` input variable ``stream_fldFileName_cultivar_gdds``. (Note that ``stream_meshfile_cropcal``, ``stream_year_first_cropcal``, and ``stream_year_last_cropcal``---see above---are all also required.) +The heat unit accumulation required for a crop to reach maturity (and thus be harvested) is typically determined by a formula with crop-specific parameters that are specified on the parameter file. However, geographically- and temporally-varying maps of maturity requirement (in units of degree-days) can also be specified using the ``user_nl_clm`` input variable ``stream_fldFileName_cultivar_gdds``. (Note that ``stream_meshfile_cropcal``, ``stream_year_first_cropcal_cultivar_gdds``, and ``stream_year_last_cropcal_cultivar_gdds``---see above---are all also required.) Generating maturity requirements -------------------------------- @@ -52,8 +52,8 @@ In a GDD-generating run, crops are planted on the specified sowing dates and are stream_fldFileName_swindow_start = '/path/to/sowing_date_file.nc' stream_fldFileName_swindow_end = '/path/to/sowing_date_file.nc' stream_meshfile_cropcal = '/path/to/mesh_file.nc' - stream_year_first_cropcal = YEAR - stream_year_last_cropcal = YEAR + stream_year_first_cropcal_swindows = YEAR + stream_year_last_cropcal_swindows = YEAR ! Special settings for "GDD-generating" run generate_crop_gdds = .true. diff --git a/python/ctsm/crop_calendars/check_rxboth_run.py b/python/ctsm/crop_calendars/check_rxboth_run.py index a1014b5e66..fa4affd220 100644 --- a/python/ctsm/crop_calendars/check_rxboth_run.py +++ b/python/ctsm/crop_calendars/check_rxboth_run.py @@ -156,7 +156,16 @@ def main(argv): any_bad = any_bad or gdds_not_obeyed if any_bad: - raise RuntimeError("Unexpected behavior in rxboth run") + msg = "\n ".join( + [ + "Unexpected behavior in rxboth run:", + f"any_bad_import_output: {any_bad_import_output}", + f"any_bad_check_const_vars: {any_bad_check_const_vars}", + f"sdate_not_obeyed: {sdate_not_obeyed}", + f"gdds_not_obeyed: {gdds_not_obeyed}", + ] + ) + raise RuntimeError(msg) if __name__ == "__main__": diff --git a/python/ctsm/crop_calendars/cropcal_module.py b/python/ctsm/crop_calendars/cropcal_module.py index b87d26816f..719d352665 100644 --- a/python/ctsm/crop_calendars/cropcal_module.py +++ b/python/ctsm/crop_calendars/cropcal_module.py @@ -345,7 +345,7 @@ def import_output( my_vars, year_1=None, year_n=None, - my_vegtypes=utils.define_mgdcrop_list(), + my_vegtypes=utils.define_mgdcrop_list_withgrasses(), sdates_rx_ds=None, gdds_rx_ds=None, verbose=False, @@ -425,7 +425,7 @@ def import_output( # Check that e.g., GDDACCUM <= HUI for var_list in [["GDDACCUM", "HUI"], ["SYEARS", "HYEARS"]]: if all(v in this_ds_gs for v in var_list): - any_bad = check_v0_le_v1( + any_bad = any_bad or check_v0_le_v1( this_ds_gs, var_list, both_nan_ok=True, throw_error=throw_errors ) diff --git a/python/ctsm/crop_calendars/cropcal_utils.py b/python/ctsm/crop_calendars/cropcal_utils.py index 00ed2413d2..584046edee 100644 --- a/python/ctsm/crop_calendars/cropcal_utils.py +++ b/python/ctsm/crop_calendars/cropcal_utils.py @@ -207,7 +207,24 @@ def is_each_vegtype(this_vegtypelist, this_filter, this_method): return [is_this_vegtype(x, this_filter, this_method) for x in this_vegtypelist] -def define_mgdcrop_list(): +def define_crop_list(): + """ + List (strings) of managed crops in CLM. + """ + notcrop_list = [ + "tree", + "c3_arctic_grass", + "c3_non-arctic_grass", + "c4_grass", + "shrub", + "not_vegetated", + ] + defined_pftlist = define_pftlist() + is_crop = is_each_vegtype(defined_pftlist, notcrop_list, "notok_contains") + return [defined_pftlist[i] for i, x in enumerate(is_crop) if x] + + +def define_mgdcrop_list_nograsses(): """ List (strings) of managed crops in CLM. """ @@ -217,6 +234,24 @@ def define_mgdcrop_list(): return [defined_pftlist[i] for i, x in enumerate(is_crop) if x] +def define_mgdcrop_list_withgrasses(): + """ + List (strings) of managed crops in CLM. + """ + notcrop_list = [ + "tree", + "c3_arctic_grass", + "c3_non-arctic_grass", + "c4_grass", + "shrub", + "unmanaged", + "not_vegetated", + ] + defined_pftlist = define_pftlist() + is_crop = is_each_vegtype(defined_pftlist, notcrop_list, "notok_contains") + return [defined_pftlist[i] for i, x in enumerate(is_crop) if x] + + def vegtype_str2int(vegtype_str, vegtype_mainlist=None): """ Convert list of vegtype strings to integer index equivalents. diff --git a/python/ctsm/crop_calendars/generate_gdd20_baseline.py b/python/ctsm/crop_calendars/generate_gdd20_baseline.py new file mode 100644 index 0000000000..13668fc850 --- /dev/null +++ b/python/ctsm/crop_calendars/generate_gdd20_baseline.py @@ -0,0 +1,330 @@ +""" +Generate stream_fldFileName_gdd20_baseline file from CTSM outputs +""" + +import sys +import argparse +import os +import datetime as dt +import numpy as np +import xarray as xr +import cftime + +# -- add python/ctsm to path (needed if we want to run generate_gdd20_baseline stand-alone) +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir) +sys.path.insert(1, _CTSM_PYTHON) + +# pylint: disable=wrong-import-position +from ctsm.crop_calendars.import_ds import import_ds +import ctsm.crop_calendars.cropcal_utils as utils +from ctsm.crop_calendars.grid_one_variable import grid_one_variable +from ctsm.crop_calendars.cropcal_module import MISSING_RX_GDD_VAL + +GRIDDING_VAR_LIST = ["patches1d_ixy", "patches1d_jxy", "lat", "lon"] +STREAM_YEAR = 2000 # The year specified for stream_yearFirst and stream_yearLast in the call of +# shr_strdata_init_from_inline() for sdat_cropcal_gdd20_baseline +MGDCROP_LIST = utils.define_crop_list() + + +def _parse_args(): + """ + Set up and parse input arguments + """ + parser = argparse.ArgumentParser( + description=( + "Given a list of CTSM history files, generate stream_fldFileName_gdd20_baseline input" + + "file from the GDD0, GDD8, and GDD10 variables." + ) + ) + + # Required + parser.add_argument( + "-i", + "--input-files", + help="Space-separated string of CTSM history files", + type=str, + required=True, + ) + parser.add_argument( + "-o", + "--output-file", + help="Path to which output file should be saved", + type=str, + required=True, + ) + parser.add_argument( + "-a", + "--author", + help=( + "String to be saved in author attribute of output files." + + "E.g., 'Author Name (authorname@ucar.edu)'" + ), + type=str, + required=True, + ) + + # Optional + parser.add_argument( + "--overwrite", + help="Overwrite existing output file, if any", + action="store_true", + ) + parser.add_argument( + "-y1", + "--first-year", + help=("First calendar year to include"), + type=int, + required=False, + ) + parser.add_argument( + "-yN", + "--last-year", + help=("Last calendar year to include"), + type=int, + required=False, + ) + parser.add_argument( + "-v", + "--variable", + help=("Which type of variable should be processed?"), + required=False, + default="GDDBX", + choices=["GDDBX", "GDDB20"], + ) + + # Get arguments + args = parser.parse_args(sys.argv[1:]) + + # Check arguments + if os.path.exists(args.output_file) and not args.overwrite: + raise FileExistsError("Output file exists but --overwrite is not specified") + + # Get and check input files + args.input_files = args.input_files.split(" ") + for filename in args.input_files: + if not os.path.exists(filename): + raise FileNotFoundError(f"Input file not found: {filename}") + + # Process time slice + # Assumes CESM behavior where data for e.g. 1987 is saved as 1988-01-01. + # It would be more robust, accounting for upcoming behavior (where timestamp for a year is the + # middle of that year), to do slice("YEAR1-01-03", "YEARN-01-02"), but that's not compatible + # with ctsm_pylib as of the version using python 3.7.9. See safer_timeslice() in cropcal_utils. + if args.first_year is not None: + date_1 = f"{args.first_year+1}-01-01" + else: + date_1 = "0000-01-01" + if args.last_year is not None: + date_n = f"{args.last_year+1}-01-01" + else: + date_n = "9999-12-31" + time_slice = slice(date_1, date_n) + + return args, time_slice + + +def _get_cft_list(crop_list): + """ + Given a list of strings, return a list of CFT names that contain any of those strings. + Will include both irrigated and rainfed! + + Args: + crop_list (list): List of crops to look for. + E.g.: ["corn", "cotton"] + + Returns: + cft_str_list: List of CFTs containing any of the crop names in crop_list. + E.g.: ["tropical_corn", "irrigated_tropical_corn", + "temperate_corn", "irrigated_temperate_corn", + "cotton", "irrigated_cotton"] + """ + + cft_str_list = [] + for crop_str in crop_list: + cft_str_list += [x for x in MGDCROP_LIST if crop_str in x] + return cft_str_list + + +def _get_gddn_for_cft(cft_str, variable): + """ + Given a CFT name, return the GDDN variable it uses. + + Args: + cft_str (str): E.g., "irrigated_temperate_corn" + + Returns: + str or None: Name of variable to use (e.g., "GDD8X"). If crop not yet handled, return None. + """ + + gddn = None + gddn_str = None + + gdd0_list_str = ["wheat", "cotton", "rice"] + if cft_str in _get_cft_list(gdd0_list_str): + gddn = 0 + + gdd8_list_str = ["corn", "sugarcane", "miscanthus", "switchgrass"] + if cft_str in _get_cft_list(gdd8_list_str): + gddn = 8 + + gdd10_list_str = ["soybean"] + if cft_str in _get_cft_list(gdd10_list_str): + gddn = 10 + + if gddn is not None: + gddn_str = variable.replace("B", str(gddn)) + + return gddn, gddn_str + + +def _get_output_varname(cft_str): + cft_int = utils.vegtype_str2int(cft_str)[0] + return f"gdd20bl_{cft_int}" + + +def _add_time_axis(da_in): + """ + Adds a size-1 time axis to a DataArray. Needed because CDEPS streams code requires a time axis, + even if the data in question is not supposed to vary over time. + + Args: + da_in (DataArray): xarray DataArray which needs a time axis added + + Returns: + DataArray: da_in with a new 1-step time axis + """ + this_date = np.array(cftime.DatetimeNoLeap(STREAM_YEAR, 1, 1, 0, 0, 0, 0, has_year_zero=True)) + this_date = np.expand_dims(this_date, axis=0) + da_time = xr.DataArray( + data=this_date, + dims={"time": this_date}, + ) + da_out = da_in.expand_dims(time=da_time) + return da_out + + +def setup_output_dataset(input_files, author, variable, year_args, ds_in): + """ + Set up output Dataset + """ + data_var_dict = {} + for gridding_var in GRIDDING_VAR_LIST: + data_var_dict[gridding_var] = ds_in[gridding_var] + ds_out = xr.Dataset( + data_vars=data_var_dict, + attrs={ + "author": author, + "created": dt.datetime.now().astimezone().isoformat(), + "input_year_range": f"{year_args[0]}-{year_args[1]}", + "input_variable": variable, + }, + ) + all_files_in_same_dir = len(np.unique([os.path.dirname(file) for file in input_files])) == 1 + if all_files_in_same_dir: + ds_out.attrs["input_files_dir"] = os.path.dirname(input_files[0]) + ds_out.attrs["input_files"] = ", ".join([os.path.basename(file) for file in input_files]) + else: + ds_out.attrs["input_files"] = ", ".join(input_files) + return ds_out + + +def generate_gdd20_baseline(input_files, output_file, author, time_slice, variable, year_args): + """ + Generate stream_fldFileName_gdd20_baseline file from CTSM outputs + """ + + # Define variables to process + if variable == "GDDBX": + suffix = "X" + elif variable == "GDDB20": + suffix = "20" + else: + raise ValueError(f"-v/--variable {variable} not recoginzed") + var_list_in = [] + for base_temp in [0, 8, 10]: + var_list_in.append(f"GDD{base_temp}{suffix}") + + # Get unique values and sort + input_files = list(set(input_files)) + input_files.sort() + + # Import history files and ensure they have lat/lon dims + ds_in = import_ds(input_files, var_list_in + GRIDDING_VAR_LIST, time_slice=time_slice) + if not all(x in ds_in.dims for x in ["lat", "lon"]): + raise RuntimeError("Input files must have lat and lon dimensions") + + # If needed, find mean over time + if "time" in ds_in.dims: + ds_in = ds_in.mean(dim="time", skipna=True) + + # Set up a dummy DataArray to use for crops without an assigned GDDN variable + dummy_da = xr.DataArray( + data=np.full_like(ds_in[var_list_in[0]].values, MISSING_RX_GDD_VAL), + dims=ds_in[var_list_in[0]].dims, + coords=ds_in[var_list_in[0]].coords, + ) + dummy_da = _add_time_axis(dummy_da) + + # Set up output Dataset + ds_out = setup_output_dataset(input_files, author, variable, year_args, ds_in) + + # Process all crops + encoding_dict = {} + for cft_str in MGDCROP_LIST: + cft_int = utils.vegtype_str2int(cft_str)[0] + print(f"{cft_str} ({cft_int})") + + # Which GDDN history variable does this crop use? E.g., GDD0, GDD10 + gddn, gddn_str = _get_gddn_for_cft(cft_str, variable) + + # Fill any missing values with MISSING_RX_GDD_VAL. This will mean that gddmaturity there + # never changes. + if gddn_str is None: + # Crop not handled yet? It's already filled with missing value + this_da = dummy_da + print(" dummy GDD20") + else: + this_da = ds_in[gddn_str] # Already did ds_in.mean(dim="time") above + this_da = _add_time_axis(this_da) + print(f" {gddn_str}") + this_da = this_da.fillna(MISSING_RX_GDD_VAL) + + # Add attributes of output file + if (gddn is None) != (gddn_str is None): + raise RuntimeError("gddn and gddn_str must either both be None or both be not None") + if gddn_str is None: + long_name = "Dummy GDD20" + else: + long_name = f"GDD{gddn}20" + this_da.attrs["long_name"] = long_name + f" baseline for {cft_str}" + this_da.attrs["units"] = "°C days" + + # Copy that to ds_out + var_out = _get_output_varname(cft_str) + print(f" Output variable {var_out}") + ds_out[var_out] = this_da + encoding_dict[var_out] = {"dtype": "float64"} + + # Grid, if needed + if any(x not in this_da.dims for x in ["lat", "lon"]): + ds_out[var_out] = grid_one_variable(ds_out, var_out) + + # Save + ds_out.to_netcdf(output_file, format="NETCDF3_CLASSIC", encoding=encoding_dict) + + print("Done!") + + +def main(): + """ + main() function for calling generate_gdd20_baseline.py from command line. + """ + args, time_slice = _parse_args() + generate_gdd20_baseline( + args.input_files, + args.output_file, + args.author, + time_slice, + args.variable, + [args.first_year, args.last_year], + ) diff --git a/python/ctsm/crop_calendars/generate_gdds.py b/python/ctsm/crop_calendars/generate_gdds.py index 156ebfb20e..49226e6f75 100644 --- a/python/ctsm/crop_calendars/generate_gdds.py +++ b/python/ctsm/crop_calendars/generate_gdds.py @@ -178,6 +178,7 @@ def main( mxmats, cc.get_gs_len_da, skip_crops, + outdir_figs, logger, ) diff --git a/python/ctsm/crop_calendars/generate_gdds_functions.py b/python/ctsm/crop_calendars/generate_gdds_functions.py index 8af2fdc049..14bd6b2e40 100644 --- a/python/ctsm/crop_calendars/generate_gdds_functions.py +++ b/python/ctsm/crop_calendars/generate_gdds_functions.py @@ -22,6 +22,7 @@ # pylint: disable=import-error from ctsm.crop_calendars.cropcal_figs_module import * from matplotlib.transforms import Bbox + import matplotlib.pyplot as plt warnings.filterwarnings( "ignore", @@ -63,7 +64,7 @@ def error(logger, string): raise RuntimeError(string) -def check_sdates(dates_ds, sdates_rx, logger, verbose=False): +def check_sdates(dates_ds, sdates_rx, outdir_figs, logger, verbose=False): """ Checking that input and output sdates match """ @@ -106,8 +107,28 @@ def check_sdates(dates_ds, sdates_rx, logger, verbose=False): log(logger, out_map_notnan[here][0:4]) log(logger, "diff:") log(logger, diff_map_notnan[here][0:4]) + first_diff = all_ok all_ok = False + if CAN_PLOT: + sdate_diffs_dir = os.path.join(outdir_figs, "sdate_diffs") + if first_diff: + log(logger, f"Saving sdate difference figures to {sdate_diffs_dir}") + if not os.path.exists(sdate_diffs_dir): + os.makedirs(sdate_diffs_dir) + in_map.where(~np.isnan(out_map)).plot() + plt.title(f"{vegtype_str} sdates in (masked)") + plt.savefig(os.path.join(sdate_diffs_dir, f"{vegtype_str}.in.png")) + plt.close() + out_map.plot() + plt.title(f"{vegtype_str} sdates out") + plt.savefig(os.path.join(sdate_diffs_dir, f"{vegtype_str}.out.png")) + plt.close() + diff_map.plot() + plt.title(f"{vegtype_str} sdates diff (out - in)") + plt.savefig(os.path.join(sdate_diffs_dir, f"{vegtype_str}.diff.png")) + plt.close() + if not any_found: error(logger, "No matching variables found in sdates_rx!") @@ -234,6 +255,7 @@ def import_and_process_1yr( mxmats, get_gs_len_da, skip_crops, + outdir_figs, logger, ): """ @@ -260,9 +282,9 @@ def import_and_process_1yr( # Get list of crops to include if skip_crops is not None: - crops_to_read = [c for c in utils.define_mgdcrop_list() if c not in skip_crops] + crops_to_read = [c for c in utils.define_mgdcrop_list_withgrasses() if c not in skip_crops] else: - crops_to_read = utils.define_mgdcrop_list() + crops_to_read = utils.define_mgdcrop_list_withgrasses() print(h1_filelist) dates_ds = import_ds( @@ -272,6 +294,8 @@ def import_and_process_1yr( time_slice=slice(f"{this_year}-01-01", f"{this_year}-12-31"), chunks=chunks, ) + for timestep in dates_ds["time"].values: + print(timestep) if dates_ds.dims["time"] > 1: if dates_ds.dims["time"] == 365: @@ -466,7 +490,7 @@ def import_and_process_1yr( # Import expected sowing dates. This will also be used as our template output file. imported_sdates = isinstance(sdates_rx, str) sdates_rx = import_rx_dates("s", sdates_rx, incl_patches1d_itype_veg, mxsowings, logger) - check_sdates(dates_incl_ds, sdates_rx, logger) + check_sdates(dates_incl_ds, sdates_rx, outdir_figs, logger) # Import hdates, if needed imported_hdates = isinstance(hdates_rx, str) @@ -575,6 +599,7 @@ def import_and_process_1yr( this_crop_gddaccum_da = this_crop_ds[clm_gdd_var] if save_figs: this_crop_gddharv_da = this_crop_ds["GDDHARV"] + check_gddharv = True if not this_crop_gddaccum_da.size: continue log(logger, f" {vegtype_str}...") @@ -625,11 +650,15 @@ def import_and_process_1yr( + "NaN after extracting GDDs accumulated at harvest", ) if save_figs and np.any(np.isnan(gddharv_atharv_p)): - log( - logger, - f" ❗ {np.sum(np.isnan(gddharv_atharv_p))}/{len(gddharv_atharv_p)} " - + "NaN after extracting GDDHARV", - ) + if np.all(np.isnan(gddharv_atharv_p)): + log(logger, " ❗ All GDDHARV are NaN; should only affect figure") + check_gddharv = False + else: + log( + logger, + f" ❗ {np.sum(np.isnan(gddharv_atharv_p))}/{len(gddharv_atharv_p)} " + + "NaN after extracting GDDHARV", + ) # Assign these to growing seasons based on whether gs crossed new year this_year_active_patch_indices = [ @@ -712,9 +741,15 @@ def import_and_process_1yr( ) else: error(logger, "Unexpected NaN for last season's GDD accumulation.") - if save_figs and np.any( - np.isnan( - gddharv_yp_list[var][year_index - 1, active_this_year_where_gs_lastyr_indices] + if ( + save_figs + and check_gddharv + and np.any( + np.isnan( + gddharv_yp_list[var][ + year_index - 1, active_this_year_where_gs_lastyr_indices + ] + ) ) ): if incorrectly_daily: @@ -1160,9 +1195,13 @@ def make_figures( else: error(logger, f"layout {layout} not recognized") - this_min = int(np.round(np.nanmin(gddharv_map_yx))) - this_max = int(np.round(np.nanmax(gddharv_map_yx))) - this_title = f"{run1_name} (range {this_min}–{this_max})" + gddharv_all_nan = np.all(np.isnan(gddharv_map_yx.values)) + if gddharv_all_nan: + this_title = f"{run1_name} (GDDHARV all NaN?)" + else: + this_min = int(np.round(np.nanmin(gddharv_map_yx))) + this_max = int(np.round(np.nanmax(gddharv_map_yx))) + this_title = f"{run1_name} (range {this_min}–{this_max})" make_gengdd_map( this_axis, gddharv_map_yx, @@ -1195,7 +1234,7 @@ def make_figures( ) # Difference - if layout == "3x2": + if not gddharv_all_nan and layout == "3x2": this_axis = fig.add_subplot(spec[2, 0], projection=ccrs.PlateCarree()) this_min = int(np.round(np.nanmin(gdd_map_yx))) this_max = int(np.round(np.nanmax(gdd_map_yx))) diff --git a/python/ctsm/crop_calendars/import_ds.py b/python/ctsm/crop_calendars/import_ds.py index 77a22b626b..486757492f 100644 --- a/python/ctsm/crop_calendars/import_ds.py +++ b/python/ctsm/crop_calendars/import_ds.py @@ -41,6 +41,28 @@ def compute_derived_vars(ds_in, var): return ds_in +def manual_mfdataset(filelist, my_vars, my_vegtypes, time_slice): + """ + Opening a list of files with Xarray's open_mfdataset requires dask. This function is a + workaround for Python environments that don't have dask. + """ + ds_out = None + for filename in filelist: + ds_in = xr.open_dataset(filename) + ds_in = mfdataset_preproc(ds_in, my_vars, my_vegtypes, time_slice) + if ds_out is None: + ds_out = ds_in + else: + ds_out = xr.concat( + [ds_out, ds_in], + data_vars="minimal", + compat="override", + coords="all", + dim="time", + ) + return ds_out + + def mfdataset_preproc(ds_in, vars_to_import, vegtypes_to_import, time_slice): """ Function to drop unwanted variables in preprocessing of open_mfdataset(). @@ -221,22 +243,20 @@ def import_ds( if isinstance(filelist, list): with warnings.catch_warnings(): warnings.filterwarnings(action="ignore", category=DeprecationWarning) - if find_spec("dask") is None: - raise ModuleNotFoundError( - "You have asked xarray to import a list of files as a single Dataset using" - " open_mfdataset(), but this requires dask, which is not available.\nFile" - f" list: {filelist}" - ) - this_ds = xr.open_mfdataset( - sorted(filelist), - data_vars="minimal", - preprocess=mfdataset_preproc_closure, - compat="override", - coords="all", - concat_dim="time", - combine="nested", - chunks=chunks, - ) + dask_unavailable = find_spec("dask") is None + if dask_unavailable: + this_ds = manual_mfdataset(filelist, my_vars, my_vegtypes, time_slice) + else: + this_ds = xr.open_mfdataset( + sorted(filelist), + data_vars="minimal", + preprocess=mfdataset_preproc_closure, + compat="override", + coords="all", + concat_dim="time", + combine="nested", + chunks=chunks, + ) elif isinstance(filelist, str): this_ds = xr.open_dataset(filelist, chunks=chunks) this_ds = mfdataset_preproc(this_ds, my_vars, my_vegtypes, time_slice) diff --git a/python/ctsm/crop_calendars/interpolate_gdds.py b/python/ctsm/crop_calendars/interpolate_gdds.py index 2aa0b79997..123d40af38 100755 --- a/python/ctsm/crop_calendars/interpolate_gdds.py +++ b/python/ctsm/crop_calendars/interpolate_gdds.py @@ -6,6 +6,7 @@ import sys import argparse import logging +import re import xarray as xr # -- add python/ctsm to path (needed if we want to run this stand-alone) @@ -64,6 +65,14 @@ def _setup_process_args(): type=str, required=True, ) + parser.add_argument( + "-p", + "--variable-prefix", + help="Interpolate variables whose names start with this string", + type=str, + required=False, + default="gdd1_", + ) parser.add_argument( "--overwrite", help="If output file exists, overwrite it.", @@ -110,8 +119,12 @@ def interpolate_gdds(args): for var in ds_in: # Check variable - if "gdd1_" not in var: - raise RuntimeError(f"Unexpected variable {var} on input file") + if "lat" not in ds_in[var].dims and "lon" not in ds_in[var].dims: + print(f"Skipping variable {var} with dimensions {ds_in[var].dims}") + continue + if not re.compile("^" + args.variable_prefix).match(var): + print(f"Unexpected variable {var} on input file. Skipping.") + continue if args.dry_run: continue diff --git a/python/ctsm/run_sys_tests.py b/python/ctsm/run_sys_tests.py index f5a7ed04dc..6afc43cd19 100644 --- a/python/ctsm/run_sys_tests.py +++ b/python/ctsm/run_sys_tests.py @@ -736,13 +736,29 @@ def _check_py_env(test_attributes): # whether import is possible. # pylint: disable=import-error disable - # Check requirements for FSURDATMODIFYCTSM, if needed - if any("FSURDATMODIFYCTSM" in t for t in test_attributes): + # Check requirements for using modify_fsurdat Python module, if needed + modify_fsurdat_users = ["FSURDATMODIFYCTSM", "RXCROPMATURITY"] + if any(any(u in t for u in modify_fsurdat_users) for t in test_attributes): try: import ctsm.modify_input_files.modify_fsurdat except ModuleNotFoundError as err: raise ModuleNotFoundError("modify_fsurdat" + err_msg) from err + # Check requirements for RXCROPMATURITY, if needed + if any("RXCROPMATURITY" in t for t in test_attributes): + try: + import ctsm.crop_calendars.check_rxboth_run + except ModuleNotFoundError as err: + raise ModuleNotFoundError("check_rxboth_run" + err_msg) from err + try: + import ctsm.crop_calendars.generate_gdds + except ModuleNotFoundError as err: + raise ModuleNotFoundError("generate_gdds" + err_msg) from err + try: + import ctsm.crop_calendars.interpolate_gdds + except ModuleNotFoundError as err: + raise ModuleNotFoundError("interpolate_gdds" + err_msg) from err + # Check that list for any testmods that use modify_fates_paramfile.py testmods_to_check = ["clm-FatesColdTwoStream", "clm-FatesColdTwoStreamNoCompFixedBioGeo"] testmods = _get_testmod_list(test_attributes) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index be785e745d..47c0295593 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -235,8 +235,16 @@ def create_2d_coords(self): lons_size = self.center_lons.size # -- convert center points from 1d to 2d - self.center_lat2d = np.broadcast_to(self.center_lats[:], (lons_size, lats_size)) - self.center_lon2d = np.broadcast_to(self.center_lons[:], (lons_size, lats_size)) + try: + self.center_lat2d = np.broadcast_to(self.center_lats[:], (lons_size, lats_size)) + self.center_lon2d = np.broadcast_to(self.center_lons[:], (lons_size, lats_size)) + except ValueError: + self.center_lat2d = np.broadcast_to( + np.expand_dims(self.center_lats[:], 0), (lons_size, lats_size) + ) + self.center_lon2d = np.broadcast_to( + np.expand_dims(self.center_lons[:], 1), (lons_size, lats_size) + ) elif self.lat_dims == 2: # -- 2D lats and lons dims = self.center_lons.shape @@ -351,7 +359,7 @@ def calculate_corners(self, unit="degrees"): ) # Longitudes should stay within 0 to 360 if np.any(self.corner_lons > 360.0): - abort("Corners have longitudes greater than 360") + abort(f"Corners have longitudes greater than 360 (max: {np.max(self.corner_lons)})") if np.any(self.corner_lons < 0.0): logger.warning( "Corners have longitudes less than zero -- %s %s", diff --git a/src/biogeochem/CNDVType.F90 b/src/biogeochem/CNDVType.F90 index 19a0f64f7d..fb6b3d9753 100644 --- a/src/biogeochem/CNDVType.F90 +++ b/src/biogeochem/CNDVType.F90 @@ -439,7 +439,7 @@ subroutine UpdateAccVars(this, bounds, t_a10_patch, t_ref2m_patch) use shr_const_mod , only : SHR_CONST_CDAY, SHR_CONST_TKFRZ use clm_time_manager , only : get_step_size, get_nstep, get_curr_date use pftconMod , only : ndllf_dcd_brl_tree - use accumulMod , only : update_accum_field, extract_accum_field, accumResetVal + use accumulMod , only : update_accum_field, extract_accum_field, markreset_accum_field ! ! !ARGUMENTS: class(dgvs_type) , intent(inout) :: this @@ -489,25 +489,34 @@ subroutine UpdateAccVars(this, bounds, t_a10_patch, t_ref2m_patch) ! Accumulate and extract AGDDTW (gdd base twmax, which is 23 deg C ! for boreal woody patches) + ! SSR 2024-06-07: Don't wrap this do-loop in an "if it's not time to reset." Behavior would + ! be identical for now, but if "missed update" behavior is fixed (see ESCOMP/CTSM#2585), you + ! would end up updating AGDDTW with uninitialized values. do p = begp,endp rbufslp(p) = max(0._r8, & (t_a10_patch(p) - SHR_CONST_TKFRZ - dgv_ecophyscon%twmax(ndllf_dcd_brl_tree)) & * dtime/SHR_CONST_CDAY) - if (month==1 .and. day==1 .and. secs==int(dtime)) rbufslp(p) = accumResetVal end do + if (month==1 .and. day==1 .and. secs==int(dtime)) then + ! Reset annually + call markreset_accum_field('AGDDTW') + end if call update_accum_field ('AGDDTW', rbufslp, nstep) call extract_accum_field ('AGDDTW', this%agddtw_patch, nstep) ! Accumulate and extract AGDD + ! SSR 2024-06-07: Don't wrap this do-loop in an "if it's not time to reset." Behavior would + ! be identical for now, but if "missed update" behavior is fixed (see ESCOMP/CTSM#2585), you + ! would end up updating AGDD with uninitialized values. do p = begp,endp rbufslp(p) = max(0.0_r8, & (t_ref2m_patch(p) - (SHR_CONST_TKFRZ + 5.0_r8)) * dtime/SHR_CONST_CDAY) - ! - ! Fix (for bug 1858) from Sam Levis to reset the annual AGDD variable - ! - if (month==1 .and. day==1 .and. secs==int(dtime)) rbufslp(p) = accumResetVal end do + if (month==1 .and. day==1 .and. secs==int(dtime)) then + ! Reset annually + call markreset_accum_field('AGDD') + end if call update_accum_field ('AGDD', rbufslp, nstep) call extract_accum_field ('AGDD', this%agdd_patch, nstep) diff --git a/src/biogeochem/CNPhenologyMod.F90 b/src/biogeochem/CNPhenologyMod.F90 index 9c7e5a331a..2ce7f029da 100644 --- a/src/biogeochem/CNPhenologyMod.F90 +++ b/src/biogeochem/CNPhenologyMod.F90 @@ -149,6 +149,9 @@ module CNPhenologyMod logical, public :: generate_crop_gdds = .false. ! If true, harvest the day before next sowing logical, public :: use_mxmat = .true. ! If true, ignore crop maximum growing season length + ! For use with adapt_cropcal_rx_cultivar_gdds .true. + real(r8), parameter :: min_gdd20_baseline = 0._r8 ! If gdd20_baseline_patch is ≤ this, do not consider baseline. + ! Constants for seasonal decidious leaf onset and offset logical, private :: onset_thresh_depends_on_veg = .false. ! If onset threshold depends on vegetation type integer, public, parameter :: critical_daylight_constant = 1 @@ -2024,7 +2027,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_swindows, use_cropcal_rx_cultivar_gdds, use_cropcal_streams + use clm_varctl , only : use_cropcal_rx_swindows ! ! !ARGUMENTS: integer , intent(in) :: num_pcropp ! number of prog crop patches in filter @@ -2189,7 +2192,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & cnveg_state_inst%gddmaturity_thisyr(p,s) = -1._r8 crop_inst%gddaccum_thisyr_patch(p,s) = -1._r8 crop_inst%hui_thisyr_patch(p,s) = -1._r8 - crop_inst%sowing_reason_perharv_patch = -1._r8 + crop_inst%sowing_reason_perharv_patch(p,s) = -1._r8 crop_inst%harvest_reason_thisyr_patch(p,s) = -1._r8 do k = repr_grain_min, repr_grain_max cnveg_carbonflux_inst%repr_grainc_to_food_perharv_patch(p,s,k) = 0._r8 @@ -2458,7 +2461,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , & do_harvest = .true. fake_harvest = .true. harvest_reason = HARVEST_REASON_SOWNBADDEC31 - else if (use_cropcal_streams .and. do_plant .and. .not. did_plant) then + else if (use_cropcal_rx_swindows .and. do_plant .and. .not. did_plant) then ! Today was supposed to be the planting day, but the previous crop still hasn't been harvested. do_harvest = .true. harvest_reason = HARVEST_REASON_SOWTODAY @@ -2775,7 +2778,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_cultivar_gdds, use_cropcal_streams + use clm_varctl , only : use_cropcal_rx_cultivar_gdds, adapt_cropcal_rx_cultivar_gdds use clm_varcon , only : c13ratio, c14ratio use clm_varpar , only : mxsowings use pftconMod , only : ntmp_corn, nswheat, nwwheat, ntmp_soybean @@ -2806,6 +2809,7 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & ! LOCAL VARAIBLES: integer s ! growing season index integer k ! grain pool index + real(r8) gdd20 ! GDD*20 value for this crop type real(r8) gdd_target ! cultivar GDD target this growing season real(r8) this_sowing_reason ! number representing sowing reason(s) logical did_rx_gdds ! did this patch use a prescribed harvest requirement? @@ -2884,42 +2888,71 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & endif endif + ! which GDD*20 variable does this crop use? + if (ivt(p) == ntmp_soybean .or. ivt(p) == nirrig_tmp_soybean .or. & + ivt(p) == ntrp_soybean .or. ivt(p) == nirrig_trp_soybean) then + gdd20 = gdd1020(p) + else if (ivt(p) == ntmp_corn .or. ivt(p) == nirrig_tmp_corn .or. & + ivt(p) == ntrp_corn .or. ivt(p) == nirrig_trp_corn .or. & + ivt(p) == nsugarcane .or. ivt(p) == nirrig_sugarcane .or. & + ivt(p) == nmiscanthus .or. ivt(p) == nirrig_miscanthus .or. & + ivt(p) == nswitchgrass .or. ivt(p) == nirrig_switchgrass) then + gdd20 = gdd820(p) + else if (ivt(p) == nswheat .or. ivt(p) == nirrig_swheat .or. & + ivt(p) == nwwheat .or. ivt(p) == nirrig_wwheat .or. & + ivt(p) == ncotton .or. ivt(p) == nirrig_cotton .or. & + ivt(p) == nrice .or. ivt(p) == nirrig_rice) then + gdd20 = gdd020(p) + else + write(iulog, *) 'ERROR: PlantCrop(): unrecognized ivt for gdd20: ',ivt(p) + call endrun(msg="Stopping") + end if + ! set GDD target did_rx_gdds = .false. if (use_cropcal_rx_cultivar_gdds .and. crop_inst%rx_cultivar_gdds_thisyr_patch(p,sowing_count(p)) .ge. 0._r8) then gddmaturity(p) = crop_inst%rx_cultivar_gdds_thisyr_patch(p,sowing_count(p)) did_rx_gdds = .true. + if (adapt_cropcal_rx_cultivar_gdds .and. crop_inst%gdd20_baseline_patch(p) > min_gdd20_baseline) then + gddmaturity(p) = gddmaturity(p) * gdd20 / crop_inst%gdd20_baseline_patch(p) + !TODO Sam Rabin: Set maximum and minimum gddmaturity + end if else if (ivt(p) == nwwheat .or. ivt(p) == nirrig_wwheat) then gddmaturity(p) = hybgdd(ivt(p)) else if (ivt(p) == ntmp_soybean .or. ivt(p) == nirrig_tmp_soybean .or. & - ivt(p) == ntrp_soybean .or. ivt(p) == nirrig_trp_soybean) then - gddmaturity(p) = min(gdd1020(p), hybgdd(ivt(p))) - end if - if (ivt(p) == ntmp_corn .or. ivt(p) == nirrig_tmp_corn .or. & + ivt(p) == ntrp_soybean .or. ivt(p) == nirrig_trp_soybean .or. & + ivt(p) == nswheat .or. ivt(p) == nirrig_swheat .or. & + ivt(p) == ncotton .or. ivt(p) == nirrig_cotton .or. & + ivt(p) == nrice .or. ivt(p) == nirrig_rice) then + gddmaturity(p) = min(gdd20, hybgdd(ivt(p))) + else if (ivt(p) == ntmp_corn .or. ivt(p) == nirrig_tmp_corn .or. & ivt(p) == ntrp_corn .or. ivt(p) == nirrig_trp_corn .or. & ivt(p) == nsugarcane .or. ivt(p) == nirrig_sugarcane .or. & ivt(p) == nmiscanthus .or. ivt(p) == nirrig_miscanthus .or. & ivt(p) == nswitchgrass .or. ivt(p) == nirrig_switchgrass) then - gddmaturity(p) = max(950._r8, min(gdd820(p)*0.85_r8, hybgdd(ivt(p)))) + gddmaturity(p) = max(950._r8, min(gdd20*0.85_r8, hybgdd(ivt(p)))) if (do_plant_normal) then gddmaturity(p) = max(950._r8, min(gddmaturity(p)+150._r8, 1850._r8)) end if - end if - if (ivt(p) == nswheat .or. ivt(p) == nirrig_swheat .or. & - ivt(p) == ncotton .or. ivt(p) == nirrig_cotton .or. & - ivt(p) == nrice .or. ivt(p) == nirrig_rice) then - gddmaturity(p) = min(gdd020(p), hybgdd(ivt(p))) + else + write(iulog, *) 'ERROR: PlantCrop(): unrecognized ivt for GDD target: ',ivt(p) + call endrun(msg="Stopping") end if endif - if (use_cropcal_streams .and. gddmaturity(p) < min_gddmaturity) then - if (did_rx_gdds) then - write(iulog,*) 'Some patch with ivt ',ivt(p),' has rx gddmaturity',gddmaturity(p),'; using min_gddmaturity instead (',min_gddmaturity,')' - endif - gddmaturity(p) = min_gddmaturity - endif + if (gddmaturity(p) < min_gddmaturity) then + if (use_cropcal_rx_cultivar_gdds .or. generate_crop_gdds) then + if (did_rx_gdds) then + write(iulog,*) 'Some patch with ivt ',ivt(p),' has rx gddmaturity ',gddmaturity(p),'; using min_gddmaturity instead (',min_gddmaturity,')' + end if + gddmaturity(p) = min_gddmaturity + else + write(iulog, *) 'ERROR: PlantCrop(): gddmaturity < minimum for ivt ',ivt(p) + call endrun(msg="Stopping") + end if + end if ! Initialize allocation coefficients. ! Because crops have no live carbon pools when planted but not emerged, this shouldn't diff --git a/src/biogeochem/CropType.F90 b/src/biogeochem/CropType.F90 index 760aa21b76..0f650a4a9f 100644 --- a/src/biogeochem/CropType.F90 +++ b/src/biogeochem/CropType.F90 @@ -52,6 +52,9 @@ module CropType 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 :: gdd20_baseline_patch (:) ! GDD20 baseline for this patch (ddays) [patch] + real(r8), pointer :: gdd20_season_start_patch(:) ! gdd20 season start date for this patch (day of year) [patch]. Real to enable history field. + real(r8), pointer :: gdd20_season_end_patch (:) ! gdd20 season end date for this patch (day of year) [patch]. Real to enable history field. 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] @@ -235,6 +238,9 @@ subroutine InitAllocate(this, bounds) 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%gdd20_baseline_patch(begp:endp)) ; this%gdd20_baseline_patch(:) = spval + allocate(this%gdd20_season_start_patch(begp:endp)); this%gdd20_season_start_patch(:) = spval + allocate(this%gdd20_season_end_patch(begp:endp)) ; this%gdd20_season_end_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 @@ -358,6 +364,19 @@ subroutine InitHistory(this, bounds) avgflag='I', long_name='Reason for each crop harvest; should only be output annually', & ptr_patch=this%harvest_reason_thisyr_patch, default='inactive') + this%gdd20_baseline_patch(begp:endp) = spval + call hist_addfld1d (fname='GDD20_BASELINE', units='ddays', & + avgflag='I', long_name='Baseline mean growing-degree days accumulated during accumulation period (from input)', & + ptr_patch=this%gdd20_baseline_patch, default='inactive') + this%gdd20_season_start_patch(begp:endp) = spval + call hist_addfld1d (fname='GDD20_SEASON_START', units='day of year', & + avgflag='I', long_name='Start of the GDD20 accumulation season (from input)', & + ptr_patch=this%gdd20_season_start_patch, default='inactive') + this%gdd20_season_end_patch(begp:endp) = spval + call hist_addfld1d (fname='GDD20_SEASON_END', units='day of year', & + avgflag='I', long_name='End of the GDD20 accumulation season (from input)', & + ptr_patch=this%gdd20_season_end_patch, default='inactive') + end subroutine InitHistory subroutine InitCold(this, bounds) @@ -772,7 +791,7 @@ subroutine CropUpdateAccVars(this, bounds, t_ref2m_patch, t_soisno_col) ! Should only be called if use_crop is true. ! ! !USES: - use accumulMod , only : update_accum_field, extract_accum_field, accumResetVal + use accumulMod , only : update_accum_field, extract_accum_field, markreset_accum_field use shr_const_mod , only : SHR_CONST_CDAY, SHR_CONST_TKFRZ use clm_time_manager , only : get_step_size, get_nstep use clm_varpar , only : nlevsno, nlevmaxurbgrnd @@ -848,7 +867,8 @@ subroutine CropUpdateAccVars(this, bounds, t_ref2m_patch, t_soisno_col) rbufslp(p) = rbufslp(p) * this%vf_patch(p) end if else - rbufslp(p) = accumResetVal + call markreset_accum_field('HUI', p) + call markreset_accum_field('GDDACCUM', p) end if end do call update_accum_field ('HUI', rbufslp, nstep) @@ -872,7 +892,7 @@ subroutine CropUpdateAccVars(this, bounds, t_ref2m_patch, t_soisno_col) rbufslp(p) = rbufslp(p) * this%vf_patch(p) end if else - rbufslp(p) = accumResetVal + call markreset_accum_field('GDDTSOI', p) end if end do call update_accum_field ('GDDTSOI', rbufslp, nstep) diff --git a/src/biogeophys/EnergyFluxType.F90 b/src/biogeophys/EnergyFluxType.F90 index 2e709596a1..16929d9708 100644 --- a/src/biogeophys/EnergyFluxType.F90 +++ b/src/biogeophys/EnergyFluxType.F90 @@ -988,7 +988,7 @@ subroutine UpdateAccVars (this, bounds) ! ! USES use clm_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 accumulMod , only : update_accum_field, extract_accum_field use abortutils , only : endrun ! ! !ARGUMENTS: diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index ab310650c8..707218cc27 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -8,6 +8,7 @@ module TemperatureType use decompMod , only : bounds_type use abortutils , only : endrun use clm_varctl , only : use_cndv, iulog, use_luna, use_crop, use_biomass_heat_storage + use clm_varctl , only : flush_gdd20 use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevurb, nlevmaxurbgrnd use clm_varcon , only : spval, ispval use GridcellType , only : grc @@ -129,6 +130,7 @@ module TemperatureType procedure, public :: InitAccBuffer procedure, public :: InitAccVars procedure, public :: UpdateAccVars + procedure, private :: UpdateAccVars_CropGDDs end type temperature_type @@ -596,9 +598,7 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) call hist_addfld1d (fname='GDD0', units='ddays', & avgflag='A', long_name='Growing degree days base 0C from planting', & ptr_patch=this%gdd0_patch, default='inactive') - end if - if (use_crop) then this%gdd8_patch(begp:endp) = spval call hist_addfld1d (fname='GDD8', units='ddays', & avgflag='A', long_name='Growing degree days base 8C from planting', & @@ -609,6 +609,21 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) avgflag='A', long_name='Growing degree days base 10C from planting', & ptr_patch=this%gdd10_patch, default='inactive') + this%gdd0_patch(begp:endp) = spval + call hist_addfld1d (fname='GDD0X', units='ddays', & + avgflag='X', long_name='Growing degree days base 0C from planting, max', & + ptr_patch=this%gdd0_patch, default='inactive') + + this%gdd8_patch(begp:endp) = spval + call hist_addfld1d (fname='GDD8X', units='ddays', & + avgflag='X', long_name='Growing degree days base 8C from planting, max', & + ptr_patch=this%gdd8_patch, default='inactive') + + this%gdd10_patch(begp:endp) = spval + call hist_addfld1d (fname='GDD10X', units='ddays', & + avgflag='X', long_name='Growing degree days base 10C from planting, max', & + ptr_patch=this%gdd10_patch, default='inactive') + this%gdd020_patch(begp:endp) = spval call hist_addfld1d (fname='GDD020', units='ddays', & avgflag='A', long_name='Twenty year average of growing degree days base 0C from planting', & @@ -895,6 +910,7 @@ subroutine Restart(this, bounds, ncid, flag, is_simple_buildtemp, is_prog_buildt ! !LOCAL VARIABLES: integer :: j,c ! indices logical :: readvar ! determine if variable is on initial file + integer :: idata !----------------------------------------------------------------------- call restartvar(ncid=ncid, flag=flag, varname='T_SOISNO', xtype=ncd_double, & @@ -1357,22 +1373,133 @@ subroutine InitAccVars(this, bounds) end subroutine InitAccVars - !----------------------------------------------------------------------- - subroutine UpdateAccVars (this, bounds) + subroutine UpdateAccVars_CropGDDs(this, rbufslp, begp, endp, month, day, secs, dtime, nstep, basetemp_int, gddx_patch, crop_inst) ! ! USES use shr_const_mod , only : SHR_CONST_CDAY, SHR_CONST_TKFRZ + use accumulMod , only : update_accum_field, extract_accum_field, markreset_accum_field + use clm_time_manager , only : is_doy_in_interval, get_curr_calday + use pftconMod , only : npcropmin + use CropType, only : crop_type + ! + ! !ARGUMENTS + class(temperature_type) :: this + real(r8), intent(inout), pointer, dimension(:) :: rbufslp ! temporary single level - pft level + integer, intent(in) :: begp, endp + integer, intent(in) :: month, day, secs, dtime, nstep + integer, intent(in) :: basetemp_int ! Crop base temperature. Integer to avoid possible float weirdness + real(r8), intent(inout), pointer, dimension(:) :: gddx_patch ! E.g., gdd0_patch + type(crop_type), intent(inout) :: crop_inst + ! + ! !LOCAL VARIABLES + real(r8) :: basetemp_r8 ! real(r8) version of basetemp for arithmetic + real(r8) :: max_accum ! Maximum daily accumulation + character(8) :: field_name ! E.g., GDD0 + character(32) :: format_string + integer :: p + logical :: in_accumulation_season + real(r8) :: lat ! latitude + integer :: gdd20_season_start, gdd20_season_end + integer :: jday ! Julian day of year (1, ..., 366) + logical :: stream_gdd20_seasons_tt ! Local derivation of this to avoid circular dependency + + associate( & + gdd20_season_starts => crop_inst%gdd20_season_start_patch, & + gdd20_season_ends => crop_inst%gdd20_season_end_patch & + ) + + basetemp_r8 = real(basetemp_int, r8) + + ! SSR 2024-06-13: This should probably be _prev_. Keeping it _curr_ for now for consistency with + ! parent subroutine UpdateAccVars(), which uses get_curr_date() to get the month/day/etc. values + ! that are passed into this subroutine. + jday = int(get_curr_calday()) + + ! Get maximum daily accumulation + if (basetemp_int == 0) then + ! SSR 2024-05-31: I'm not sure why this was different for base temp 0, but I'm keeping it as I refactor into UpdateAccVars_CropGDDs() + max_accum = 26._r8 + else + max_accum = 30._r8 + end if + + ! Get field name + if (basetemp_int < 10) then + format_string = "(A3,I1)" + else if (basetemp_int < 100) then + format_string = "(A3,I2)" + else + format_string = "(A3,I3)" + end if + write(field_name, format_string) "GDD",basetemp_int + + stream_gdd20_seasons_tt = any(gdd20_season_starts(begp:endp) > 0.5_r8) .and. any(gdd20_season_starts(begp:endp) < 366.5_r8) + + do p = begp,endp + + ! Avoid unnecessary calculations over inactive points + if (.not. patch%active(p)) then + cycle + end if + + ! Is this patch in its gdd20 accumulation season? + ! First, check based on latitude. This will be fallback if read-in gdd20 accumulation season is invalid. + lat = grc%latdeg(patch%gridcell(p)) + in_accumulation_season = & + ((month > 3 .and. month < 10) .and. lat >= 0._r8) .or. & + ((month > 9 .or. month < 4) .and. lat < 0._r8) + ! Replace with read-in gdd20 accumulation season, if needed and valid + ! (If these aren't being read in or they're invalid, they'll be -1) + if (stream_gdd20_seasons_tt .and. patch%itype(p) >= npcropmin) then + gdd20_season_start = int(gdd20_season_starts(p)) + gdd20_season_end = int(gdd20_season_ends(p)) + if (gdd20_season_start >= 1 .and. gdd20_season_end >= 1) then + if (gdd20_season_start > 366 .or. gdd20_season_end > 366) then + write(iulog,*) 'invalid gdd20 season!' + write(iulog,*) ' start: ',gdd20_season_start + write(iulog,*) ' end: ',gdd20_season_end + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + in_accumulation_season = is_doy_in_interval( & + gdd20_season_start, gdd20_season_end, jday) + end if + end if + + if (month==1 .and. day==1 .and. secs==dtime) then + call markreset_accum_field(field_name, p) + else if (in_accumulation_season) then + rbufslp(p) = max(0._r8, min(max_accum, & + this%t_ref2m_patch(p)-(SHR_CONST_TKFRZ + basetemp_r8))) * dtime/SHR_CONST_CDAY + else + rbufslp(p) = 0._r8 ! keeps gdd unchanged outside accumulation season + end if + end do + + ! Save + call update_accum_field (trim(field_name), rbufslp, nstep) + call extract_accum_field (trim(field_name), gddx_patch, nstep) + + end associate + end subroutine UpdateAccVars_CropGDDs + + !----------------------------------------------------------------------- + subroutine UpdateAccVars (this, bounds, crop_inst) + ! + ! USES + use shr_const_mod , only : SHR_CONST_TKFRZ use clm_time_manager , only : get_step_size, get_nstep, is_end_curr_day, get_curr_date, is_end_curr_year - use accumulMod , only : update_accum_field, extract_accum_field, accumResetVal + use accumulMod , only : update_accum_field, extract_accum_field, markreset_accum_field use CNSharedParamsMod, only : upper_soil_layer + use CropType , only : crop_type ! ! !ARGUMENTS: class(temperature_type) :: this type(bounds_type) , intent(in) :: bounds + type(crop_type), intent(inout) :: crop_inst ! ! !LOCAL VARIABLES: - integer :: m,g,l,c,p ! indices + integer :: m,l,c,p ! indices integer :: ier ! error status integer :: dtime ! timestep size [seconds] integer :: nstep ! timestep number @@ -1392,6 +1519,7 @@ subroutine UpdateAccVars (this, bounds) dtime = get_step_size() nstep = get_nstep() + ! SSR 2024-06-13: This should probably be changed to _prev_ call get_curr_date (year, month, day, secs) ! Allocate needed dynamic memory for single level pft field @@ -1535,69 +1663,25 @@ subroutine UpdateAccVars (this, bounds) call update_accum_field ('TDM10', rbufslp, nstep) call extract_accum_field ('TDM10', this%t_a10min_patch, nstep) - - ! Accumulate and extract GDD0 - - do p = begp,endp - ! Avoid unnecessary calculations over inactive points - if (patch%active(p)) then - g = patch%gridcell(p) - if (month==1 .and. day==1 .and. secs==dtime) then - rbufslp(p) = accumResetVal ! reset gdd - else if (( month > 3 .and. month < 10 .and. grc%latdeg(g) >= 0._r8) .or. & - ((month > 9 .or. month < 4) .and. grc%latdeg(g) < 0._r8) ) then - rbufslp(p) = max(0._r8, min(26._r8, this%t_ref2m_patch(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) - call extract_accum_field ('GDD0', this%gdd0_patch, nstep) + call this%UpdateAccVars_CropGDDs(rbufslp, begp, endp, month, day, secs, dtime, nstep, 0, this%gdd0_patch, crop_inst) ! Accumulate and extract GDD8 - - do p = begp,endp - ! Avoid unnecessary calculations over inactive points - if (patch%active(p)) then - g = patch%gridcell(p) - if (month==1 .and. day==1 .and. secs==dtime) then - rbufslp(p) = accumResetVal ! reset gdd - else if (( month > 3 .and. month < 10 .and. grc%latdeg(g) >= 0._r8) .or. & - ((month > 9 .or. month < 4) .and. grc%latdeg(g) < 0._r8) ) then - rbufslp(p) = max(0._r8, min(30._r8, & - this%t_ref2m_patch(p)-(SHR_CONST_TKFRZ + 8._r8))) * 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 ('GDD8', rbufslp, nstep) - call extract_accum_field ('GDD8', this%gdd8_patch, nstep) + call this%UpdateAccVars_CropGDDs(rbufslp, begp, endp, month, day, secs, dtime, nstep, 8, this%gdd8_patch, crop_inst) ! Accumulate and extract GDD10 - - do p = begp,endp - ! Avoid unnecessary calculations over inactive points - if (patch%active(p)) then - g = patch%gridcell(p) - if (month==1 .and. day==1 .and. secs==dtime) then - rbufslp(p) = accumResetVal ! reset gdd - else if (( month > 3 .and. month < 10 .and. grc%latdeg(g) >= 0._r8) .or. & - ((month > 9 .or. month < 4) .and. grc%latdeg(g) < 0._r8) ) then - rbufslp(p) = max(0._r8, min(30._r8, & - this%t_ref2m_patch(p)-(SHR_CONST_TKFRZ + 10._r8))) * 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 ('GDD10', rbufslp, nstep) - call extract_accum_field ('GDD10', this%gdd10_patch, nstep) + call this%UpdateAccVars_CropGDDs(rbufslp, begp, endp, month, day, secs, dtime, nstep, 10, this%gdd10_patch, crop_inst) ! Accumulate and extract running 20-year means if (is_end_curr_year()) then + ! Flush, if needed + if (flush_gdd20) then + write(iulog, *) 'Flushing GDD20 variables' + call markreset_accum_field('GDD020') + call markreset_accum_field('GDD820') + call markreset_accum_field('GDD1020') + flush_gdd20 = .false. + end if call update_accum_field ('GDD020', this%gdd0_patch, nstep) call extract_accum_field ('GDD020', this%gdd020_patch, nstep) call update_accum_field ('GDD820', this%gdd8_patch, nstep) diff --git a/src/cpl/share_esmf/cropcalStreamMod.F90 b/src/cpl/share_esmf/cropcalStreamMod.F90 index 0ea63f2c6d..5587641c21 100644 --- a/src/cpl/share_esmf/cropcalStreamMod.F90 +++ b/src/cpl/share_esmf/cropcalStreamMod.F90 @@ -14,7 +14,9 @@ module cropcalStreamMod use decompMod , only : bounds_type use abortutils , only : endrun use clm_varctl , only : iulog + use clm_varctl , only : use_crop use clm_varctl , only : use_cropcal_rx_swindows, use_cropcal_rx_cultivar_gdds, use_cropcal_streams + use clm_varctl , only : adapt_cropcal_rx_cultivar_gdds use clm_varpar , only : mxpft use clm_varpar , only : mxsowings use perf_mod , only : t_startf, t_stopf @@ -35,14 +37,26 @@ module cropcalStreamMod integer, allocatable :: g_to_ig(:) ! Array matching gridcell index to data index 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 + type(shr_strdata_type) :: sdat_cropcal_cultivar_gdds ! maturity requirement input data stream + type(shr_strdata_type) :: sdat_cropcal_gdd20_baseline ! GDD20 baseline input data stream + type(shr_strdata_type) :: sdat_cropcal_gdd20_season_start ! gdd20 season start input data stream + type(shr_strdata_type) :: sdat_cropcal_gdd20_season_end ! gdd20 season end input data stream character(len=CS), allocatable :: stream_varnames_sdate(:) ! used for both start and end dates character(len=CS), allocatable :: stream_varnames_cultivar_gdds(:) + character(len=CS), allocatable :: stream_varnames_gdd20_baseline(:) + character(len=CS), allocatable :: stream_varnames_gdd20_season_enddate(:) ! start uses stream_varnames_sdate integer :: ncft ! Number of crop functional types (excl. generic crops) 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=CL) :: stream_fldFileName_gdd20_baseline ! GDD20 baseline stream filename to read + logical :: cropcals_rx ! Used only for setting input files in namelist; does nothing in code, but needs to be here so namelist read doesn't crash + logical :: cropcals_rx_adapt ! Used only for setting input files in namelist; does nothing in code, but needs to be here so namelist read doesn't crash + logical :: stream_gdd20_seasons ! Read start and end dates for gdd20 seasons from streams instead of using hemisphere-specific values + logical :: allow_invalid_gdd20_season_inputs ! Fall back on hemisphere "warm periods" in cases of invalid values in stream_fldFileName_gdd20_season_start and _end? + character(len=CL) :: stream_fldFileName_gdd20_season_start ! Stream filename to read for start of gdd20 season + character(len=CL) :: stream_fldFileName_gdd20_season_end ! Stream filename to read for end of gdd20 season character(len=*), parameter :: sourcefile = & __FILE__ @@ -67,9 +81,12 @@ subroutine cropcal_init(bounds) ! ! !LOCAL VARIABLES: integer :: i,n,ivt ! index - integer :: stream_year_first_cropcal ! first year in crop calendar streams to use - integer :: stream_year_last_cropcal ! last year in crop calendar streams to use - integer :: model_year_align_cropcal ! align stream_year_first_cropcal with + integer :: stream_year_first_cropcal_swindows ! first year in sowing window streams to use + integer :: stream_year_last_cropcal_swindows ! last year in sowing window streams to use + integer :: model_year_align_cropcal_swindows ! alignment year for sowing window streams + integer :: stream_year_first_cropcal_cultivar_gdds ! first year in cultivar gdd stream to use + integer :: stream_year_last_cropcal_cultivar_gdds ! last year in cultivar gdd stream to use + integer :: model_year_align_cropcal_cultivar_gdds ! alignment year for cultivar gdd stream integer :: nu_nml ! unit for namelist file integer :: nml_error ! namelist i/o error flag character(len=CL) :: stream_meshfile_cropcal ! crop calendar stream meshfile @@ -83,32 +100,54 @@ subroutine cropcal_init(bounds) ! deal with namelist variables here in init ! namelist /cropcal_streams/ & - stream_year_first_cropcal, & - stream_year_last_cropcal, & - model_year_align_cropcal, & + stream_year_first_cropcal_swindows, & + stream_year_last_cropcal_swindows, & + model_year_align_cropcal_swindows, & + stream_year_first_cropcal_cultivar_gdds, & + stream_year_last_cropcal_cultivar_gdds, & + model_year_align_cropcal_cultivar_gdds, & allow_invalid_swindow_inputs, & stream_fldFileName_swindow_start, & stream_fldFileName_swindow_end, & stream_fldFileName_cultivar_gdds, & - stream_meshfile_cropcal + stream_fldFileName_gdd20_baseline, & + stream_meshfile_cropcal, & + cropcals_rx, & + cropcals_rx_adapt, & + stream_gdd20_seasons, & + allow_invalid_gdd20_season_inputs, & + stream_fldFileName_gdd20_season_start, & + stream_fldFileName_gdd20_season_end ! Default values for namelist - 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 + stream_year_first_cropcal_swindows = 1 ! first year in sowing window streams to use + stream_year_last_cropcal_swindows = 1 ! last year in sowing window streams to use + model_year_align_cropcal_swindows = 1 ! alignment year for sowing window streams + stream_year_first_cropcal_cultivar_gdds = 1 ! first year in cultivar gdd stream to use + stream_year_last_cropcal_cultivar_gdds = 1 ! last year in cultivar gdd stream to use + model_year_align_cropcal_cultivar_gdds = 1 ! alignment year for cultivar gdd stream allow_invalid_swindow_inputs = .false. stream_meshfile_cropcal = '' stream_fldFileName_swindow_start = '' stream_fldFileName_swindow_end = '' stream_fldFileName_cultivar_gdds = '' + stream_fldFileName_gdd20_baseline = '' + stream_gdd20_seasons = .false. + allow_invalid_gdd20_season_inputs = .false. + stream_fldFileName_gdd20_season_start = '' + stream_fldFileName_gdd20_season_end = '' ! Will need modification to work with mxsowings > 1 ncft = mxpft - npcropmin + 1 ! Ignores generic crops allocate(stream_varnames_sdate(ncft)) allocate(stream_varnames_cultivar_gdds(ncft)) + allocate(stream_varnames_gdd20_baseline(ncft)) + allocate(stream_varnames_gdd20_season_enddate(ncft)) do n = 1,ncft ivt = npcropmin + n - 1 write(stream_varnames_sdate(n),'(a,i0)') "sdate1_",ivt write(stream_varnames_cultivar_gdds(n),'(a,i0)') "gdd1_",ivt + write(stream_varnames_gdd20_baseline(n),'(a,i0)') "gdd20bl_",ivt + write(stream_varnames_gdd20_season_enddate(n),'(a,i0)') "hdate1_",ivt end do ! Read cropcal_streams namelist @@ -125,29 +164,47 @@ subroutine cropcal_init(bounds) end if close(nu_nml) endif - 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_year_first_cropcal_swindows , mpicom) + call shr_mpi_bcast(stream_year_last_cropcal_swindows , mpicom) + call shr_mpi_bcast(model_year_align_cropcal_swindows , mpicom) + call shr_mpi_bcast(stream_year_first_cropcal_cultivar_gdds, mpicom) + call shr_mpi_bcast(stream_year_last_cropcal_cultivar_gdds , mpicom) + call shr_mpi_bcast(model_year_align_cropcal_cultivar_gdds , 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_fldFileName_gdd20_baseline, mpicom) call shr_mpi_bcast(stream_meshfile_cropcal , mpicom) + call shr_mpi_bcast(stream_gdd20_seasons, mpicom) + call shr_mpi_bcast(allow_invalid_gdd20_season_inputs, mpicom) + call shr_mpi_bcast(stream_fldFileName_gdd20_season_start, mpicom) + call shr_mpi_bcast(stream_fldFileName_gdd20_season_end, mpicom) if (masterproc) then write(iulog,*) write(iulog,*) 'cropcal_stream settings:' - 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,i8)') ' stream_year_first_cropcal_swindows = ',stream_year_first_cropcal_swindows + write(iulog,'(a,i8)') ' stream_year_last_cropcal_swindows = ',stream_year_last_cropcal_swindows + write(iulog,'(a,i8)') ' model_year_align_cropcal_swindows = ',model_year_align_cropcal_swindows + write(iulog,'(a,i8)') ' stream_year_first_cropcal_cultivar_gdds = ',stream_year_first_cropcal_cultivar_gdds + write(iulog,'(a,i8)') ' stream_year_last_cropcal_cultivar_gdds = ',stream_year_last_cropcal_cultivar_gdds + write(iulog,'(a,i8)') ' model_year_align_cropcal_cultivar_gdds = ',model_year_align_cropcal_cultivar_gdds 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_fldFileName_gdd20_baseline = ',trim(stream_fldFileName_gdd20_baseline) write(iulog,'(a,a)' ) ' stream_meshfile_cropcal = ',trim(stream_meshfile_cropcal) + write(iulog,'(a,l1)') ' stream_gdd20_seasons = ',stream_gdd20_seasons + write(iulog,'(a,l1)') ' allow_invalid_gdd20_season_inputs = ',allow_invalid_gdd20_season_inputs + write(iulog,'(a,a)' ) ' stream_fldFileName_gdd20_season_start = ',stream_fldFileName_gdd20_season_start + write(iulog,'(a,a)' ) ' stream_fldFileName_gdd20_season_end = ',stream_fldFileName_gdd20_season_end do n = 1,ncft write(iulog,'(a,a)' ) ' stream_varnames_sdate = ',trim(stream_varnames_sdate(n)) write(iulog,'(a,a)' ) ' stream_varnames_cultivar_gdds = ',trim(stream_varnames_cultivar_gdds(n)) + write(iulog,'(a,a)' ) ' stream_varnames_gdd20_season_enddate = ',trim(stream_varnames_gdd20_season_enddate(n)) + write(iulog,'(a,a)' ) ' stream_varnames_gdd20_baseline = ',trim(stream_varnames_gdd20_baseline(n)) end do write(iulog,*) endif @@ -155,9 +212,12 @@ subroutine cropcal_init(bounds) ! 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_swindows .or. use_cropcal_rx_cultivar_gdds + adapt_cropcal_rx_cultivar_gdds = stream_fldFileName_gdd20_baseline /= '' + use_cropcal_streams = .false. ! Will be set to true if any file is read if (use_cropcal_rx_swindows) then + use_cropcal_streams = .true. + ! 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, & @@ -172,9 +232,9 @@ subroutine cropcal_init(bounds) stream_filenames = (/trim(stream_fldFileName_swindow_start)/), & 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_yearFirst = stream_year_first_cropcal_swindows, & + stream_yearLast = stream_year_last_cropcal_swindows, & + stream_yearAlign = model_year_align_cropcal_swindows, & stream_offset = cropcal_offset, & stream_taxmode = 'extend', & stream_dtlimit = 1.0e30_r8, & @@ -199,9 +259,9 @@ subroutine cropcal_init(bounds) 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_yearFirst = stream_year_first_cropcal_swindows, & + stream_yearLast = stream_year_last_cropcal_swindows, & + stream_yearAlign = model_year_align_cropcal_swindows, & stream_offset = cropcal_offset, & stream_taxmode = 'extend', & stream_dtlimit = 1.0e30_r8, & @@ -216,6 +276,7 @@ subroutine cropcal_init(bounds) ! Initialize the cdeps data type sdat_cropcal_cultivar_gdds ! NOTE: stream_dtlimit 1.5 didn't work for some reason if (use_cropcal_rx_cultivar_gdds) then + use_cropcal_streams = .true. call shr_strdata_init_from_inline(sdat_cropcal_cultivar_gdds, & my_task = iam, & logunit = iulog, & @@ -228,9 +289,9 @@ subroutine cropcal_init(bounds) stream_filenames = (/trim(stream_fldFileName_cultivar_gdds)/), & stream_fldlistFile = stream_varnames_cultivar_gdds, & stream_fldListModel = stream_varnames_cultivar_gdds, & - stream_yearFirst = stream_year_first_cropcal, & - stream_yearLast = stream_year_last_cropcal, & - stream_yearAlign = model_year_align_cropcal, & + stream_yearFirst = stream_year_first_cropcal_cultivar_gdds,& + stream_yearLast = stream_year_last_cropcal_cultivar_gdds, & + stream_yearAlign = model_year_align_cropcal_cultivar_gdds, & stream_offset = cropcal_offset, & stream_taxmode = 'extend', & stream_dtlimit = 1.0e30_r8, & @@ -242,6 +303,109 @@ subroutine cropcal_init(bounds) end if end if + ! Initialize the cdeps data type sdat_cropcal_gdd20_baseline + ! NOTE: Hard-coded to one particular year because it should NOT vary over time. Note that the + ! particular year chosen doesn't matter. Users can base their file on whatever baseline they + ! want; they just need to put 2000 on the time axis. + if (adapt_cropcal_rx_cultivar_gdds) then + use_cropcal_streams = .true. + call shr_strdata_init_from_inline(sdat_cropcal_gdd20_baseline, & + 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 = 'nn', & + stream_filenames = (/trim(stream_fldFileName_gdd20_baseline)/), & + stream_fldlistFile = stream_varnames_gdd20_baseline, & + stream_fldListModel = stream_varnames_gdd20_baseline, & + stream_yearFirst = 2000, & + stream_yearLast = 2000, & + stream_yearAlign = 2000, & + stream_offset = cropcal_offset, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = cropcal_tintalgo, & + stream_name = 'GDD20 baseline 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 + end if + + if (stream_gdd20_seasons) then + use_cropcal_streams = .true. + + ! Initialize the cdeps data type sdat_cropcal_gdd20_season_start + ! NOTE: Hard-coded to one particular year because it should NOT vary over time. Note that the + ! particular year chosen doesn't matter. + call shr_strdata_init_from_inline(sdat_cropcal_gdd20_season_start, & + 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_gdd20_season_start)/), & + stream_fldlistFile = stream_varnames_sdate, & + stream_fldListModel = stream_varnames_sdate, & + stream_yearFirst = 2000, & + stream_yearLast = 2000, & + stream_yearAlign = 2000, & + stream_offset = cropcal_offset, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = cropcal_tintalgo, & + stream_name = 'gdd20 season 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_gdd20_season_end + ! NOTE: Hard-coded to one particular year because it should NOT vary over time. Note that the + ! particular year chosen doesn't matter. + call shr_strdata_init_from_inline(sdat_cropcal_gdd20_season_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_gdd20_season_end)/), & + stream_fldlistFile = stream_varnames_gdd20_season_enddate, & + stream_fldListModel = stream_varnames_gdd20_season_enddate, & + stream_yearFirst = 2000, & + stream_yearLast = 2000, & + stream_yearAlign = 2000, & + stream_offset = cropcal_offset, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = cropcal_tintalgo, & + stream_name = 'gdd20 season 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 + + end if + + if (masterproc) then + write(iulog,*) + write(iulog,*) 'cropcal_stream DERIVED settings:' + write(iulog,'(a,l1)') ' use_cropcal_rx_swindows = ',use_cropcal_rx_swindows + write(iulog,'(a,l1)') ' use_cropcal_rx_cultivar_gdds = ',use_cropcal_rx_cultivar_gdds + write(iulog,'(a,l1)') ' adapt_cropcal_rx_cultivar_gdds = ',adapt_cropcal_rx_cultivar_gdds + write(iulog,'(a,l1)') ' use_cropcal_streams = ',use_cropcal_streams + write(iulog,*) + endif + end subroutine cropcal_init !================================================================ @@ -258,6 +422,7 @@ subroutine cropcal_advance( bounds ) ! ! !LOCAL VARIABLES: integer :: g, ig ! Indices + integer :: begg, endg ! gridcell bounds integer :: year ! year (0, ...) for nstep+1 integer :: mon ! month (1, ..., 12) for nstep+1 integer :: day ! day of month (1, ..., 31) for nstep+1 @@ -266,6 +431,9 @@ subroutine cropcal_advance( bounds ) integer :: rc !----------------------------------------------------------------------- + begg = bounds%begg + endg = bounds%endg + call get_curr_date(year, mon, day, sec) mcdate = year*10000 + mon*100 + day if (use_cropcal_rx_swindows) then @@ -285,10 +453,31 @@ subroutine cropcal_advance( bounds ) end if end if + ! The following should not have an associated time axis, but still need to be here + ! - GDD20 baseline values + ! - GDD20 season start dates + ! - GDD20 season end dates + if (adapt_cropcal_rx_cultivar_gdds) then + call shr_strdata_advance(sdat_cropcal_gdd20_baseline, 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 + end if + if (stream_gdd20_seasons) then + call shr_strdata_advance(sdat_cropcal_gdd20_season_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_gdd20_season_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 + end if + if ( .not. allocated(g_to_ig) )then - allocate (g_to_ig(bounds%begg:bounds%endg) ) + allocate (g_to_ig(begg:endg) ) ig = 0 - do g = bounds%begg,bounds%endg + do g = begg,endg ig = ig+1 g_to_ig(g) = ig end do @@ -298,7 +487,7 @@ end subroutine cropcal_advance !================================================================ - subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) + subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, crop_inst) ! ! Interpolate data stream information for crop calendars. ! @@ -314,6 +503,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_pcropp ! number of prog. crop patches in filter integer , intent(in) :: filter_pcropp(:) ! filter for prognostic crop patches + logical , intent(in) :: init ! is this being called as initialization? type(crop_type) , intent(inout) :: crop_inst ! ! !LOCAL VARIABLES: @@ -321,39 +511,48 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) integer :: nc, fp integer :: dayspyr integer :: n, g - integer :: lsize integer :: rc + integer :: begg, endg integer :: begp, endp real(r8), pointer :: dataptr1d_swindow_start(:) real(r8), pointer :: dataptr1d_swindow_end (:) real(r8), pointer :: dataptr1d_cultivar_gdds(:) + real(r8), pointer :: dataptr1d_gdd20_baseline(:) + real(r8), pointer :: dataptr1d_gdd20_season_start(:) + real(r8), pointer :: dataptr1d_gdd20_season_end (:) real(r8), pointer :: dataptr2d_swindow_start(:,:) real(r8), pointer :: dataptr2d_swindow_end (:,:) real(r8), pointer :: dataptr2d_cultivar_gdds(:,:) + real(r8), pointer :: dataptr2d_gdd20_baseline(:,:) + real(r8), pointer :: dataptr2d_gdd20_season_start(:,:) + real(r8), pointer :: dataptr2d_gdd20_season_end (:,:) !----------------------------------------------------------------------- associate( & - starts => crop_inst%rx_swindow_starts_thisyr_patch, & - ends => crop_inst%rx_swindow_ends_thisyr_patch & + swindow_starts => crop_inst%rx_swindow_starts_thisyr_patch, & + swindow_ends => crop_inst%rx_swindow_ends_thisyr_patch, & + gdd20_season_starts => crop_inst%gdd20_season_start_patch, & + gdd20_season_ends => crop_inst%gdd20_season_end_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__) + begg = bounds%begg + endg = bounds%endg + SHR_ASSERT_FL( (lbound(g_to_ig,1) <= begg ), sourcefile, __LINE__) + SHR_ASSERT_FL( (ubound(g_to_ig,1) >= endg ), sourcefile, __LINE__) ! Get pointer for stream data that is time and spatially interpolate to model time and grid ! Place all data from each type into a temporary 2d array - lsize = bounds%endg - bounds%begg + 1 begp = bounds%begp - endp= bounds%endp + endp = bounds%endp 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 + allocate(dataptr2d_swindow_start(begg:endg, ncft)) + dataptr2d_swindow_start(begg:endg,:) = -1._r8 + allocate(dataptr2d_swindow_end (begg:endg, ncft)) + dataptr2d_swindow_end(begg:endg,:) = -1._r8 if (use_cropcal_rx_swindows) then ! Starting with npcropmin will skip generic crops do n = 1, ncft @@ -369,7 +568,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) 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 + do g = begg, endg ! 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 & @@ -392,8 +591,8 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) n = ivt - npcropmin + 1 ! vegetated pft ig = g_to_ig(patch%gridcell(p)) - starts(p,1) = dataptr2d_swindow_start(ig,n) - ends(p,1) = dataptr2d_swindow_end (ig,n) + swindow_starts(p,1) = dataptr2d_swindow_start(ig,n) + swindow_ends(p,1) = dataptr2d_swindow_end (ig,n) else write(iulog,'(a,i0)') 'cropcal_interp(), prescribed sowing windows: Crop patch has ivt ',ivt call ESMF_Finalize(endflag=ESMF_END_ABORT) @@ -402,23 +601,23 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! 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(begp:endp,2:mxsowings) <= ends(begp:endp,1:mxsowings-1) .and. & - ends(begp:endp,2:mxsowings) >= 1)) then + if (any(swindow_ends(begp:endp,2:mxsowings) <= swindow_ends(begp:endp,1:mxsowings-1) .and. & + swindow_ends(begp:endp,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(begp:endp,:) < 1 .or. ends(begp:endp,:) < 1)) then + if (any(swindow_starts(begp:endp,:) < 1 .or. swindow_ends(begp:endp,:) < 1)) then ! Fail if not allowing fallback to paramfile sowing windows - if ((.not. allow_invalid_swindow_inputs) .and. any(all(starts(begp:endp,:) < 1, dim=2) .and. patch%wtgcell > 0._r8 .and. patch%itype >= npcropmin)) then + if ((.not. allow_invalid_swindow_inputs) .and. any(all(swindow_starts(begp:endp,:) < 1, dim=2) .and. patch%wtgcell(begp:endp) > 0._r8 .and. patch%itype(begp:endp) >= 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 + if (ivt == patch%itype(p) .and. patch%wtgcell(p) > 0._r8 .and. all(swindow_starts(p,:) < 1)) then write(iulog, *) ' ',pftname(ivt),' (',ivt,')' exit ! Stop looking for patches of this type end if @@ -427,7 +626,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) 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(begp:endp,:) >= 1 .and. ends(begp:endp,:) < 1) .or. (starts(begp:endp,:) < 1 .and. ends(begp:endp,:) >= 1))) then + else if (any((swindow_starts(begp:endp,:) >= 1 .and. swindow_ends(begp:endp,:) < 1) .or. (swindow_starts(begp:endp,:) < 1 .and. swindow_ends(begp:endp,:) >= 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 @@ -437,7 +636,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) deallocate(dataptr2d_swindow_start) deallocate(dataptr2d_swindow_end) - allocate(dataptr2d_cultivar_gdds(lsize, ncft)) + allocate(dataptr2d_cultivar_gdds(begg:endg, ncft)) if (use_cropcal_rx_cultivar_gdds) then ! Read prescribed cultivar GDDs from input files ! Starting with npcropmin will skip generic crops @@ -450,7 +649,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! 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 + do g = begg, endg ! If read-in value is invalid, have PlantCrop() set gddmaturity to PFT-default value. if (dataptr1d_cultivar_gdds(g) < 0 .or. dataptr1d_cultivar_gdds(g) > 1000000._r8) then @@ -478,8 +677,8 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! vegetated pft ig = g_to_ig(patch%gridcell(p)) - if (ig > lsize) then - write(iulog,'(a,i0,a,i0,a)') 'ig (',ig,') > lsize (',lsize,')' + if (ig < begg .or. ig > endg) then + write(iulog,'(a,i0,a,i0,a)') 'ig (',ig,') < begg (',begg,') or > endg (',endg,')' call ESMF_Finalize(endflag=ESMF_END_ABORT) end if @@ -490,11 +689,147 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif end do - write(iulog,*) 'cropcal_interp(): Reading cultivar_gdds file DONE' end if ! use_cropcal_rx_cultivar_gdds deallocate(dataptr2d_cultivar_gdds) + allocate(dataptr2d_gdd20_baseline(begg:endg, ncft)) + if (adapt_cropcal_rx_cultivar_gdds) then + ! Read GDD20 baselines from input files + ! Starting with npcropmin will skip generic crops + do n = 1, ncft + call dshr_fldbun_getFldPtr(sdat_cropcal_gdd20_baseline%pstrm(1)%fldbun_model, trim(stream_varnames_gdd20_baseline(n)), & + fldptr1=dataptr1d_gdd20_baseline, 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 = begg, endg + dataptr2d_gdd20_baseline(g,n) = dataptr1d_gdd20_baseline(g) + end do + end do + + ! Set gdd20_baseline_patch for each gridcell/patch combination + do fp = 1, num_pcropp + p = filter_pcropp(fp) + + ivt = patch%itype(p) + ! Will skip generic crops + if (ivt >= npcropmin) then + n = ivt - npcropmin + 1 + + if (n > ncft) then + write(iulog,'(a,i0,a,i0,a)') 'n (',n,') > ncft (',ncft,')' + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! vegetated pft + ig = g_to_ig(patch%gridcell(p)) + + if (ig < begg .or. ig > endg) then + write(iulog,'(a,i0,a,i0,a)') 'ig (',ig,') < begg (',begg,') or > endg (',endg,')' + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + crop_inst%gdd20_baseline_patch(p) = dataptr2d_gdd20_baseline(ig,n) + + else + write(iulog,'(a,i0)') 'cropcal_interp(), rx_gdd20_baseline: Crop patch has ivt ',ivt + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + end do + end if ! adapt_cropcal_rx_cultivar_gdds + + deallocate(dataptr2d_gdd20_baseline) + + + ! Read prescribed gdd20 season start dates from input files + allocate(dataptr2d_gdd20_season_start(begg:endg, ncft)) + dataptr2d_gdd20_season_start(begg:endg,:) = -1._r8 + allocate(dataptr2d_gdd20_season_end (begg:endg, ncft)) + dataptr2d_gdd20_season_end(begg:endg,:) = -1._r8 + if (stream_gdd20_seasons) then + ! Starting with npcropmin will skip generic crops + do n = 1, ncft + call dshr_fldbun_getFldPtr(sdat_cropcal_gdd20_season_start%pstrm(1)%fldbun_model, trim(stream_varnames_sdate(n)), & + fldptr1=dataptr1d_gdd20_season_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_gdd20_season_end%pstrm(1)%fldbun_model, trim(stream_varnames_gdd20_season_enddate(n)), & + fldptr1=dataptr1d_gdd20_season_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 = begg, endg + + ! If read-in value is invalid, set to -1. Will be handled later in this subroutine. + if (dataptr1d_gdd20_season_start(g) <= 0 .or. dataptr1d_gdd20_season_start(g) > 366 & + .or. dataptr1d_gdd20_season_start(g) /= dataptr1d_gdd20_season_start(g)) then + dataptr1d_gdd20_season_start(g) = -1 + end if + if (dataptr1d_gdd20_season_end(g) <= 0 .or. dataptr1d_gdd20_season_end(g) > 366 & + .or. dataptr1d_gdd20_season_end(g) /= dataptr1d_gdd20_season_end(g)) then + dataptr1d_gdd20_season_end (g) = -1 + end if + + dataptr2d_gdd20_season_start(g,n) = dataptr1d_gdd20_season_start(g) + dataptr2d_gdd20_season_end (g,n) = dataptr1d_gdd20_season_end (g) + end do + end do + + ! Set gdd20 season for each gridcell/patch combination + do fp = 1, num_pcropp + p = filter_pcropp(fp) + ivt = patch%itype(p) + ! Will skip generic crops + if (ivt >= npcropmin) then + n = ivt - npcropmin + 1 + ! vegetated pft + ig = g_to_ig(patch%gridcell(p)) + + gdd20_season_starts(p) = real(dataptr2d_gdd20_season_start(ig,n), r8) + gdd20_season_ends(p) = real(dataptr2d_gdd20_season_end (ig,n), r8) + else + write(iulog,'(a,i0)') 'cropcal_interp(), gdd20 seasons: Crop patch has ivt ',ivt + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + end do + + ! Handle invalid gdd20 season values + if (any(gdd20_season_starts(begp:endp) < 1._r8 .or. gdd20_season_ends(begp:endp) < 1._r8)) then + ! Fail if not allowing fallback to paramfile sowing windows. Only need to check for + ! values < 1 because values outside [1, 366] are set to -1 above. + if ((.not. allow_invalid_gdd20_season_inputs) .and. any(gdd20_season_starts(begp:endp) < 1._r8 .and. patch%wtgcell(begp:endp) > 0._r8 .and. patch%itype(begp:endp) >= npcropmin)) then + write(iulog, *) 'At least one crop in one gridcell has invalid gdd20 season start and/or end date(s). To ignore and fall back to paramfile sowing windows for such crop-gridcells, set allow_invalid_gdd20_season_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. gdd20_season_starts(p) < 1._r8) 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 gdd20 season start date is given without an end date (or vice versa) + else if (any((gdd20_season_starts(begp:endp) >= 1._r8 .and. gdd20_season_ends(begp:endp) < 1._r8) .or. (gdd20_season_starts(begp:endp) < 1._r8 .and. gdd20_season_ends(begp:endp) >= 1._r8))) then + write(iulog, *) 'Every gdd20 season start date must have a corresponding end date.' + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + end if + + end if ! stream_gdd20_seasons + deallocate(dataptr2d_gdd20_season_start) + deallocate(dataptr2d_gdd20_season_end) + + end associate end subroutine cropcal_interp diff --git a/src/main/accumulMod.F90 b/src/main/accumulMod.F90 index e328632501..0c2462c1f8 100644 --- a/src/main/accumulMod.F90 +++ b/src/main/accumulMod.F90 @@ -40,6 +40,8 @@ module accumulMod public :: print_accum_fields ! Print info about accumulator fields public :: extract_accum_field ! Extracts the current value of an accumulator field public :: update_accum_field ! Update the current value of an accumulator field + public :: markreset_accum_field ! Mark (value(s) in) an accumulator field as needing to be reset + public :: get_accum_reset ! Get reset array of accumulator public :: clean_accum_fields ! Deallocate space and reset accum fields list interface extract_accum_field @@ -76,6 +78,7 @@ module accumulMod logical, pointer :: active(:)!whether each point (patch, column, etc.) is active real(r8) :: initval !initial value of accumulated field real(r8), pointer :: val(:,:) !accumulated field + logical , pointer :: reset(:,:) !whether accumulated field needs to be reset integer :: period !field accumulation period (in model time steps) logical :: scale_by_thickness ! true/false flag to scale vertically interpolated variable by soil thickness or not character(len=128) :: old_name !previous name of variable (may be present in restart files) @@ -85,6 +88,8 @@ module accumulMod ! different reset points for different levels. integer, pointer :: nsteps(:,:)!number of steps each point has accumulated, since last reset time + integer, pointer :: ndays_reset_shifted(:,:)!accumulated number of days that resetting timeavg field has moved it out of sync with model timestep + ! NOTE(wjs, 2017-12-03) We should convert this to fully object-oriented (with ! inheritance / polymorphism). For now, in the interest of time, I'm going with a ! semi-object-oriented solution of using procedure pointers. @@ -112,8 +117,6 @@ subroutine update_accum_field_interface(this, level, nstep, field) end subroutine update_accum_field_interface end interface - real(r8), parameter, public :: accumResetVal = -99999._r8 ! used to do an annual reset ( put in for bug 1858) - integer, parameter :: ACCTYPE_TIMEAVG = 1 integer, parameter :: ACCTYPE_RUNMEAN = 2 integer, parameter :: ACCTYPE_RUNACCUM = 3 @@ -295,9 +298,15 @@ subroutine init_accum_field (name, units, desc, & end if accum(nf)%val(beg1d:end1d,1:numlev) = init_value + allocate(accum(nf)%reset(beg1d:end1d,numlev)) + accum(nf)%reset(beg1d:end1d,1:numlev) = .false. + allocate(accum(nf)%nsteps(beg1d:end1d,numlev)) accum(nf)%nsteps(beg1d:end1d,1:numlev) = 0 + allocate(accum(nf)%ndays_reset_shifted(beg1d:end1d,numlev)) + accum(nf)%ndays_reset_shifted(beg1d:end1d,1:numlev) = 0 + end subroutine init_accum_field !------------------------------------------------------------------------ @@ -464,6 +473,7 @@ subroutine extract_accum_field_timeavg(this, level, nstep, field) ! !LOCAL VARIABLES: integer :: begi,endi !subgrid beginning,ending indices integer :: k, kf + integer :: effective_nstep ! Timestep accounting for resets character(len=*), parameter :: subname = 'extract_accum_field_basic' !----------------------------------------------------------------------- @@ -472,17 +482,15 @@ subroutine extract_accum_field_timeavg(this, level, nstep, field) endi = this%end1d SHR_ASSERT_FL((size(field) == endi-begi+1), sourcefile, __LINE__) - if (mod(nstep,this%period) == 0) then - do k = begi, endi - kf = k - begi + 1 + do k = begi, endi + kf = k - begi + 1 + effective_nstep = nstep - this%ndays_reset_shifted(k,level) + if (mod(effective_nstep,this%period) == 0) then field(kf) = this%val(k,level) - end do - else - do k = begi, endi - kf = k - begi + 1 + else field(kf) = spval - end do - end if + end if + end do end subroutine extract_accum_field_timeavg @@ -572,6 +580,8 @@ subroutine update_accum_field_timeavg(this, level, nstep, field) ! !LOCAL VARIABLES: integer :: begi,endi !subgrid beginning,ending indices integer :: k, kf + logical :: time_to_reset + integer :: effective_nstep ! Timestep accounting for resets character(len=*), parameter :: subname = 'update_accum_field_timeavg' !----------------------------------------------------------------------- @@ -583,14 +593,21 @@ subroutine update_accum_field_timeavg(this, level, nstep, field) ! time average field: reset every accumulation period; normalize at end of ! accumulation period - if ((mod(nstep,this%period) == 1 .or. this%period == 1) .and. (nstep /= 0))then - do k = begi,endi - if (this%active(k)) then - this%val(k,level) = 0._r8 - this%nsteps(k,level) = 0 + do k = begi,endi + effective_nstep = nstep - this%ndays_reset_shifted(k,level) + time_to_reset = (mod(effective_nstep,this%period) == 1 .or. this%period == 1) .and. effective_nstep /= 0 + if (this%active(k) .and. (time_to_reset .or. this%reset(k,level))) then + if (this%reset(k,level) .and. .not. time_to_reset) then + this%ndays_reset_shifted(k,level) = this%ndays_reset_shifted(k,level) + this%nsteps(k,level) end if - end do - end if + this%val(k,level) = this%initval + this%nsteps(k,level) = 0 + this%reset(k,level) = .false. + end if + end do + + ! Ignore reset requests that occurred when patch was inactive + this%reset(begi:endi,level) = .false. do k = begi,endi if (this%active(k)) then @@ -600,13 +617,12 @@ subroutine update_accum_field_timeavg(this, level, nstep, field) end if end do - if (mod(nstep,this%period) == 0) then - do k = begi,endi - if (this%active(k)) then - this%val(k,level) = this%val(k,level) / this%nsteps(k,level) - end if - end do - end if + do k = begi,endi + effective_nstep = nstep - this%ndays_reset_shifted(k,level) + if (this%active(k) .and. mod(effective_nstep,this%period) == 0) then + this%val(k,level) = this%val(k,level) / this%nsteps(k,level) + end if + end do end subroutine update_accum_field_timeavg @@ -636,6 +652,11 @@ subroutine update_accum_field_runmean(this, level, nstep, field) do k = begi,endi if (this%active(k)) then + if (this%reset(k,level)) then + this%nsteps(k,level) = 0 + this%val(k,level) = this%initval + this%reset(k,level) = .false. + end if kf = k - begi + 1 this%nsteps(k,level) = this%nsteps(k,level) + 1 ! Cap nsteps at 'period' - partly to avoid overflow, but also because it doesn't @@ -649,7 +670,10 @@ subroutine update_accum_field_runmean(this, level, nstep, field) ((accper-1)*this%val(k,level) + field(kf)) / accper end if end do - + + ! SSR 2024-06-05: Note that, unlike other accumulator types, runmean preserves reset requests + ! that occurred when the patch was inactive. + end subroutine update_accum_field_runmean !----------------------------------------------------------------------- @@ -680,9 +704,12 @@ subroutine update_accum_field_runaccum(this, level, nstep, field) do k = begi,endi if (this%active(k)) then kf = k - begi + 1 - if (nint(field(kf)) == -99999) then + if (this%reset(k,level)) then + ! SSR 2024-06-05: Note that, unlike the other accumulator types, runaccum can not + ! reset AND update in the same call of its update_accum_field subroutine. this%val(k,level) = 0._r8 - this%nsteps(k,level) = 0 + this%nsteps(k,level) = this%initval + this%reset(k,level) = .false. else this%val(k,level) = & min(max(this%val(k,level) + field(kf), 0._r8), 99999._r8) @@ -691,9 +718,85 @@ subroutine update_accum_field_runaccum(this, level, nstep, field) end if end do + ! Ignore reset requests that occurred when patch was inactive + this%reset(begi:endi,level) = .false. + end subroutine update_accum_field_runaccum + !----------------------------------------------------------------------- + subroutine markreset_accum_field(name, kf, level) + ! + ! !DESCRIPTION: + ! Mark accumulator values as needing to be reset. Note that resetting happens in + ! update_accum_field subroutines. + ! + ! !ARGUMENTS: + character(len=*), intent(in) :: name ! field name + integer, optional, intent(in) :: kf ! point index to update (in field, not accumulator's val) + integer, optional, intent(in) :: level ! level index to update (in accumulator's val; 1 for a 1-d field) + ! + ! !LOCAL VARIABLES: + integer :: nf ! field index within the accum(nf) array + integer :: begi, endi ! subgrid beginning, ending indices + integer :: k ! index in accumulator's beg1d:end1d + integer :: numlev ! number of levels in this accumulator + + character(len=*), parameter :: subname = 'markreset_accum_field' + !----------------------------------------------------------------------- + + call find_field(field_name=name, caller_name=subname, field_index=nf) + begi = accum(nf)%beg1d + endi = accum(nf)%end1d + numlev = accum(nf)%numlev + + if (present(kf)) then + k = kf + begi - 1 + SHR_ASSERT_FL(k >= begi .and. k <= endi, sourcefile, __LINE__) + if (present(level)) then + ! Reset one level of a single point + accum(nf)%reset(k,level) = .true. + else + ! Reset all levels of a single point + accum(nf)%reset(k,1:numlev) = .true. + end if + else if (present(level)) then + ! Reset one level of all points + accum(nf)%reset(begi:endi,level) = .true. + else + ! Reset all levels of all points + accum(nf)%reset(begi:endi,1:numlev) = .true. + end if + + end subroutine markreset_accum_field + + + !----------------------------------------------------------------------- + function get_accum_reset(name) result(reset) + ! + ! !DESCRIPTION: + ! Get reset array for an accumulator + ! + ! !ARGUMENTS: + character(len=*), intent(in) :: name ! field name + ! + ! !RESULT: + logical, pointer :: reset(:,:) + ! + ! !LOCAL VARIABLES: + integer :: nf ! field index within the accum(nf) array + + character(len=*), parameter :: subname = 'get_reset' + !----------------------------------------------------------------------- + + call find_field(field_name=name, caller_name=subname, field_index=nf) + + allocate(reset(size(accum(nf)%reset, dim=1), size(accum(nf)%reset, dim=2))) + reset(:,:) = accum(nf)%reset + + end function get_accum_reset + + !------------------------------------------------------------------------ subroutine accumulRest( ncid, flag ) ! @@ -786,9 +889,15 @@ subroutine clean_accum_fields if (associated(accum(i)%val)) then deallocate(accum(i)%val) end if + if (associated(accum(i)%reset)) then + deallocate(accum(i)%reset) + end if if (associated(accum(i)%nsteps)) then deallocate(accum(i)%nsteps) end if + if (associated(accum(i)%ndays_reset_shifted)) then + deallocate(accum(i)%ndays_reset_shifted) + end if end do naccflds = 0 diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index c1845d0de6..53b1e55f7c 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -475,7 +475,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro ! When crop calendar streams are being used ! NOTE: This call needs to happen outside loops over nclumps (as streams are not threadsafe) - if (use_cropcal_streams .and. is_beg_curr_year()) then + if (use_crop .and. use_cropcal_streams .and. is_beg_curr_year()) then call cropcal_advance( bounds_proc ) end if @@ -1073,11 +1073,12 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro frictionvel_inst, photosyns_inst, drydepvel_inst) call t_stopf('depvel') - if (use_cropcal_streams .and. is_beg_curr_year()) then + if (use_crop .and. use_cropcal_streams .and. is_beg_curr_year()) then ! ============================================================================ ! Update crop calendars ! ============================================================================ - call cropcal_interp(bounds_clump, filter_inactive_and_active(nc)%num_pcropp, filter_inactive_and_active(nc)%pcropp, crop_inst) + call cropcal_interp(bounds_clump, filter_inactive_and_active(nc)%num_pcropp, & + filter_inactive_and_active(nc)%pcropp, .false., crop_inst) end if ! ============================================================================ @@ -1372,7 +1373,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call atm2lnd_inst%UpdateAccVars(bounds_proc) - call temperature_inst%UpdateAccVars(bounds_proc) + call temperature_inst%UpdateAccVars(bounds_proc, crop_inst) call canopystate_inst%UpdateAccVars(bounds_proc) diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 9bf0cc59a2..0c0570c577 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -659,19 +659,21 @@ subroutine initialize2(ni,nj) end if ! Initialize crop calendars - call t_startf('init_cropcal') - call cropcal_init(bounds_proc) - if (use_cropcal_streams) then - call cropcal_advance( bounds_proc ) - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) - do nc = 1,nclumps - call get_clump_bounds(nc, bounds_clump) - call cropcal_interp(bounds_clump, filter_inactive_and_active(nc)%num_pcropp, & - filter_inactive_and_active(nc)%pcropp, crop_inst) - end do - !$OMP END PARALLEL DO + if (use_crop) then + call t_startf('init_cropcal') + call cropcal_init(bounds_proc) + if (use_cropcal_streams) then + call cropcal_advance( bounds_proc ) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + do nc = 1,nclumps + call get_clump_bounds(nc, bounds_clump) + call cropcal_interp(bounds_clump, filter_inactive_and_active(nc)%num_pcropp, & + filter_inactive_and_active(nc)%pcropp, .true., crop_inst) + end do + !$OMP END PARALLEL DO + end if + call t_stopf('init_cropcal') end if - call t_stopf('init_cropcal') ! Initialize active history fields. ! This is only done if not a restart run. If a restart run, then this diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 0d62407fd9..73d528fbc5 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -402,6 +402,8 @@ module clm_varctl logical, public :: use_cropcal_streams = .false. logical, public :: use_cropcal_rx_swindows = .false. logical, public :: use_cropcal_rx_cultivar_gdds = .false. + logical, public :: adapt_cropcal_rx_cultivar_gdds = .false. + logical, public :: flush_gdd20 = .false. !---------------------------------------------------------- ! biomass heat storage switch diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 95fbd3b4fc..787b827605 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -305,7 +305,7 @@ subroutine control_init(dtime) use_lch4, use_nitrif_denitrif, use_extralakelayers, & use_vichydro, use_cn, use_cndv, use_crop, use_fertilizer, & use_grainproduct, use_snicar_frc, use_vancouver, use_mexicocity, use_noio, & - use_nguardrail, crop_residue_removal_frac + use_nguardrail, crop_residue_removal_frac, flush_gdd20 ! SNICAR namelist /clm_inparm/ & @@ -717,6 +717,7 @@ subroutine control_spmd() call mpi_bcast (use_cndv, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_nguardrail, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_crop, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (flush_gdd20, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_fertilizer, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_grainproduct, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (crop_residue_removal_frac, 1, MPI_REAL8, 0, mpicom, ier) @@ -1000,6 +1001,7 @@ subroutine control_print () write(iulog,*) ' use_cn = ', use_cn write(iulog,*) ' use_cndv = ', use_cndv write(iulog,*) ' use_crop = ', use_crop + write(iulog,*) ' flush_gdd20 = ', flush_gdd20 write(iulog,*) ' use_fertilizer = ', use_fertilizer write(iulog,*) ' use_grainproduct = ', use_grainproduct write(iulog,*) ' crop_residue_removal_frac = ', crop_residue_removal_frac diff --git a/src/main/test/accumul_test/test_accumul.pf b/src/main/test/accumul_test/test_accumul.pf index 15d2e7403e..423a0aea18 100644 --- a/src/main/test/accumul_test/test_accumul.pf +++ b/src/main/test/accumul_test/test_accumul.pf @@ -5,7 +5,7 @@ module test_accumul use funit use accumulMod use unittestSubgridMod - use unittestSimpleSubgridSetupsMod, only : setup_single_veg_patch + use unittestSimpleSubgridSetupsMod, only : setup_single_veg_patch, setup_n_veg_patches use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varcon, only : spval use PatchType, only : patch @@ -94,12 +94,12 @@ contains numlev = nlev, & subgrid_type = 'pft', & init_value = l_init_value, & - type2d = 'irrelevant', & ! type2d just needed for restart + type2d = 'patch', & ! Irrelevant for one-patch fields, but some tests have more than one scale_by_thickness = .false.) end subroutine init_ml_patch_field subroutine update_and_extract_sl_patch_field(this, fieldname, values, val_output, & - pactive, timestep_start) + pactive, timestep_start, reset) ! Calls update_accum_field once for each value in 'values', assuming that the values ! come once per timestep. For the first call, all input values are set equal to ! values(1); for the second call, all input values are set equal to values(2); etc. @@ -122,17 +122,24 @@ contains ! If present, this specifies the starting nstep value. If absent, we start with 1. integer, optional, intent(in) :: timestep_start + ! If present, same size as 'values', indicating whether each should be associated with a reset + logical, optional, intent(in) :: reset(:) + integer :: n_timesteps integer :: timestep integer :: timestep_offset real(r8), pointer :: vals_input(:) real(r8), pointer :: vals_output(:) logical, allocatable :: l_pactive(:) ! local version of pactive + logical, allocatable :: l_reset(:) ! local version of reset n_timesteps = size(values) if (present(pactive)) then @assertEqual(n_timesteps, size(pactive)) end if + if (present(reset)) then + @assertEqual(n_timesteps, size(reset)) + end if allocate(l_pactive(n_timesteps)) if (present(pactive)) then @@ -147,11 +154,21 @@ contains timestep_offset = 0 end if + allocate(l_reset(n_timesteps)) + if (present(reset)) then + l_reset(:) = reset(:) + else + l_reset(:) = .false. + end if + allocate(vals_input(bounds%begp:bounds%endp)) allocate(vals_output(bounds%begp:bounds%endp)) do timestep = 1, n_timesteps vals_input(:) = values(timestep) patch%active(bounds%begp) = l_pactive(timestep) + if (l_reset(timestep)) then + call markreset_accum_field(fieldname) + end if call update_accum_field(fieldname, vals_input, timestep+timestep_offset) end do call extract_accum_field(fieldname, vals_output, n_timesteps+timestep_offset) @@ -257,6 +274,54 @@ contains @assertEqual(expected, val_output, tolerance=tol) end subroutine timeavg_basic + @Test + subroutine timeavg_reset1(this) + ! Test reset of timeavg field 1 day after period start + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: accum_period = 3 + real(r8), parameter :: values(accum_period+1) = [11._r8, 12._r8, 13._r8, 14._r8] + logical, parameter :: reset(accum_period+1) = [.false., .true., .false., .false.] + real(r8) :: val_output + real(r8) :: expected + + ! Setup + call setup_single_veg_patch(pft_type=1) + call this%init_sl_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period) + + ! Exercise + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, reset=reset) + + ! Verify + expected = sum(values(2:4))/accum_period + @assertEqual(expected, val_output, tolerance=tol) + end subroutine timeavg_reset1 + + @Test + subroutine timeavg_reset2(this) + ! Test reset of timeavg field 2 days after period start + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: accum_period = 3 + real(r8), parameter :: values(accum_period+2) = [11._r8, 12._r8, 13._r8, 14._r8, 15._r8] + logical, parameter :: reset(accum_period+2) = [.false., .false., .true., .false., .false.] + real(r8) :: val_output + real(r8) :: expected + + ! Setup + call setup_single_veg_patch(pft_type=1) + call this%init_sl_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period) + + ! Exercise + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, reset=reset) + + ! Verify + expected = sum(values(3:5))/accum_period + @assertEqual(expected, val_output, tolerance=tol) + end subroutine timeavg_reset2 + @Test subroutine timeavg_wrongTime(this) ! Test a timeavg field when it's the wrong time for producing an average @@ -303,6 +368,32 @@ contains @assertEqual(expected, val_output, tolerance=tol) end subroutine timeavg_onlyLatestPeriod + @Test + subroutine timeavg_onlyLatestPeriod_redundantReset(this) + ! Manually requesting a reset when one was going to happen anyway should have no effect. + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: accum_period = 3 + real(r8), parameter :: values(accum_period*2) = & + [11._r8, 12._r8, 13._r8, 21._r8, 22._r8, 23._r8] + logical, parameter :: reset(accum_period*2) = & + [.false., .false., .false., .true., .false., .false.] + real(r8) :: val_output + real(r8) :: expected + + ! Setup + call setup_single_veg_patch(pft_type=1) + call this%init_sl_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period) + + ! Exercise + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, reset=reset) + + ! Verify + expected = sum(values(accum_period+1:2*accum_period))/accum_period + @assertEqual(expected, val_output, tolerance=tol) + end subroutine timeavg_onlyLatestPeriod_redundantReset + @Test subroutine timeavg_newlyActive(this) ! For timeavg: If a point becomes active in the middle of a period, then it should @@ -505,6 +596,54 @@ contains @assertEqual(expected_ts5, val_output, tolerance=tol) end subroutine runmean_afterPeriod + @Test + subroutine runmean_afterPeriod_reset(this) + ! Test runmean accumulation after accum_period is reached, with a reset + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: accum_period = 3 + real(r8), parameter :: values(5) = [11._r8, 22._r8, 43._r8, 110._r8, 17._r8] + logical, parameter :: reset(5) = [.false., .false., .false., .true., .false.] + real(r8) :: val_output + real(r8) :: expected_ts5 + + ! Setup + call setup_single_veg_patch(pft_type=1) + call this%init_sl_patch_field(name=fieldname, accum_type='runmean', & + accum_period = accum_period) + + ! Exercise + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, reset=reset) + + ! Verify + expected_ts5 = (values(4) + values(5)) / 2._r8 + @assertEqual(expected_ts5, val_output, tolerance=tol) + end subroutine runmean_afterPeriod_reset + + @Test + subroutine runmean_afterPeriod_resetWhileInactive(this) + ! Test runmean accumulation after accum_period is reached, with a reset while the patch was inactive. Unlike the other accumulator types, runmean should preserve this reset request and apply it when the patch is active again. This may or may not be the ideal behavior; we can change this if some other + ! behavior would be better in this situation. + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: accum_period = 3 + real(r8), parameter :: values(5) = [11._r8, 22._r8, 43._r8, 110._r8, 17._r8] + logical, parameter :: pactive(5) = [.true., .true., .true., .false., .true.] + logical, parameter :: reset(5) = [.false., .false., .false., .true., .false.] + real(r8) :: val_output + + ! Setup + call setup_single_veg_patch(pft_type=1) + call this%init_sl_patch_field(name=fieldname, accum_type='runmean', & + accum_period = accum_period) + + ! Exercise + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, pactive=pactive, reset=reset) + + ! Verify + @assertEqual(values(5), val_output, tolerance=tol) + end subroutine runmean_afterPeriod_resetWhileInactive + @Test subroutine runmean_newlyActive(this) ! For runmean: If a point recently became active, its running mean should only @@ -599,7 +738,8 @@ contains class(TestAccumul), intent(inout) :: this character(len=*), parameter :: fieldname = 'foo' integer, parameter :: accum_period = 3 ! irrelevant for this type - real(r8), parameter :: values(5) = [11._r8, 12._r8, accumResetVal, 13._r8, 24._r8] + real(r8), parameter :: values(5) = [11._r8, 12._r8, -99999._r8, 13._r8, 24._r8] + logical , parameter :: reset(5) = [.false., .false., .true., .false., .false.] real(r8) :: val_output real(r8) :: expected @@ -609,7 +749,7 @@ contains accum_period = accum_period) ! Exercise - call this%update_and_extract_sl_patch_field(fieldname, values, val_output) + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, reset=reset) ! Verify expected = sum(values(4:5)) @@ -626,7 +766,8 @@ contains class(TestAccumul), intent(inout) :: this character(len=*), parameter :: fieldname = 'foo' integer, parameter :: accum_period = 3 ! irrelevant for this type - real(r8), parameter :: values(5) = [11._r8, accumResetVal, 12._r8, 13._r8, 24._r8] + real(r8), parameter :: values(5) = [11._r8, -99999._r8, 12._r8, 13._r8, 24._r8] + logical , parameter :: reset(5) = [.false., .true., .false., .false., .false.] logical, parameter :: pactive(5) = [.false., .false., .false., .true., .true.] real(r8) :: val_output real(r8) :: expected @@ -649,7 +790,7 @@ contains ! Test runaccum with a point that starts active, becomes inactive, then later becomes ! active again. ! - ! Should ignore values and accumResetVal in the inactive steps. + ! Should ignore values and reset request in the inactive steps. ! ! Also, should continue where it left off - i.e., including the values accumulated ! when it was first active. This may or may not be the ideal behavior; we can change @@ -657,7 +798,8 @@ contains class(TestAccumul), intent(inout) :: this character(len=*), parameter :: fieldname = 'foo' integer, parameter :: accum_period = 3 ! irrelevant for this type - real(r8), parameter :: values(5) = [11._r8, accumResetVal, 12._r8, 17._r8, 24._r8] + real(r8), parameter :: values(5) = [11._r8, -99999._r8, 12._r8, 17._r8, 24._r8] + logical , parameter :: reset(5) = [.false., .true., .false., .false., .false.] logical, parameter :: pactive(5) = [.true., .false., .false., .true., .true.] real(r8) :: val_output real(r8) :: expected @@ -668,7 +810,7 @@ contains accum_period = accum_period) ! Exercise - call this%update_and_extract_sl_patch_field(fieldname, values, val_output, pactive=pactive) + call this%update_and_extract_sl_patch_field(fieldname, values, val_output, pactive=pactive, reset=reset) ! Verify expected = values(1) + values(4) + values(5) @@ -718,4 +860,150 @@ contains @assertEqual(expected2, val_output2, tolerance=tol) end subroutine multipleFields + ! ------------------------------------------------------------------------ + ! Tests of markreset_accum_field() + ! ------------------------------------------------------------------------ + + @Test + subroutine markreset_nopoints_nolevels(this) + ! Make sure that NOT calling markreset_accum_field() means no point-levels are marked for reset. + ! + ! Note that type of accumulator and values don't matter. + + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: npatches = 2 + integer, parameter :: nlevels = 5 + real(r8), parameter :: pwtcol(npatches) = [0.5, 0.5] + integer, parameter :: pft_type(npatches) = [1, 2] + integer, parameter :: accum_period = 3 + integer, parameter :: expected(npatches, nlevels) = transpose(reshape([ & + 0, 0, 0, 0, 0, & + 0, 0, 0, 0, 0 & + ], [nlevels, npatches])) + logical :: reset_output(npatches, nlevels) + integer :: reset_output_int(npatches, nlevels) + integer :: l, p + + ! Setup + call setup_n_veg_patches(pwtcol, pft_type) + call this%init_ml_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period, nlev=nlevels) + + ! Exercise + reset_output = get_accum_reset(fieldname) + + ! Verify + reset_output_int = merge(1, 0, reset_output) + @assertEqual(expected, reset_output_int) + end subroutine markreset_nopoints_nolevels + + @Test + subroutine markreset_1point_nolevels(this) + ! Make sure that calling markreset_accum_field() with kf but no level works right (marks + ! all levels of that patch for reset). + ! Note that type of accumulator and values don't matter. + + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: npatches = 2 + integer, parameter :: nlevels = 5 + real(r8), parameter :: pwtcol(npatches) = [0.5, 0.5] + integer, parameter :: pft_type(npatches) = [1, 2] + integer, parameter :: accum_period = 3 + integer, parameter :: expected(npatches, nlevels) = transpose(reshape([ & + 1, 1, 1, 1, 1, & + 0, 0, 0, 0, 0 & + ], [nlevels, npatches])) + logical :: reset_output(npatches, nlevels) + integer :: reset_output_int(npatches, nlevels) + integer :: l, p + + ! Setup + call setup_n_veg_patches(pwtcol, pft_type) + call this%init_ml_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period, nlev=nlevels) + + ! Exercise + call markreset_accum_field(fieldname, kf=1) + reset_output = get_accum_reset(fieldname) + + ! Verify + reset_output_int = merge(1, 0, reset_output) + @assertEqual(expected, reset_output_int) + end subroutine markreset_1point_nolevels + + @Test + subroutine markreset_allpoints_1level(this) + ! Make sure that calling markreset_accum_field() with level but no kf works right (marks that + ! level for reset for all points). + ! + ! Note that type of accumulator and values don't matter. + + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: npatches = 2 + integer, parameter :: nlevels = 5 + real(r8), parameter :: pwtcol(npatches) = [0.5, 0.5] + integer, parameter :: pft_type(npatches) = [1, 2] + integer, parameter :: accum_period = 3 + integer, parameter :: expected(npatches, nlevels) = transpose(reshape([ & + 0, 0, 1, 0, 0, & + 0, 0, 1, 0, 0 & + ], [nlevels, npatches])) + logical :: reset_output(npatches, nlevels) + integer :: reset_output_int(npatches, nlevels) + integer :: l, p + + ! Setup + call setup_n_veg_patches(pwtcol, pft_type) + call this%init_ml_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period, nlev=nlevels) + + ! Exercise + call markreset_accum_field(fieldname, level=3) + reset_output = get_accum_reset(fieldname) + + ! Verify + reset_output_int = merge(1, 0, reset_output) + @assertEqual(expected, reset_output_int) +end subroutine markreset_allpoints_1level + +@Test +subroutine markreset_allpoints_alllevels(this) + ! Make sure that calling markreset_accum_field() with neither kf nor level works right (marks + ! all for reset). + ! + ! Note that type of accumulator and values don't matter. + + class(TestAccumul), intent(inout) :: this + character(len=*), parameter :: fieldname = 'foo' + integer, parameter :: npatches = 2 + integer, parameter :: nlevels = 5 + real(r8), parameter :: pwtcol(npatches) = [0.5, 0.5] + integer, parameter :: pft_type(npatches) = [1, 2] + integer, parameter :: accum_period = 3 + integer, parameter :: expected(npatches, nlevels) = transpose(reshape([ & + 1, 1, 1, 1, 1, & + 1, 1, 1, 1, 1 & + ], [nlevels, npatches])) + logical :: reset_output(npatches, nlevels) + integer :: reset_output_int(npatches, nlevels) + integer :: l, p + + ! Setup + call setup_n_veg_patches(pwtcol, pft_type) + call this%init_ml_patch_field(name=fieldname, accum_type='timeavg', & + accum_period = accum_period, nlev=nlevels) + + ! Exercise + call markreset_accum_field(fieldname) + reset_output = get_accum_reset(fieldname) + + ! Verify + reset_output_int = merge(1, 0, reset_output) + @assertEqual(expected, reset_output_int) +end subroutine markreset_allpoints_alllevels + + end module test_accumul diff --git a/tools/contrib/tweak_latlons.py b/tools/contrib/tweak_latlons.py new file mode 100644 index 0000000000..2bae06d229 --- /dev/null +++ b/tools/contrib/tweak_latlons.py @@ -0,0 +1,269 @@ +""" +'Tweak' the latitude and longitude coordinates to avoid ambiguous nearest neighbors +""" +import os +import sys +import contextlib +import argparse +import numpy as np +import xarray as xr +from netCDF4 import Dataset # pylint: disable=no-name-in-module + +# -- add python/ctsm to path (needed if we want to run this stand-alone) +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python") +sys.path.insert(1, _CTSM_PYTHON) +# pylint: disable=wrong-import-position +from ctsm.mesh_maker import main as mesh_maker + +COORD_LIST = ["lat", "lon"] +COORD_DATATYPE = np.float64 + +def get_tweak(ds_in, coord_str, init_tweak): + """ + Get the tweak that will be applied to all datasets' lat/lon coordinates + """ + da = ds_in[coord_str] + coord2_orig = da.values.astype(COORD_DATATYPE) + coord2 = coord2_orig + tweak = init_tweak + coord2 += tweak + + # This is necessary if precision is lower than float64 + max_tweak = 1e-2 + while np.any(coord2 == da.values): + tweak *= 10 + if tweak > max_tweak: + raise RuntimeError(f"Tweaking by +{max_tweak} failed to 'take'") + coord2 = coord2_orig + coord2 += tweak + return tweak + +def apply_tweak(ds_in, coord_str, tweak): + # Apply tweak + da = ds_in[coord_str] + coord2 = da.values.astype(COORD_DATATYPE) + coord2 += tweak + if np.any(coord2 == da.values): + raise RuntimeError('Tweak didn''t "take"') + coord_tweak = np.full_like(coord2, tweak) + + # Ensure that no value is above maximum in input data. This is needed for mesh_maker to work. + max_coord = np.max(da.values) + where_toohigh = np.where(coord2 > max_coord) + Ntoohigh = len(where_toohigh[0]) + if Ntoohigh != 1: + raise RuntimeError( + f"Expected 1 coordinate value too high; got {Ntoohigh}" + ) + coord2[where_toohigh] = max_coord + coord_tweak[where_toohigh] = max_coord + + # Convert to DataArray + new_coords_dict = {coord_str: coord2} + da2 = xr.DataArray( + data=coord2, + coords=new_coords_dict, + dims=da.dims, + attrs=da.attrs, + ) + + # Replace coordinate in dataset + ds_in[coord_str] = da2 + + # Add a variable with the amount of the tweak + tweak_attrs = {} + if "standard_name" in da.attrs: + coord_name = da.attrs["standard_name"] + elif "long_name" in da.attrs: + coord_name = da.attrs["long_name"].replace("coordinate_", "") + else: + coord_name = coord_str + tweak_attrs["standard_name"] = coord_name + "_tweak" + tweak_attrs[ + "long_name" + ] = f"Amount {coord_name} was shifted to avoid ambiguous nearest neighbors" + if "units" in da.attrs: + tweak_attrs["units"] = da.attrs["units"] + da_tweak = xr.DataArray( + data=coord_tweak, + coords=new_coords_dict, + dims=da.dims, + attrs=tweak_attrs, + ) + tweak_name = coord_str + "_tweak" + ds_in[tweak_name] = da_tweak + + return ds_in + +def check(ds, f0_base, ds2, f_base, var): + if not np.array_equal(ds[var].values, ds2[var].values): + if not np.array_equal(ds[var].shape, ds2[var].shape): + msg = f"{var} shapes differ b/w {f0_base} ({ds[var].shape}) and {f_base} ({ds2[var].shape})" + raise RuntimeError(msg) + max_diff = np.max(np.abs(ds[var].values - ds2[var].values)) + msg = f"{var}s differ between {f0_base} and {f_base}; max = {max_diff}" + type0 = type(ds[var].values[0]) + type2 = type(ds2[var].values[0]) + if type0 != type2: + msg += f"\nTypes also differ: {type0} vs. {type2}" + raise RuntimeError(msg) + +@contextlib.contextmanager +def redirect_argv(arglist): + """ + Preserve actual arg list while giving a new one to mesh_maker + """ + argv_tmp = sys.argv[:] + sys.argv = arglist + yield + sys.argv = argv_tmp + +def main(input_files, mesh_file_in, output_files): + """ + Apply tweak to all files + """ + + # Set up + tweak_dict = {} + for coord in COORD_LIST: + tweak_dict[coord] = -np.inf + mesh_file_out = output_files[-1] + output_files = output_files[:-1] + + # Get tweaks + for file_in in input_files: + ds = xr.open_dataset(file_in) + for coord in COORD_LIST: + this_tweak = get_tweak(ds, coord, init_tweak=1e-6) + if this_tweak > tweak_dict[coord]: + tweak_dict[coord] = this_tweak + for coord in COORD_LIST: + print(f"Tweaking {coord} by {tweak_dict[coord]}") + print(" ") + + # Apply tweaks + for i, file_in in enumerate(input_files): + ds = xr.open_dataset(file_in) + + for coord in COORD_LIST: + ds = apply_tweak(ds, coord, tweak_dict[coord]) + + # Set up for save + file_out = output_files[i] + with Dataset(file_in, "r") as netcdf_file: + netcdf_format = netcdf_file.data_model + + # Make output dir, if needed + output_dir = os.path.dirname(file_out) + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + # Save + print(f"Saving {file_out}") + ds.to_netcdf(file_out, format=netcdf_format) + print("Done") + + + # Ensure all files got the same tweaks + ds = xr.open_dataset(output_files[0]) + f0_base = os.path.basename(output_files[0]) + for file_out in output_files[1:]: + ds2 = xr.open_dataset(file_out) + f_base = os.path.basename(file_out) + for coord in COORD_LIST: + check(ds, f0_base, ds2, f_base, coord) + check(ds, f0_base, ds2, f_base, coord + "_tweak") + + + # Save new mesh file + mesh_maker_args = [ + "mesh_maker", + "--input", + output_files[0], + "--output", + mesh_file_out, + "--lat", + "lat", + "--lon", + "lon", + "--overwrite", + ] + print(f"Saving {mesh_file_out}...") + with redirect_argv(mesh_maker_args): + mesh_maker() + + # Change format, if needed + with Dataset(mesh_file_in, "r") as netcdf_file: + netcdf_format_in = netcdf_file.data_model + with Dataset(mesh_file_out, "r") as netcdf_file: + netcdf_format_out = netcdf_file.data_model + if netcdf_format_in != netcdf_format_out: + mesh_file_out_tmp = mesh_file_out + ".tmp" + os.rename(mesh_file_out, mesh_file_out_tmp) + ds = xr.open_dataset(mesh_file_out_tmp) + ds.to_netcdf(mesh_file_out, format=netcdf_format_in) + os.remove(mesh_file_out_tmp) + + print("Done") + + + + +if __name__ == "__main__": + ############################### + ### Process input arguments ### + ############################### + parser = argparse.ArgumentParser( + description="'Tweak' the latitude and longitude coordinates to avoid ambiguous nearest neighbors", + ) + + # Required + parser.add_argument( + "-i", + "--input-files", + help="Comma-separated stream files whose coordinates need tweaking", + required=True, + ) + parser.add_argument( + "-m", + "--mesh-file", + help="Mesh file associated with input files", + required=True, + ) + + # Optional + parser.add_argument( + "--overwrite", + help="Overwrite any existing output files", + action="store_true", + default=False, + ) + default_output_dir = os.getcwd() + parser.add_argument( + "-o", + "--output-dir", + help=f"Directory where output files should be saved. Default is current working directory: {default_output_dir}", + default=default_output_dir, + ) + + # Get arguments + args = parser.parse_args(sys.argv[1:]) + + # Check/process input and output files + _input_files = args.input_files.split(",") + _output_files = [] + for file in _input_files + [args.mesh_file]: + if not os.path.exists(file): + raise FileNotFoundError(f"File not found: {file}") + + filename, ext = os.path.splitext(os.path.basename(file)) + output_file = os.path.join( + args.output_dir, filename + ".tweaked_latlons" + ext + ) + if os.path.exists(output_file) and not args.overwrite: + raise FileExistsError( + f"Output file exists but --overwrite not specified: {output_file}" + ) + _output_files.append(output_file) + + main(_input_files, args.mesh_file, _output_files) diff --git a/tools/crop_calendars/generate_gdd20_baseline b/tools/crop_calendars/generate_gdd20_baseline new file mode 100755 index 0000000000..a0238c8d0f --- /dev/null +++ b/tools/crop_calendars/generate_gdd20_baseline @@ -0,0 +1,19 @@ +#!/usr/bin/env python3 +""" +For description and instructions, please see README. +""" + +import os +import sys + +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), + os.pardir, + os.pardir, + 'python') +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm.crop_calendars.generate_gdd20_baseline import main + +if __name__ == "__main__": + main() + diff --git a/tools/mksurfdata_esmf/Makefile b/tools/mksurfdata_esmf/Makefile index d8bacdc5dd..c344843d06 100644 --- a/tools/mksurfdata_esmf/Makefile +++ b/tools/mksurfdata_esmf/Makefile @@ -54,7 +54,10 @@ SUBSETDATA_POINT_URBAN = $(SUBSETDATA_POINT) --include-nonveg # Subset data sites... SUBSETDATA_1X1_BRAZIL := --lat -7 --lon -55 --site 1x1_brazil SUBSETDATA_1X1_NUMAIA := --lat 40.6878 --lon 267.0228 --site 1x1_numaIA -SUBSETDATA_1X1_SMALL := --lat 40.6878 --lon 267.0228 --site 1x1_smallvilleIA \ +SUBSETDATA_1X1_SMALL_IA := --lat 40.6878 --lon 267.0228 --site 1x1_smallvilleIA \ + --dompft 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 \ + --pctpft 6.5 1.5 1.6 1.7 1.8 1.9 1.5 1.6 1.7 1.8 1.9 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 +SUBSETDATA_1X1_SMALL_BR := --lat -12.9952 --lon 305.3233 --site 1x1_cidadinhoBR \ --dompft 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 \ --pctpft 6.5 1.5 1.6 1.7 1.8 1.9 1.5 1.6 1.7 1.8 1.9 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 # NOTE: The 1850 smallvilleIA site is constructed to start with 100% natural vegetation, so we can test transition to crops @@ -112,6 +115,7 @@ all-subset : \ 1x1-smallville-present \ 1x1-smallville-1850 \ 1x1-smallville-transient \ + 1x1-cidadinho-present \ urban DEBUG: @@ -239,7 +243,7 @@ crop-global-hist-ne30 : FORCE $(SUBSETDATA_POINT_ALLLU) --create-surface $(SUBSETDATA_1X1_NUMAIA) 1x1-smallville-present : FORCE - $(SUBSETDATA_POINT) --create-surface $(SUBSETDATA_1X1_SMALL) + $(SUBSETDATA_POINT) --create-surface $(SUBSETDATA_1X1_SMALL_IA) # Note that the smallville 1850 dataset is entirely natural vegetation. This # facilitates testing a transient case that starts with no crop, and then later @@ -254,6 +258,9 @@ crop-global-hist-ne30 : FORCE $(SUBSETDATA_POINT) --create-landuse $(SUBSETDATA_1X1_SMALLTRANSIENT) ../modify_input_files/modify_smallville.sh +1x1-cidadinho-present : FORCE + $(SUBSETDATA_POINT) --create-surface $(SUBSETDATA_1X1_SMALL_BR) + # # Crop with future scenarios #