diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 2620839223..31f63ef8ae 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -47,3 +47,5 @@ aa04d1f7d86cc2503b98b7e2b2d84dbfff6c316b 6c6f57e948bfa31e60b383536cc21663fedb8b70 9660667b1267dcd4150889f5f39db540158be74a 665cf86102e09b4c4c5a140700676dca23bc55a9 +045d90f1d80f713eb3ae0ac58f6c2352937f1eb0 +753fda3ff0147837231a73c9c728dd9ce47b5997 diff --git a/Externals.cfg b/Externals.cfg index 185f412cab..49cfd25962 100644 --- a/Externals.cfg +++ b/Externals.cfg @@ -34,9 +34,9 @@ hash = 34723c2 required = True [ccs_config] -tag = ccs_config_cesm0.0.92 +hash = 75847d58bac1149adcb8d9b40181df631c273e19 protocol = git -repo_url = https://github.com/ESMCI/ccs_config_cesm.git +repo_url = https://github.com/samsrabin/ccs_config_cesm.git local_path = ccs_config required = True diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 26725b7d96..5202717d74 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1842,9 +1842,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."); } } @@ -4175,8 +4176,11 @@ sub setup_logic_cropcal_streams { # 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') ; @@ -4185,12 +4189,27 @@ sub setup_logic_cropcal_streams { $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." ); + } + } + # 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'); - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_cropcal'); if ( &value_is_true($cropcals_rx_adapt) ) { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldFileName_gdd20_baseline'); } @@ -4202,7 +4221,7 @@ sub setup_logic_cropcal_streams { 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 ( !&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) or !&string_is_undef_or_empty($mesh_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 diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index da0a523adc..70ff048aa7 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -1406,6 +1406,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 @@ -1477,11 +1479,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 + + + @@ -1689,6 +1694,8 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c .false. .false. +.false. +.false. 2000 2000 2000 @@ -1697,15 +1704,18 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c 2000 -lnd/clm2/cropdata/calendars/processed/swindow_starts_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.nc -lnd/clm2/cropdata/calendars/processed/swindow_ends_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.nc -lnd/clm2/cropdata/calendars/processed/gdds_20230829_161011.nc -share/meshes/360x720_120830_ESMFmesh_c20210507_cdf5.nc -lnd/clm2/cropdata/calendars/processed/swindow_starts_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.nc -lnd/clm2/cropdata/calendars/processed/swindow_ends_ggcmi_crop_calendar_phase3_v1.01.2000-2000.20231005_145103.nc -lnd/clm2/cropdata/calendars/processed/gdds_20230829_161011.nc -lnd/clm2/testdata/gdd20baseline.tmp_dontupload.nc -share/meshes/360x720_120830_ESMFmesh_c20210507_cdf5.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/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 50e645519c..a24f59339a 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1838,6 +1838,31 @@ Filename of input stream data for cultivar growing degree-day targets 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 b0ddb1e448..5c8aa3468d 100755 --- a/bld/unit_testers/build-namelist_test.pl +++ b/bld/unit_testers/build-namelist_test.pl @@ -1590,21 +1590,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 75fff8a0e0..965398821f 100644 --- a/cime_config/SystemTests/rxcropmaturity.py +++ b/cime_config/SystemTests/rxcropmaturity.py @@ -20,16 +20,22 @@ from CIME.SystemTests.system_tests_common import SystemTestsCommon from CIME.XML.standard_module_setup import * from CIME.SystemTests.test_utils.user_nl_utils import append_to_user_nl_files +from CIME.case import Case import shutil, glob logger = logging.getLogger(__name__) -class RXCROPMATURITY(SystemTestsCommon): +class RXCROPMATURITYSHARED(SystemTestsCommon): def __init__(self, case): # initialize an object interface to the SMS system test SystemTestsCommon.__init__(self, case) + # Is this a real RXCROPMATURITY test or not? + casebaseid = self._case.get_value("CASEBASEID") + full_test = "RXCROPMATURITY_" in casebaseid + skipgen_test = "RXCROPMATURITYSKIPGEN_" in casebaseid + # Ensure run length is at least 5 years. Minimum to produce one complete growing season (i.e., two complete calendar years) actually 4 years, but that only gets you 1 season usable for GDD generation, so you can't check for season-to-season consistency. stop_n = self._case.get_value("STOP_N") stop_option = self._case.get_value("STOP_OPTION") @@ -56,11 +62,17 @@ def __init__(self, case): f"STOP_OPTION ({stop_option_orig}) must be nsecond(s), nminute(s), " + "nhour(s), nday(s), nmonth(s), or nyear(s)" ) - elif stop_n < 5: + elif full_test and stop_n < 5: error_message = ( "RXCROPMATURITY must be run for at least 5 years; you requested " + f"{stop_n_orig} {stop_option_orig[1:]}" ) + elif skipgen_test and stop_n < 3: + # First year is discarded because crops are already in the ground at restart, and those aren't affected by the new crop calendar inputs. The second year is useable, but we need a third year so that all crops planted in the second year have a chance to finish. + error_message = ( + "RXCROPMATURITYSKIPGEN (both-forced part) must be run for at least 3 years; you requested " + + f"{stop_n_orig} {stop_option_orig[1:]}" + ) if error_message is not None: logger.error(error_message) raise RuntimeError(error_message) @@ -69,7 +81,6 @@ def __init__(self, case): self._run_Nyears = int(stop_n) # Only allow RXCROPMATURITY to be called with test cropMonthOutput - casebaseid = self._case.get_value("CASEBASEID") if casebaseid.split("-")[-1] != "cropMonthOutput": error_message = ( "Only call RXCROPMATURITY with test cropMonthOutput " @@ -81,10 +92,16 @@ def __init__(self, case): # Get files with prescribed sowing and harvest dates self._get_rx_dates() + # Get cultivar maturity requirement file to fall back on if not generating it here + self._gdds_file = None + self._fallback_gdds_file = os.path.join( + os.path.dirname(self._sdatefile), "gdds_20230829_161011.nc" + ) + # Which conda environment should we use? self._get_conda_env() - def run_phase(self): + def _run_phase(self, skip_gen=False): # Modeling this after the SSP test, we create a clone to be the case whose outputs we don't # want to be saved as baseline. @@ -113,6 +130,7 @@ def run_phase(self): 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.", " ", @@ -133,6 +151,12 @@ def run_phase(self): # Download files from the server, if needed case_gddgen.check_all_input_data() + # Copy needed file from original to gddgen directory + shutil.copyfile( + os.path.join(caseroot, ".env_mach_specific.sh"), + os.path.join(self._path_gddgen, ".env_mach_specific.sh"), + ) + # Make custom version of surface file logger.info("RXCROPMATURITY log: run fsurdat_modifier") self._run_fsurdat_modifier() @@ -146,9 +170,19 @@ def run_phase(self): # "No history files expected, set suffix=None to avoid compare error" # We *do* expect history files here, but anyway. This works. self._skip_pnl = False - self.run_indv(suffix=None, st_archive=True) - self._run_generate_gdds(case_gddgen) + # If not generating GDDs, only run a few days of this. + if skip_gen: + with Case(self._path_gddgen, read_only=False) as case: + case.set_value("STOP_N", 5) + case.set_value("STOP_OPTION", "ndays") + + self.run_indv(suffix=None, st_archive=True) + if skip_gen: + # Interpolate an existing GDD file. Needed to check obedience to GDD inputs. + self._run_interpolate_gdds() + else: + self._run_generate_gdds(case_gddgen) # ------------------------------------------------------------------- # (3) Set up and perform Prescribed Calendars run @@ -174,7 +208,7 @@ def run_phase(self): # (4) Check Prescribed Calendars run # ------------------------------------------------------------------- logger.info("RXCROPMATURITY log: output check: Prescribed Calendars") - self._run_check_rxboth_run() + self._run_check_rxboth_run(skip_gen) # Get sowing and harvest dates for this resolution. def _get_rx_dates(self): @@ -331,11 +365,16 @@ def _create_config_file_evenlysplitcrop(self): cfg_out.write("PCT_OCEAN = 0.0\n") cfg_out.write("PCT_URBAN = 0.0 0.0 0.0\n") - def _run_check_rxboth_run(self): + def _run_check_rxboth_run(self, skip_gen): output_dir = os.path.join(self._get_caseroot(), "run") - first_usable_year = self._run_startyear + 2 - last_usable_year = self._run_startyear + self._run_Nyears - 2 + + if skip_gen: + first_usable_year = self._run_startyear + 1 + last_usable_year = first_usable_year + else: + first_usable_year = self._run_startyear + 2 + last_usable_year = self._run_startyear + self._run_Nyears - 2 tool_path = os.path.join( self._ctsm_root, "python", "ctsm", "crop_calendars", "check_rxboth_run.py" @@ -357,6 +396,7 @@ def _run_check_rxboth_run(self): def _modify_user_nl_allruns(self): nl_additions = [ + "cropcals_rx = .true.", "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), @@ -398,7 +438,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( @@ -416,6 +456,30 @@ def _run_generate_gdds(self, case_gddgen): raise RuntimeError(error_message) self._gdds_file = generated_gdd_files[0] + def _run_interpolate_gdds(self): + # Save where? + self._gdds_file = os.path.join(self._get_caseroot(), "interpolated_gdds.nc") + + # It'd be much nicer to call interpolate_gdds.main(), but I can't import interpolate_gdds. + tool_path = os.path.join( + self._ctsm_root, "python", "ctsm", "crop_calendars", "interpolate_gdds.py" + ) + command = " ".join( + [ + f"python3 {tool_path}", + f"--input-file {self._fallback_gdds_file}", + f"--target-file {self._sdatefile}", + f"--output-file {self._gdds_file}", + "--overwrite", + ] + ) + stu.run_python_script( + self._get_caseroot(), + self._this_conda_env, + command, + tool_path, + ) + def _get_conda_env(self): conda_setup_commands = stu.cmds_to_setup_conda(self._get_caseroot()) @@ -442,3 +506,8 @@ def _get_flanduse_timeseries_in(self, case): if flanduse_timeseries_in: self._flanduse_timeseries_in = flanduse_timeseries_in.group(1) break + + +class RXCROPMATURITY(RXCROPMATURITYSHARED): + def run_phase(self): + self._run_phase() diff --git a/cime_config/SystemTests/rxcropmaturityskipgen.py b/cime_config/SystemTests/rxcropmaturityskipgen.py new file mode 100644 index 0000000000..409f2b9847 --- /dev/null +++ b/cime_config/SystemTests/rxcropmaturityskipgen.py @@ -0,0 +1,6 @@ +from rxcropmaturity import RXCROPMATURITYSHARED + + +class RXCROPMATURITYSKIPGEN(RXCROPMATURITYSHARED): + def run_phase(self): + self._run_phase(skip_gen=True) diff --git a/cime_config/config_tests.xml b/cime_config/config_tests.xml index c0b6afed9d..5c5bdc2892 100644 --- a/cime_config/config_tests.xml +++ b/cime_config/config_tests.xml @@ -133,6 +133,16 @@ This defines various CTSM-specific system tests $STOP_N + + As RXCROPMATURITY but don't actually generate GDDs. Allows short testing with existing GDD inputs. + 1 + FALSE + FALSE + never + $STOP_OPTION + $STOP_N + + - - - FAIL - #2460 - - - FAIL @@ -177,6 +170,13 @@ + + + FAIL + #2659 + + + diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index bf763c4775..5d3fbeb11a 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -3660,7 +3660,7 @@ - + @@ -3670,6 +3670,29 @@ + + + + + + + + + + + + + + + + + + + + + + + @@ -3729,13 +3752,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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..02ec13743f --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/GddGen/include_user_mods @@ -0,0 +1 @@ +../cropMonthOutput 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..af5fe8591e --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/include_user_mods @@ -0,0 +1 @@ +../RxCropCalsAdapt 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..dd4ac3117c --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/RxCropCalsAdaptGGCMI/user_nl_clm @@ -0,0 +1,5 @@ + +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/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/sowingWindows/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/sowingWindows/user_nl_clm deleted file mode 100644 index 545e5f6ebd..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_swindows = 2000 -stream_year_last_cropcal_swindows = 2000 diff --git a/python/ctsm/crop_calendars/check_constant_vars.py b/python/ctsm/crop_calendars/check_constant_vars.py index aa25a412fe..be7718256a 100644 --- a/python/ctsm/crop_calendars/check_constant_vars.py +++ b/python/ctsm/crop_calendars/check_constant_vars.py @@ -382,4 +382,4 @@ def check_constant_vars( raise RuntimeError("Stopping due to failed check_constant_vars().") bad_patches = np.unique(bad_patches) - return [int(p) for p in bad_patches] + return [int(p) for p in bad_patches], any_bad diff --git a/python/ctsm/crop_calendars/check_rx_obeyed.py b/python/ctsm/crop_calendars/check_rx_obeyed.py index 99b8d80bde..383ebbf1e0 100644 --- a/python/ctsm/crop_calendars/check_rx_obeyed.py +++ b/python/ctsm/crop_calendars/check_rx_obeyed.py @@ -104,6 +104,27 @@ def get_extreme_info(diff_array, rx_array, mxn, dims, gs_da, patches1d_lon, patc return round(themxn, 3), round(this_lon, 3), round(this_lat, 3), this_gs, round(this_rx) +def summarize_results(which_ds, output_var, verbose, all_ok, gdd_tolerance, diffs_eg_txt): + """ + Summarize results + """ + bad = True + if all_ok == 2: + bad = False + print(f"✅ {which_ds}: Prescribed {output_var} always obeyed") + elif all_ok == 1: + bad = False + # print(f"🟨 {which_ds}: Prescribed {output_var} *not* always obeyed, but acceptable:") + # for x in diff_str_list: print(x) + print( + f"🟨 {which_ds}: Prescribed {output_var} *not* always obeyed, but acceptable (diffs <= " + + f"{gdd_tolerance})" + ) + elif not verbose: + print(f"❌ {which_ds}: Prescribed {output_var} *not* always obeyed. E.g., {diffs_eg_txt}") + return bad + + def check_rx_obeyed( vegtype_list, rx_ds, dates_ds, which_ds, output_var, gdd_min=None, verbose=False ): @@ -114,6 +135,7 @@ def check_rx_obeyed( dates_ds, which_ds, output_var, verbose ) + diffs_eg_txt = None for vegtype_str in vegtype_list: thisveg_patches = np.where(dates_ds.patches1d_itype_veg_str == vegtype_str)[0] if thisveg_patches.size == 0: @@ -203,14 +225,6 @@ def check_rx_obeyed( else: break - if all_ok == 2: - print(f"✅ {which_ds}: Prescribed {output_var} always obeyed") - elif all_ok == 1: - # print(f"🟨 {which_ds}: Prescribed {output_var} *not* always obeyed, but acceptable:") - # for x in diff_str_list: print(x) - print( - f"🟨 {which_ds}: Prescribed {output_var} *not* always obeyed, but acceptable (diffs <= " - + f"{gdd_tolerance})" - ) - elif not verbose: - print(f"❌ {which_ds}: Prescribed {output_var} *not* always obeyed. E.g., {diffs_eg_txt}") + bad = summarize_results(which_ds, output_var, verbose, all_ok, gdd_tolerance, diffs_eg_txt) + + return bad diff --git a/python/ctsm/crop_calendars/check_rxboth_run.py b/python/ctsm/crop_calendars/check_rxboth_run.py index ae4decde30..fa4affd220 100644 --- a/python/ctsm/crop_calendars/check_rxboth_run.py +++ b/python/ctsm/crop_calendars/check_rxboth_run.py @@ -75,6 +75,8 @@ def main(argv): "HARVEST_REASON_PERHARV", ] + any_bad = False + annual_outfiles = glob.glob(os.path.join(args.directory, "*.clm2.h1.*.nc")) # These should be constant in a Prescribed Calendars (rxboth) run, as long as the inputs were @@ -85,13 +87,19 @@ def main(argv): "rx_gdds_file": args.rx_gdds_file, } - case["ds"] = cc.import_output( + case["ds"], any_bad_import_output = cc.import_output( annual_outfiles, my_vars=my_vars, year_1=args.first_usable_year, year_n=args.last_usable_year, + throw_errors=False, + ) + any_bad = any_bad or any_bad_import_output + + _, any_bad_check_const_vars = check_constant_vars( + case["ds"], case, ignore_nan=True, verbose=True, throw_error=True ) - check_constant_vars(case["ds"], case, ignore_nan=True, verbose=True, throw_error=True) + any_bad = any_bad or any_bad_check_const_vars # Import GGCMI sowing and harvest dates, and check sims casename = "Prescribed Calendars" @@ -128,15 +136,16 @@ def main(argv): # Check if case["rx_sdates_file"]: - check_rx_obeyed( + sdate_not_obeyed = check_rx_obeyed( case["ds"].vegtype_str.values, case["rx_sdates_ds"].isel(time=0), case["ds"], casename, "SDATES", ) + any_bad = any_bad or sdate_not_obeyed if case["rx_gdds_file"]: - check_rx_obeyed( + gdds_not_obeyed = check_rx_obeyed( case["ds"].vegtype_str.values, case["rx_gdds_ds"].isel(time=0), case["ds"], @@ -144,6 +153,19 @@ def main(argv): "GDDHARV", gdd_min=gdd_min, ) + any_bad = any_bad or gdds_not_obeyed + + if any_bad: + 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 3fe6942f94..719d352665 100644 --- a/python/ctsm/crop_calendars/cropcal_module.py +++ b/python/ctsm/crop_calendars/cropcal_module.py @@ -13,6 +13,8 @@ from ctsm.crop_calendars.cropcal_constants import DEFAULT_GDD_MIN from ctsm.crop_calendars.import_ds import import_ds +MISSING_RX_GDD_VAL = -1 + def check_and_trim_years(year_1, year_n, ds_in): """ @@ -168,8 +170,10 @@ def check_v0_le_v1(this_ds, var_list, msg_txt=" ", both_nan_ok=False, throw_erro if both_nan_ok: gdd_lt_hui = gdd_lt_hui | (np.isnan(this_ds[var0]) & np.isnan(this_ds[var1])) if np.all(gdd_lt_hui): + any_bad = False print(f"✅{msg_txt}{var0} always <= {var1}") else: + any_bad = True msg = f"❌{msg_txt}{var0} *not* always <= {var1}" gdd_lt_hui_vals = gdd_lt_hui.values patch_index = np.where(~gdd_lt_hui_vals)[0][0] @@ -185,6 +189,7 @@ def check_v0_le_v1(this_ds, var_list, msg_txt=" ", both_nan_ok=False, throw_erro print(msg) else: raise RuntimeError(msg) + return any_bad def get_gs_len_da(this_da): @@ -231,6 +236,13 @@ def import_max_gs_length(paramfile_dir, my_clm_ver, my_clm_subver): return mxmat_dict +def unexpected_negative_rx_gdd(data_array): + """ + Return True if there's a negative value not matching the designated missing value + """ + return np.any((data_array.values < 0) & (data_array.values != MISSING_RX_GDD_VAL)) + + def import_rx_dates(var_prefix, date_infile, dates_ds, set_neg1_to_nan=True): """ Import prescribed sowing/harvest dates @@ -260,19 +272,19 @@ def import_rx_dates(var_prefix, date_infile, dates_ds, set_neg1_to_nan=True): v_new = var.replace(var_prefix, "gs") this_ds = this_ds.rename({var: v_new}) - # Set -1 prescribed GDD values to NaN. Only warn the first time. + # Set GDD values matching MISSING_RX_GDD_VAL to NaN. Only warn the first time. if ( set_neg1_to_nan and var_prefix == "gdd" and v_new != var and np.any(this_ds[v_new].values < 0) ): - if np.any((this_ds[v_new].values < 0) & (this_ds[v_new].values != -1)): + if unexpected_negative_rx_gdd(this_ds[v_new]): raise RuntimeError(f"Unexpected negative value in {var}") if not did_warn: - print("Setting -1 rx GDD values to NaN") + print(f"Setting {MISSING_RX_GDD_VAL} rx GDD values to NaN") did_warn = True - this_ds[v_new] = this_ds[v_new].where(this_ds[v_new] != -1) + this_ds[v_new] = this_ds[v_new].where(this_ds[v_new] != MISSING_RX_GDD_VAL) return this_ds @@ -333,14 +345,17 @@ 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, + throw_errors=True, ): """ Import CLM output """ + any_bad = False + # Import this_ds = import_ds(filename, my_vars=my_vars, my_vegtypes=my_vegtypes) @@ -410,7 +425,9 @@ 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): - check_v0_le_v1(this_ds_gs, var_list, both_nan_ok=True, throw_error=True) + any_bad = any_bad or check_v0_le_v1( + this_ds_gs, var_list, both_nan_ok=True, throw_error=throw_errors + ) # Check that prescribed calendars were obeyed if sdates_rx_ds: @@ -438,7 +455,7 @@ def import_output( raise RuntimeError("How to get NHARVEST_DISCREP for NHARVESTS > 2?") this_ds_gs["NHARVEST_DISCREP"] = (this_ds_gs["NHARVESTS"] == 2).astype(int) - return this_ds_gs + return this_ds_gs, any_bad def handle_zombie_crops(this_ds): 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 index d28bfda0e7..13668fc850 100644 --- a/python/ctsm/crop_calendars/generate_gdd20_baseline.py +++ b/python/ctsm/crop_calendars/generate_gdd20_baseline.py @@ -17,12 +17,13 @@ # 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 -VAR_LIST_IN = ["GDD0", "GDD8", "GDD10"] -VAR_LIST_IN = [x + "20" for x in VAR_LIST_IN] # TODO: Delete this once using the right variables -MISSING_FILL = -1 # Something negative to ensure that gddmaturity never changes (see PlantCrop) +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(): @@ -68,6 +69,28 @@ def _parse_args(): 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:]) @@ -76,7 +99,28 @@ def _parse_args(): if os.path.exists(args.output_file) and not args.overwrite: raise FileExistsError("Output file exists but --overwrite is not specified") - return args + # 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): @@ -95,14 +139,13 @@ def _get_cft_list(crop_list): "cotton", "irrigated_cotton"] """ - mgdcrop_list = utils.define_mgdcrop_list() cft_str_list = [] for crop_str in crop_list: - cft_str_list += [x for x in mgdcrop_list if crop_str in x] + 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): +def _get_gddn_for_cft(cft_str, variable): """ Given a CFT name, return the GDDN variable it uses. @@ -110,28 +153,28 @@ def _get_gddn_for_cft(cft_str): cft_str (str): E.g., "irrigated_temperate_corn" Returns: - str or None: Name of variable to use (e.g., "GDD8"). If crop isn't yet handled, return None. + 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 = "GDD0" + gddn = 0 gdd8_list_str = ["corn", "sugarcane", "miscanthus", "switchgrass"] if cft_str in _get_cft_list(gdd8_list_str): - gddn = "GDD8" + gddn = 8 gdd10_list_str = ["soybean"] if cft_str in _get_cft_list(gdd10_list_str): - gddn = "GDD10" + gddn = 10 - # TODO: Delete this once using the right variables if gddn is not None: - gddn += "20" + gddn_str = variable.replace("B", str(gddn)) - return gddn + return gddn, gddn_str def _get_output_varname(cft_str): @@ -160,19 +203,53 @@ def _add_time_axis(da_in): return da_out -def generate_gdd20_baseline(input_files, output_file, author): +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 """ - # Get input file list - input_files = input_files.split(sep=" ") + # 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) + 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") @@ -182,41 +259,43 @@ def generate_gdd20_baseline(input_files, output_file, author): # Set up a dummy DataArray to use for crops without an assigned GDDN variable dummy_da = xr.DataArray( - data=MISSING_FILL * np.ones_like(ds_in[VAR_LIST_IN[0]].values), - dims=ds_in[VAR_LIST_IN[0]].dims, - coords=ds_in[VAR_LIST_IN[0]].coords, + 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 - ds_out = xr.Dataset( - data_vars=None, - attrs={ - "author": author, - "created": dt.datetime.now().astimezone().isoformat(), - }, - ) - for cft_str in utils.define_mgdcrop_list(): + 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 = _get_gddn_for_cft(cft_str) + gddn, gddn_str = _get_gddn_for_cft(cft_str, variable) - # Fill any missing values with MISSING_FILL. This will mean that gddmaturity in these cells + # Fill any missing values with MISSING_RX_GDD_VAL. This will mean that gddmaturity there # never changes. - if gddn is None: - # Crop not handled yet? Fill it entirely with missing value + if gddn_str is None: + # Crop not handled yet? It's already filled with missing value this_da = dummy_da - long_name = "Dummy GDD20" print(" dummy GDD20") else: - this_da = ds_in[gddn].fillna(MISSING_FILL) + this_da = ds_in[gddn_str] # Already did ds_in.mean(dim="time") above this_da = _add_time_axis(this_da) - long_name = gddn - print(f" {gddn}") + print(f" {gddn_str}") + this_da = this_da.fillna(MISSING_RX_GDD_VAL) - # Add attributes + # 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" @@ -224,9 +303,14 @@ def generate_gdd20_baseline(input_files, output_file, author): 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) + ds_out.to_netcdf(output_file, format="NETCDF3_CLASSIC", encoding=encoding_dict) print("Done!") @@ -235,9 +319,12 @@ def main(): """ main() function for calling generate_gdd20_baseline.py from command line. """ - args = _parse_args() + 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 new file mode 100755 index 0000000000..123d40af38 --- /dev/null +++ b/python/ctsm/crop_calendars/interpolate_gdds.py @@ -0,0 +1,167 @@ +""" +Interpolate a maturity requirement (GDD) file +""" + +import os +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) +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir) +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm import ctsm_logging # pylint: disable=wrong-import-position +from ctsm.crop_calendars.cropcal_module import ( # pylint: disable=wrong-import-position + unexpected_negative_rx_gdd, +) + +logger = logging.getLogger(__name__) + +OUTPUT_FORMAT = "NETCDF3_CLASSIC" + + +def _file_missing(filepath, descriptor): + if not os.path.exists(filepath) and not os.path.exists(os.path.realpath(filepath)): + raise FileNotFoundError(f"{descriptor} not found: {filepath}") + return os.path.realpath(filepath) + + +def _setup_process_args(): + """Process input arguments + + Returns: + argparse.ArgumentParser: Arguments/options + """ + + # set up logging allowing user control + ctsm_logging.setup_logging_pre_config() + + parser = argparse.ArgumentParser( + description=("Interpolate a maturity requirement (GDD) file."), + ) + + # Define arguments + parser.add_argument( + "-i", + "--input-file", + help="Maturity requirement file to interpolate", + type=lambda x: _file_missing(x, "Input file"), + required=True, + ) + parser.add_argument( + "-t", + "--target-file", + help="File to whose coordinates the interpolation should take place", + type=lambda x: _file_missing(x, "Target file"), + required=True, + ) + parser.add_argument( + "-o", + "--output-file", + help="Where to save interpolated result", + 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.", + action="store_true", + required=False, + ) + parser.add_argument( + "--dry-run", + help="Check arguments but do not run.", + action="store_true", + required=False, + ) + ctsm_logging.add_logging_args(parser) + + # Get arguments + args = parser.parse_args(sys.argv[1:]) + ctsm_logging.process_logging_args(args) + + # Process arguments + if os.path.exists(os.path.realpath(args.output_file)) and not args.overwrite: + raise FileExistsError(f"Output file exists but --overwrite not given: {args.output_file}") + args.output_file = os.path.realpath(args.output_file) + + # Make directory for output file, if needed + if not args.dry_run: + parent_dir = os.path.dirname(args.output_file) + if not os.path.exists(parent_dir): + os.makedirs(parent_dir) + + return args + + +def interpolate_gdds(args): + """ + Do the interpolation and save + """ + + # Open inputs + ds_in = xr.open_dataset(args.input_file) + ds_target = xr.open_dataset(args.target_file) + + # Interpolate + ds_out = xr.Dataset() + for var in ds_in: + + # Check variable + 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 + + # Interpolate + da_out = ds_in[var].interp_like( + ds_target, + method="nearest", + kwargs={"fill_value": "extrapolate"}, # Otherwise you get NaNs at edges + ) + + if unexpected_negative_rx_gdd(da_out): + raise RuntimeError("Unexpected negative value") + + # Add to dataset + ds_out[var] = da_out + + # Finish up + ds_out.attrs["original"] = args.input_file + ds_out.attrs["interpolation_target"] = args.target_file + ds_out.attrs["interpolation_script"] = os.path.basename(__file__) + if not args.dry_run: + ds_out.to_netcdf(args.output_file, format=OUTPUT_FORMAT) + else: + print("Dry run looks good!") + + +def main(): + """ + Description + ----------- + Calls function that interpolates a maturity requirement (GDD) file. + """ + + args = _setup_process_args() + + interpolate_gdds(args) + + +if __name__ == "__main__": + main() diff --git a/python/ctsm/run_sys_tests.py b/python/ctsm/run_sys_tests.py index de93081504..b3a3b78379 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/CNPhenologyMod.F90 b/src/biogeochem/CNPhenologyMod.F90 index cfeee0b867..1fb13d32eb 100644 --- a/src/biogeochem/CNPhenologyMod.F90 +++ b/src/biogeochem/CNPhenologyMod.F90 @@ -1811,7 +1811,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_streams + use clm_varctl , only : use_cropcal_rx_swindows ! ! !ARGUMENTS: integer , intent(in) :: num_pcropp ! number of prog crop patches in filter @@ -1976,7 +1976,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 @@ -2245,7 +2245,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 @@ -2562,7 +2562,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, adapt_cropcal_rx_cultivar_gdds + 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 @@ -2699,7 +2699,7 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & 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 SSR: Set maximum and minimum gddmaturity + !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)) @@ -2726,12 +2726,17 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, & 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 f35c2fffe6..0f650a4a9f 100644 --- a/src/biogeochem/CropType.F90 +++ b/src/biogeochem/CropType.F90 @@ -53,6 +53,8 @@ module CropType 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] @@ -237,6 +239,8 @@ subroutine InitAllocate(this, bounds) 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 @@ -360,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) diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 0144f1ef18..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, 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,78 +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 - - ! 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 GDD0 with uninitialized values. - 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 - call markreset_accum_field('GDD0', p) ! 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 - - ! 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 GDD8 with uninitialized values. - 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 - call markreset_accum_field('GDD8', p) ! 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 - - ! 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 GDD10 with uninitialized values. - 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 - call markreset_accum_field('GDD10', p) ! 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 413ddbefce..85218a0e1c 100644 --- a/src/cpl/share_esmf/cropcalStreamMod.F90 +++ b/src/cpl/share_esmf/cropcalStreamMod.F90 @@ -38,9 +38,12 @@ module cropcalStreamMod type(shr_strdata_type) :: sdat_cropcal_swindow_end ! sowing window end 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 @@ -49,6 +52,10 @@ module cropcalStreamMod 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__ @@ -105,7 +112,11 @@ subroutine cropcal_init(bounds) stream_fldFileName_gdd20_baseline, & stream_meshfile_cropcal, & cropcals_rx, & - cropcals_rx_adapt + 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_swindows = 1 ! first year in sowing window streams to use @@ -120,16 +131,22 @@ subroutine cropcal_init(bounds) 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 @@ -158,6 +175,10 @@ subroutine cropcal_init(bounds) 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,*) @@ -174,9 +195,14 @@ subroutine cropcal_init(bounds) 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,*) @@ -186,9 +212,11 @@ subroutine cropcal_init(bounds) use_cropcal_rx_swindows = stream_fldFileName_swindow_start /= '' use_cropcal_rx_cultivar_gdds = stream_fldFileName_cultivar_gdds /= '' adapt_cropcal_rx_cultivar_gdds = stream_fldFileName_gdd20_baseline /= '' - use_cropcal_streams = use_cropcal_rx_swindows .or. use_cropcal_rx_cultivar_gdds + 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, & @@ -247,6 +275,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, & @@ -278,6 +307,7 @@ subroutine cropcal_init(bounds) ! 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, & @@ -286,7 +316,7 @@ subroutine cropcal_init(bounds) model_mesh = mesh, & stream_meshfile = trim(stream_meshfile_cropcal), & stream_lev_dimname = 'null', & - stream_mapalgo = 'bilinear', & + stream_mapalgo = 'nn', & stream_filenames = (/trim(stream_fldFileName_gdd20_baseline)/), & stream_fldlistFile = stream_varnames_gdd20_baseline, & stream_fldListModel = stream_varnames_gdd20_baseline, & @@ -304,6 +334,77 @@ subroutine cropcal_init(bounds) 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 !================================================================ @@ -347,7 +448,26 @@ subroutine cropcal_advance( bounds ) end if end if - ! GDD20 baseline values do not have an associated time axis and thus will not be advanced here + ! 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) ) @@ -393,15 +513,21 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, crop_inst) 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__) @@ -459,8 +585,8 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, 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) @@ -469,23 +595,23 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, 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 @@ -494,7 +620,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, 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 @@ -562,7 +688,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, crop_inst) deallocate(dataptr2d_cultivar_gdds) allocate(dataptr2d_gdd20_baseline(lsize, ncft)) - if (adapt_cropcal_rx_cultivar_gdds .and. init) then + if (adapt_cropcal_rx_cultivar_gdds) then ! Read GDD20 baselines from input files ! Starting with npcropmin will skip generic crops do n = 1, ncft @@ -613,6 +739,91 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, init, crop_inst) deallocate(dataptr2d_gdd20_baseline) + ! Read prescribed gdd20 season start dates from input files + allocate(dataptr2d_gdd20_season_start(lsize, ncft)) + dataptr2d_gdd20_season_start(:,:) = -1._r8 + allocate(dataptr2d_gdd20_season_end (lsize, ncft)) + dataptr2d_gdd20_season_end(:,:) = -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 = 1,lsize + + ! 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/clm_driver.F90 b/src/main/clm_driver.F90 index 2a3cef4b8b..e660ab9d8d 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -1374,7 +1374,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_varctl.F90 b/src/main/clm_varctl.F90 index fe67380fd2..678f386f23 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -392,6 +392,7 @@ module clm_varctl 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 46d9e9958a..634128798b 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -302,7 +302,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/ & @@ -711,6 +711,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) @@ -983,6 +984,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/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/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 #