diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index e14ce11599..b093c298e9 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -11,11 +11,15 @@ b88e1cd1b28e3609684c79a2ec0e88f26cfc362b b771971e3299c4fa56534b93421f7a2b9c7282fd 9de88bb57ea9855da408cbec1dc8acb9079eda47 8bc4688e52ea23ef688e283698f70a44388373eb +4ee49e3e516ca7dee5df378f65664f93a7db4415 +0207bc98dd5c75cd69a0e788bc53e41093712f5c +e4d38681df23ccca0ae29581a45f8362574e0630 0a5a9e803b56ec1bbd6232eff1c99dbbeef25eb7 810cb346f05ac1aabfff931ab1a2b7b584add241 5933b0018f8e29413e30dda9b906370d147bad45 # Ran SystemTests and python/ctsm through black python formatter 5364ad66eaceb55dde2d3d598fe4ce37ac83a93c 8056ae649c1b37f5e10aaaac79005d6e3a8b2380 +0bc3f00115d86d026a977918661c93779b3b19f9 540b256d1f3382f4619d7b0877c32d54ce5c40b6 8a168bb0895f4f2421608dd2589398e13a6663e6 diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 65cef8b77e..5769af0ea2 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -613,7 +613,7 @@ sub process_namelist_user_input { process_namelist_commandline_infile($opts, $definition, $nl, $envxml_ref); # Apply the commandline options and make sure the user didn't change it above - process_namelist_commandline_options($opts, $nl_flags, $definition, $defaults, $nl, $physv); + process_namelist_commandline_options($opts, $nl_flags, $definition, $defaults, $nl, $envxml_ref, $physv); # The last two process command line arguments for usr_name and use_case # They require that process_namelist_commandline_options was called before this @@ -634,10 +634,10 @@ sub process_namelist_commandline_options { # Obtain default values for the following build-namelist input arguments # : res, mask, ssp_rcp, sim_year, sim_year_range, and clm_accelerated_spinup. - my ($opts, $nl_flags, $definition, $defaults, $nl, $physv) = @_; + my ($opts, $nl_flags, $definition, $defaults, $nl, $envxml_ref, $physv) = @_; setup_cmdl_chk_res($opts, $defaults); - setup_cmdl_resolution($opts, $nl_flags, $definition, $defaults); + setup_cmdl_resolution($opts, $nl_flags, $definition, $defaults, $envxml_ref); setup_cmdl_mask($opts, $nl_flags, $definition, $defaults, $nl); setup_cmdl_configuration_and_structure($opts, $nl_flags, $definition, $defaults, $nl); setup_cmdl_bgc($opts, $nl_flags, $definition, $defaults, $nl); @@ -668,7 +668,7 @@ sub setup_cmdl_chk_res { } sub setup_cmdl_resolution { - my ($opts, $nl_flags, $definition, $defaults) = @_; + my ($opts, $nl_flags, $definition, $defaults, $envxml_ref) = @_; my $var = "res"; my $val; @@ -686,16 +686,30 @@ sub setup_cmdl_resolution { $val = "e_string( $nl_flags->{'res'} ); if ( ! $definition->is_valid_value( $var, $val ) ) { my @valid_values = $definition->get_valid_values( $var ); - if ( ! defined($opts->{'clm_usr_name'}) || $nl_flags->{'res'} ne $opts->{'clm_usr_name'} ) { + if ( $nl_flags->{'res'} ne "CLM_USRDAT" ) { $log->fatal_error("$var has a value ($val) that is NOT valid. Valid values are: @valid_values"); } } } + if ( $nl_flags->{'res'} eq "CLM_USRDAT" ) { + if ( ! defined($opts->{'clm_usr_name'}) ) { + $log->fatal_error("Resolution is CLM_USRDAT, but --clm_usr_name option is NOT set, and it is required for CLM_USRDAT resolutions"); + } + } + # # For NEON sites - if ($nl_flags->{'res'} =~ /NEON/) { - $nl_flags->{'neon'} = ".true." - } else { - $nl_flags->{'neon'} = ".false." + # + $nl_flags->{'neon'} = ".false."; + $nl_flags->{'neonsite'} = ""; + if ( $nl_flags->{'res'} eq "CLM_USRDAT" ) { + if ( $opts->{'clm_usr_name'} eq "NEON" ) { + $nl_flags->{'neon'} = ".true."; + $nl_flags->{'neonsite'} = $envxml_ref->{'NEONSITE'}; + $log->verbose_message( "This is a NEON site with NEONSITE = " . $nl_flags->{'neonsite'} ); + } + } + if ( ! &value_is_true( $nl_flags->{'neon'} ) ) { + $log->verbose_message( "This is NOT a NEON site" ); } } @@ -1572,7 +1586,7 @@ sub process_namelist_inline_logic { setup_logic_grainproduct($opts, $nl_flags, $definition, $defaults, $nl, $physv); setup_logic_soilstate($opts, $nl_flags, $definition, $defaults, $nl); setup_logic_demand($opts, $nl_flags, $definition, $defaults, $nl); - setup_logic_surface_dataset($opts, $nl_flags, $definition, $defaults, $nl); + setup_logic_surface_dataset($opts, $nl_flags, $definition, $defaults, $nl, $envxml_ref); setup_logic_dynamic_subgrid($opts, $nl_flags, $definition, $defaults, $nl); if ( remove_leading_and_trailing_quotes($nl_flags->{'clm_start_type'}) ne "branch" ) { setup_logic_initial_conditions($opts, $nl_flags, $definition, $defaults, $nl, $physv); @@ -2031,6 +2045,13 @@ sub setup_logic_snicar_methods { sub setup_logic_snow { my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_thermal_cond_method' ); + + my $var = $nl->get_value('snow_thermal_cond_method'); + if ( $var ne "'Jordan1991'" && $var ne "'Sturm1997'" ) { + $log->fatal_error("$var is incorrect entry for the namelist variable snow_thermal_cond_method; expected Jordan1991 or Sturm1997"); + } + my $numrad_snw = $nl->get_value('snicar_numrad_snw'); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'fsnowoptics', 'snicar_numrad_snw' => $numrad_snw); @@ -2290,6 +2311,7 @@ sub setup_logic_demand { $settings{'use_lch4'} = $nl_flags->{'use_lch4'}; $settings{'use_nitrif_denitrif'} = $nl_flags->{'use_nitrif_denitrif'}; $settings{'use_crop'} = $nl_flags->{'use_crop'}; + $settings{'neon'} = $nl_flags->{'neon'}; my $demand = $nl->get_value('clm_demand'); if (defined($demand)) { @@ -2342,7 +2364,7 @@ sub setup_logic_surface_dataset { # consistent with it # MUST BE AFTER: setup_logic_demand which is where flanduse_timeseries is set # - my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + my ($opts, $nl_flags, $definition, $defaults, $nl, $xmlvar_ref) = @_; $nl_flags->{'flanduse_timeseries'} = "null"; my $flanduse_timeseries = $nl->get_value('flanduse_timeseries'); @@ -2367,26 +2389,42 @@ sub setup_logic_surface_dataset { if ( ! &value_is_true($nl_flags->{'use_fates'}) ) { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, 'hgrid'=>$nl_flags->{'res'}, 'ssp_rcp'=>$nl_flags->{'ssp_rcp'}, + 'neon'=>$nl_flags->{'neon'}, 'neonsite'=>$nl_flags->{'neonsite'}, 'sim_year'=>$nl_flags->{'sim_year'}, 'irrigate'=>".true.", 'use_vichydro'=>$nl_flags->{'use_vichydro'}, - 'use_crop'=>".true.", 'glc_nec'=>$nl_flags->{'glc_nec'}, 'nofail'=>1); + 'use_crop'=>".true.", 'glc_nec'=>$nl_flags->{'glc_nec'}, 'use_fates'=>$nl_flags->{'use_fates'}, 'nofail'=>1); } # If didn't find the crop version check for the exact match - if ( ! defined($nl->get_value($var) ) ) { + my $fsurdat = $nl->get_value($var); + if ( ! defined($fsurdat) ) { if ( ! &value_is_true($nl_flags->{'use_fates'}) ) { $log->verbose_message( "Crop version of $var NOT found, searching for an exact match" ); } add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, 'hgrid'=>$nl_flags->{'res'}, 'ssp_rcp'=>$nl_flags->{'ssp_rcp'}, 'use_vichydro'=>$nl_flags->{'use_vichydro'}, - 'sim_year'=>$nl_flags->{'sim_year'}, 'irrigate'=>$nl_flags->{'irrigate'}, + 'sim_year'=>$nl_flags->{'sim_year'}, 'irrigate'=>$nl_flags->{'irrigate'}, 'use_fates'=>$nl_flags->{'use_fates'}, + 'neon'=>$nl_flags->{'neon'}, 'neonsite'=>$nl_flags->{'neonsite'}, 'use_crop'=>$nl_flags->{'use_crop'}, 'glc_nec'=>$nl_flags->{'glc_nec'}, 'nofail'=>1 ); - if ( ! defined($nl->get_value($var) ) ) { + if ( ! defined($fsurdat) ) { $log->verbose_message( "Exact match of $var NOT found, searching for version with irrigate true" ); } add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, 'hgrid'=>$nl_flags->{'res'}, 'ssp_rcp'=>$nl_flags->{'ssp_rcp'}, 'use_vichydro'=>$nl_flags->{'use_vichydro'}, - 'sim_year'=>$nl_flags->{'sim_year'}, 'irrigate'=>".true.", + 'sim_year'=>$nl_flags->{'sim_year'}, 'irrigate'=>".true.", 'use_fates'=>$nl_flags->{'use_fates'}, + 'neon'=>$nl_flags->{'neon'}, 'neonsite'=>$nl_flags->{'neonsite'}, 'use_crop'=>$nl_flags->{'use_crop'}, 'glc_nec'=>$nl_flags->{'glc_nec'} ); } + # + # Expand the XML variables for NEON cases so that NEONSITE will be used + # + if ( &value_is_true($nl_flags->{'neon'}) ) { + my $fsurdat = $nl->get_value($var); + my $newval = SetupTools::expand_xml_var( $fsurdat, $xmlvar_ref ); + if ( $newval ne $fsurdat ) { + my $group = $definition->get_group_name($var); + $nl->set_variable_value($group, $var, $newval); + $log->verbose_message( "This is a NEON site and the fsurdat file selected is: $newval" ); + } + } } #------------------------------------------------------------------------------- diff --git a/bld/env_run.xml b/bld/env_run.xml index 8bf59d0911..f3b7467168 100644 --- a/bld/env_run.xml +++ b/bld/env_run.xml @@ -9,5 +9,6 @@ Sample env_run.xml file that allows build-namelist to be run for testing in this --> + diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 2e57391df7..ec3bc2f2b4 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -445,6 +445,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case). 1.e9 SwensonLawrence2012 +Jordan1991 + +lnd/clm2/surfdata_map/NEON/16PFT_mixed/surfdata_1x1_NEON_${NEONSITE}_hist_16pfts_Irrig_CMIP6_simyr2000_c230120.nc + +lnd/clm2/surfdata_map/NEON/surfdata_1x1_NEON_${NEONSITE}_hist_78pfts_CMIP6_simyr2000_c230601.nc + + + + + lnd/clm2/surfdata_map/landuse.timeseries_0.9x1.25_hist_16pfts_Irrig_CMIP6_simyr1850-2015_c170824.nc -flanduse_timeseries -flanduse_timeseries +null +null +flanduse_timeseries +flanduse_timeseries diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index bca1beca6a..93306b32d5 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2845,6 +2845,11 @@ NiuYang2007: Niu and Yang 2007 SwensonLawrence2012: Swenson and Lawrence 2012 + +Parameterization to use for snow thermal conductivity + + diff --git a/bld/namelist_files/use_cases/2018-PD_transient.xml b/bld/namelist_files/use_cases/2018-PD_transient.xml index d838efbd00..96f14207ad 100644 --- a/bld/namelist_files/use_cases/2018-PD_transient.xml +++ b/bld/namelist_files/use_cases/2018-PD_transient.xml @@ -1,8 +1,12 @@ + + -Simulate transient land-use, and aerosol deposition changes from 2018 to current day with a mix of historical data, and future scenario data +Simulate transient Nitrogen-deposition, aerosol deposition, urban, and fire related (pop-density, lightning) changes from 2018 to current day with a mix of historical data, and future scenario data +Simulate transient Nitrogen-deposition, aerosol deposition, urban, and fire related (pop-density, lightning) changes from 2018 to current day with a mix of historical data, and future scenario data +Simulate transient urban and aerosol deposition changes from 2018 to current day with a mix of historical data, and future scenario data diff --git a/bld/namelist_files/use_cases/2018_control.xml b/bld/namelist_files/use_cases/2018_control.xml index e5e572d749..28554074c4 100644 --- a/bld/namelist_files/use_cases/2018_control.xml +++ b/bld/namelist_files/use_cases/2018_control.xml @@ -1,5 +1,8 @@ + + + Conditions to simulate 2018 land-use diff --git a/bld/unit_testers/build-namelist_test.pl b/bld/unit_testers/build-namelist_test.pl index da4201d68f..4c717a58bf 100755 --- a/bld/unit_testers/build-namelist_test.pl +++ b/bld/unit_testers/build-namelist_test.pl @@ -42,7 +42,7 @@ sub make_env_run { my %settings = @_; # Set default settings - my %env_vars = ( DIN_LOC_ROOT=>"MYDINLOCROOT", GLC_TWO_WAY_COUPLING=>"FALSE" ); + my %env_vars = ( DIN_LOC_ROOT=>"MYDINLOCROOT", GLC_TWO_WAY_COUPLING=>"FALSE", NEONSITE=>"" ); # Set any settings that came in from function call foreach my $item ( keys(%settings) ) { $env_vars{$item} = $settings{$item}; @@ -381,7 +381,7 @@ sub cat_and_create_namelistinfile { "JORN", "LAJA", "MOAB", "OAES", "OSBS", "SCBI", "SOAP", "STER", "TOOL", "UNDE", "YELL" ) { - &make_env_run(); + &make_env_run( NEONSITE=>"$site" ); # # Concatonate default usermods and specific sitetogether expanding env variables while doing that # @@ -400,7 +400,7 @@ sub cat_and_create_namelistinfile { # # Now run the site # - my $options = "-res CLM_USRDAT -clm_usr_name NEON -no-megan -bgc bgc -sim_year 2018 -infile $namelistfile"; + my $options = "--res CLM_USRDAT --clm_usr_name NEON --no-megan --bgc bgc --use_case 2018_control --infile $namelistfile"; eval{ system( "$bldnml -envxml_dir . $options > $tempfile 2>&1 " ); }; is( $@, '', "options: $options" ); $cfiles->checkfilesexist( "$options", $mode ); diff --git a/cime_config/buildnml b/cime_config/buildnml index e239f0ec58..84e1581406 100755 --- a/cime_config/buildnml +++ b/cime_config/buildnml @@ -136,7 +136,6 @@ def buildnml(case, caseroot, compname): clmusr = "" if lnd_grid == "CLM_USRDAT": clm_usrdat_name = case.get_value("CLM_USRDAT_NAME") - lnd_grid = clm_usrdat_name clmusr = " -clm_usr_name %s " % clm_usrdat_name # Write warning about initial condition data if "NEON" in clm_usrdat_name and clm_force_coldstart == "off": diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index b696fadcf2..8e72abb8eb 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -57,13 +57,6 @@ - - - FAIL - #2236 - - - diff --git a/cime_config/usermods_dirs/NEON/defaults/shell_commands b/cime_config/usermods_dirs/NEON/defaults/shell_commands index 437201f2b2..1f5427ca98 100644 --- a/cime_config/usermods_dirs/NEON/defaults/shell_commands +++ b/cime_config/usermods_dirs/NEON/defaults/shell_commands @@ -32,7 +32,7 @@ else fi # If needed for SP simulations: & set history file variables -if [[ $compset =~ .*CLM[0-9]+%.*SP.* ]]; then +if [[ $compset =~ .*CLM[0-9]+%[^_]*SP.* ]]; then if [[ $TEST != "TRUE" ]]; then ./xmlchange STOP_OPTION=nyears fi diff --git a/doc/ChangeLog b/doc/ChangeLog index 66fc04a752..9405e19922 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,168 @@ =============================================================== +Tag name: ctsm5.1.dev152 +Originator(s): multiple (tking (Teagan King); slevis (Sam Levis); AdrienDams (Adrien Damseaux); afoster (Adrianna Foster); samrabin (Sam Rabin); ekluzek (Erik Kluzek); wwieder (Will Wieder); sacks (Bill Sacks); a few others listed below) +Date: Tue Nov 14 17:09:43 MST 2023 +One-line Summary: Mv tools to /python and add tests; add snow_thermal_cond_method; a few fixes / refactors + +Purpose and description of changes +---------------------------------- + +#2156 tking, slevis +Move the following scripts to /python/ctsm/site_and_regional +and make wrapper scripts for them in /tools/site_and_regional: +- run_neon.py +- neon_surf_wrapper.py +- modify_singlept_site_neon.py + +Add unit testing for: +- iso_utils +- modify_singlept_site_neon +- neon_surf_wrapper +- run_neon + +Add system testing for: +- modify_singlept_site_neon +- run_neon + +#2148 Adrien Damseaux (AWI, Germany), Victoria Dutch, Leanne Wake +Add namelist option snow_thermal_cond_method to select between Jordan (1991) (default) and +Sturm et al. (1997). Sturm option described for single point runs by Dutch et al. (2022). + +#2233 afoster, sacks +Fix a compiler error (for GNU 13.2) within cropcalStreamMod. +Simple fix was to change whole-array assignments/references for the starts and ends arrays to specifically +reference bounds (begp and endp). + +#2235 srabin, wwieder +Refactor ssp_anomaly_forcing script to make it easier to read and more amenable to future development. +- Adds --output-dir option; default ./anomaly_forcing reproduces previous behavior +- Makes synonyms for options with hyphens replacing underscores + +#2237 srabin +Add the following fields to restart files: +- repr_grainc_to_seed_perharv_patch +- swindow_starts_thisyr_patch +- swindow_ends_thisyr_patch + +#2044 ekluzek +More confined regular expression for NEON and a few simple fixes. + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[ ] clm5_1 + +[ ] clm5_0 + +[ ] ctsm5_0-nwp + +[ ] clm4_5 + + +Bugs fixed or introduced +------------------------ +CTSM issues fixed (include CTSM Issue #): +Closes #2156 Fixes #1441 +Closes #2148 +Closes #2233 Fixes #2232 +Closes #2235 +Closes #2237 Fixes #2236 +Closes #2044 Fixes #2039 Fixes #2103 Fixes #2028 Fixes #1506 Fixes #1499 + +Known Issues: +pylint errors from previous work remain in this tag. + +Notes of particular relevance for users +--------------------------------------- +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): +#2156 New wrapper scripts don't have .py suffixes. +#2148 New namelist option snow_thermal_cond_method as described above. +#2133 None +#2135 New --output-dir option; default ./anomaly_forcing reproduces previous behavior. +Also makes synonyms for options with hyphens replacing underscores. +#2137 None +#2044 None + +Notes of particular relevance for developers: +--------------------------------------------- +Changes to tests or testing: +#2156 Numerous changes were made to include new tests. +README.md for testing was updated to clarify that arguments should be used. + +Testing summary: +---------------- +[Remove any lines that don't apply.] + + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + cheyenne - + + tools-tests (test/tools) (if tools have been changed): + + cheyenne - + + python testing (if python code has changed; see instructions in python/README.md; document testing done): + + (any machine) - cheyenne OK (pylint suggestions from previous work remain) + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + cheyenne ---- OK + izumi ------- OK, the following PASS/FAILs are expected: + +PASS ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-cropMonthlyNoinitial COMPARE_base_rest (UNEXPECTED: expected FAIL) +FAIL ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-cropMonthlyNoinitial BASELINE ctsm5.1.dev151: DIFF + +FAIL SMS_Ld10_D_Mmpi-serial.CLM_USRDAT.I1PtClm51Bgc.izumi_nag.clm-default--clm-NEON-NIWO BASELINE ctsm5.1.dev151: DIFF +FAIL SMS_Ld10_D_Mmpi-serial.CLM_USRDAT.I1PtClm51Bgc.izumi_nag.clm-NEON-MOAB--clm-PRISM BASELINE ctsm5.1.dev151: DIFF + + fates tests: (give name of baseline if different from CTSM tagname, normally fates baselines are fates--) + cheyenne ---- + izumi ------- + + any other testing (give details below): + +If the tag used for baseline comparisons was NOT the previous tag, note that here: + +Answer changes +-------------- +Changes answers relative to baseline: +#2156 NO +#2148 NO +#2233 NO +#2235 NO, adds attributes to write_climo files' dimension variables +#2237 ONLY Smallville "no initial" restarts; specifically, this previously +failing (COMPARE_base_rest) aux_clm test +ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-cropMonthlyNoinitial +now differs from the baseline as follows: + SUMMARY of cprnc: + A total number of 76 fields were compared + and 3 had differences in fill patterns + A total number of 2 fields could not be analyzed + diff_test: the two files seem to be DIFFERENT +#2044 ONLY the NEON tests listed above due to the one-line change in +cime_config/usermods_dirs/NEON/defaults/shell_commands in #2044 + +Other details +------------- +Pull Requests that document the changes (include PR ids): + https://github.com/ESCOMP/ctsm/pull/2156 + https://github.com/ESCOMP/ctsm/pull/2148 + https://github.com/ESCOMP/ctsm/pull/2233 + https://github.com/ESCOMP/ctsm/pull/2235 + https://github.com/ESCOMP/ctsm/pull/2237 + https://github.com/ESCOMP/ctsm/pull/2044 + +=============================================================== +=============================================================== Tag name: ctsm5.1.dev151 Originator(s): rgknox (Ryan Knox,LAWRENCE BERKELEY NATIONAL LABORATORY,510-495-2153) Date: Sat Nov 11 16:53:01 MST 2023 diff --git a/doc/ChangeSum b/doc/ChangeSum index 9ef20790d9..6af88be0e4 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.1.dev152 multiple 11/14/2023 Mv tools to /python and add tests; add snow_thermal_cond_method; a few fixes / refactors ctsm5.1.dev151 rgknox 11/11/2023 Fixes to FATES long run restarts ctsm5.1.dev150 rgknox 11/06/2023 FATES API fix to support future fates npp-fixation coupling, and urgent coupling fixes with E3SM. ctsm5.1.dev149 samrabin 11/03/2023 Rearrange leaf/stem "harvest" and fix soil gas diffusivity diff --git a/python/README.md b/python/README.md index cf3b893084..c40f55c6c7 100644 --- a/python/README.md +++ b/python/README.md @@ -47,7 +47,8 @@ thing, but support different options: 2. via `./run_ctsm_py_tests` You can specify various arguments to this; run `./run_ctsm_py_tests - -h` for details + -h` for details. Please specify either --unit or --sys rather than + not including any arguments. In any configuration where you run the system tests, you need to first execute `module load nco`. diff --git a/python/ctsm/run_ctsm_py_tests.py b/python/ctsm/run_ctsm_py_tests.py index 0542dc41cb..8b39d69afa 100644 --- a/python/ctsm/run_ctsm_py_tests.py +++ b/python/ctsm/run_ctsm_py_tests.py @@ -45,7 +45,10 @@ def main(description): def _commandline_args(description): - """Parse and return command-line arguments""" + """Parse and return command-line arguments + Note that run_ctsm_py_tests is not intended to be + used without argument specifications + """ parser = argparse.ArgumentParser( description=description, formatter_class=argparse.RawTextHelpFormatter ) diff --git a/tools/site_and_regional/modify_singlept_site_neon.py b/python/ctsm/site_and_regional/modify_singlept_site_neon.py similarity index 86% rename from tools/site_and_regional/modify_singlept_site_neon.py rename to python/ctsm/site_and_regional/modify_singlept_site_neon.py index e135760a48..1013ba944c 100755 --- a/tools/site_and_regional/modify_singlept_site_neon.py +++ b/python/ctsm/site_and_regional/modify_singlept_site_neon.py @@ -31,81 +31,34 @@ """ # TODO (NS) # --[] If subset file not found run subset_data.py -# --[] List of valid neon sites for all scripts come from one place. # --[] Download files only when available. # Import libraries from __future__ import print_function +import argparse +from datetime import date +from getpass import getuser +import glob +import logging import os import sys -import glob -import argparse import requests -import logging import numpy as np import pandas as pd import xarray as xr from packaging import version -from datetime import date -from getpass import getuser - +from ctsm.path_utils import path_to_ctsm_root myname = getuser() # -- valid neon sites -valid_neon_sites = [ - "ABBY", - "BARR", - "BART", - "BLAN", - "BONA", - "CLBJ", - "CPER", - "DCFS", - "DEJU", - "DELA", - "DSNY", - "GRSM", - "GUAN", - "HARV", - "HEAL", - "JERC", - "JORN", - "KONA", - "KONZ", - "LAJA", - "LENO", - "MLBS", - "MOAB", - "NIWO", - "NOGP", - "OAES", - "ONAQ", - "ORNL", - "OSBS", - "PUUM", - "RMNP", - "SCBI", - "SERC", - "SJER", - "SOAP", - "SRER", - "STEI", - "STER", - "TALL", - "TEAK", - "TOOL", - "TREE", - "UKFS", - "UNDE", - "WOOD", - "WREF", - "YELL", -] +valid_neon_sites = glob.glob( + os.path.join(path_to_ctsm_root(), "cime_config", "usermods_dirs", "NEON", "[!d]*") +) def get_parser(): @@ -153,14 +106,14 @@ def get_parser(): parser.add_argument( "--inputdata-dir", help=""" - Directory to write updated single point surface dataset. + Directory containing standard input files from CESM input data such as the surf_soildepth_file. [default: %(default)s] """, action="store", dest="inputdatadir", type=str, required=False, - default="/glade/p/cesmdata/cseg/inputdata" + default="/glade/p/cesmdata/cseg/inputdata", ) parser.add_argument( "-d", @@ -233,10 +186,7 @@ def get_neon(neon_dir, site_name): print("Download finished successfully for", site_name) elif response.status_code == 404: sys.exit( - "Data for this site " - + site_name - + " was not available on the neon server:" - + url + "Data for this site " + site_name + " was not available on the neon server:" + url ) print("Download exit status code: ", response.status_code) @@ -270,12 +220,12 @@ def find_surffile(surf_dir, site_name, pft_16): """ if pft_16: - sf_name = "surfdata_1x1_NEON_"+site_name+"*hist_16pfts_Irrig_CMIP6_simyr2000_*.nc" + sf_name = "surfdata_1x1_NEON_" + site_name + "*hist_16pfts_Irrig_CMIP6_simyr2000_*.nc" else: - sf_name = "surfdata_1x1_NEON_" +site_name+"*hist_78pfts_CMIP6_simyr2000_*.nc" + sf_name = "surfdata_1x1_NEON_" + site_name + "*hist_78pfts_CMIP6_simyr2000_*.nc" - print (os.path.join(surf_dir , sf_name)) - surf_file = sorted(glob.glob(os.path.join(surf_dir , sf_name))) + print(os.path.join(surf_dir, sf_name)) + surf_file = sorted(glob.glob(os.path.join(surf_dir, sf_name))) if len(surf_file) > 1: print("The following files found :", *surf_file, sep="\n- ") @@ -287,10 +237,14 @@ def find_surffile(surf_dir, site_name, pft_16): surf_file = surf_file[0] else: sys.exit( - "Surface data for this site " + str(site_name) + " was not found:" + str(surf_dir) + str(sf_name) + - "." + - "\n" + - "Please run ./subset_data.py for this site." + "Surface data for this site " + + str(site_name) + + " was not found:" + + str(surf_dir) + + str(sf_name) + + "." + + "\n" + + "Please run ./subset_data.py for this site." ) return surf_file @@ -298,7 +252,7 @@ def find_surffile(surf_dir, site_name, pft_16): def find_soil_structure(args, surf_file): """ Function for finding surface dataset soil - strucutre using surface data metadata. + structure using surface data metadata. In CLM surface data, soil layer information is in a file from surface data metadata @@ -324,10 +278,8 @@ def find_soil_structure(args, surf_file): print("------------") # print (f1.attrs["Soil_texture_raw_data_file_name"]) - clm_input_dir = os.path.join( args.inputdatadir, "lnd/clm2/rawdata/" ) - surf_soildepth_file = os.path.join( - clm_input_dir, f1.attrs["Soil_texture_raw_data_file_name"] - ) + clm_input_dir = os.path.join(args.inputdatadir, "lnd/clm2/rawdata/") + surf_soildepth_file = os.path.join(clm_input_dir, f1.attrs["Soil_texture_raw_data_file_name"]) if os.path.exists(surf_soildepth_file): print( @@ -345,9 +297,7 @@ def find_soil_structure(args, surf_file): else: sys.exit( - "Cannot find soil structure file : " - + surf_soildepth_file - + "for the surface dataset." + "Cannot find soil structure file : " + surf_soildepth_file + "for the surface dataset." ) return soil_bot, soil_top @@ -356,7 +306,7 @@ def find_soil_structure(args, surf_file): def update_metadata(nc, surf_file, neon_file, zb_flag): """ Function for updating modified surface dataset - metadat for neon sites. + metadata for neon sites. Args: nc (xr Dataset): netcdf file including updated neon surface data @@ -492,9 +442,9 @@ def download_file(url, fname): elif response.status_code == 404: print("File " + fname + "was not available on the neon server:" + url) except Exception as err: - print ('The server could not fulfill the request.') - print ('Something went wrong in downloading', fname) - print ('Error code:', err.code) + print("The server could not fulfill the request.") + print("Something went wrong in downloading", fname) + print("Error code:", err.code) def fill_interpolate(f2, var, method): @@ -526,7 +476,6 @@ def fill_interpolate(f2, var, method): def main(): - args = get_parser().parse_args() # -- debugging option @@ -536,10 +485,13 @@ def main(): # Check if pandas is a recent enough version pdvers = pd.__version__ if version.parse(pdvers) < version.parse("1.1.0"): - sys.exit("The pandas version in your python environment is too old, update to a newer version of pandas (>=1.1.0): version=%s", pdvers ) - + sys.exit( + """The pandas version in your python environment is too old, + update to a newer version of pandas (>=1.1.0): version=%s""", + pdvers, + ) - file_time = check_neon_time() + # file_time = check_neon_time() # -- specify site from which to extract data site_name = args.site_name @@ -550,12 +502,9 @@ def main(): # -- directory structure current_dir = os.getcwd() - parent_dir = os.path.dirname(current_dir) clone_dir = os.path.abspath(os.path.join(__file__, "../../../..")) neon_dir = os.path.join(clone_dir, "neon_surffiles") - print("Present Directory", current_dir) - # -- download neon data if needed neon_file = get_neon(neon_dir, site_name) @@ -575,18 +524,11 @@ def main(): # better suggestion by WW to write dzsoi to neon surface dataset # This todo needs to go to the subset_data - # TODO Will: if I sum them up , are they 3.5? (m) YES - print("soil_top:", soil_top) - print("soil_bot:", soil_bot) - print("Sum of soil top depths :", sum(soil_top)) - print("Sum of soil bottom depths :", sum(soil_bot)) - soil_top = np.cumsum(soil_top) soil_bot = np.cumsum(soil_bot) soil_mid = 0.5 * (soil_bot - soil_top) + soil_top # print ("Cumulative sum of soil bottom depths :", sum(soil_bot)) - obs_top = df["biogeoTopDepth"] / 100 obs_bot = df["biogeoBottomDepth"] / 100 # -- Mapping surface dataset and neon soil levels @@ -635,12 +577,10 @@ def main(): # -- Check to make sure the rounded oc is not higher than carbon_tot. # -- Use carbon_tot if estimated_oc is bigger than carbon_tot. - if estimated_oc > carbon_tot: - estimated_oc = carbon_tot + estimated_oc = min(estimated_oc, carbon_tot) layer_depth = ( - df["biogeoBottomDepth"][bin_index[soil_lev]] - - df["biogeoTopDepth"][bin_index[soil_lev]] + df["biogeoBottomDepth"][bin_index[soil_lev]] - df["biogeoTopDepth"][bin_index[soil_lev]] ) # f2["ORGANIC"][soil_lev] = estimated_oc * bulk_den / 0.58 @@ -709,9 +649,7 @@ def main(): print("Updated : ", f2.PCT_CROP.values) print("Updating PCT_NAT_PFT") - #print (f2.PCT_NAT_PFT) print(f2.PCT_NAT_PFT.values[0]) - #f2.PCT_NAT_PFT.values[0] = [[100.0]] print(f2.PCT_NAT_PFT[0].values) out_dir = args.out_dir @@ -729,13 +667,4 @@ def main(): print(f2.attrs) f2.to_netcdf(path=wfile, mode="w", format="NETCDF3_64BIT") - print( - "Successfully updated surface data file for neon site(" - + site_name - + "):\n - " - + wfile - ) - - -if __name__ == "__main__": - main() + print("Successfully updated surface data file for neon site(" + site_name + "):\n - " + wfile) diff --git a/python/ctsm/site_and_regional/neon_surf_wrapper.py b/python/ctsm/site_and_regional/neon_surf_wrapper.py new file mode 100755 index 0000000000..a2cc619a29 --- /dev/null +++ b/python/ctsm/site_and_regional/neon_surf_wrapper.py @@ -0,0 +1,227 @@ +#! /usr/bin/env python3 +""" +|------------------------------------------------------------------| +|--------------------- Instructions -----------------------------| +|------------------------------------------------------------------| +This script is a simple wrapper for neon sites that performs the +following: + 1) For neon sites, subset surface dataset from global dataset + (i.e. ./subset_data.py ) + 2) Download neon and update the created surface dataset + based on the downloaded neon data. + (i.e. modify_singlept_site_neon.py) + +Instructions for running using conda python environments: + +../../py_env_create +conda activate ctsm_pylib + +""" +# TODO +# Automatic downloading of missing files if they are missing +# -[ ] Download neon sites and dom pft file +# -[ ] Make sure verbose works for printing out commands running + +# Import libraries +from __future__ import print_function + +import os +import logging +import argparse +import subprocess +import tqdm +import pandas as pd + + +def get_parser(): + """ + Get parser object for this script. + """ + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + + parser.print_usage = parser.print_help + + parser.add_argument( + "-v", + "--verbose", + help="Verbose mode will print more information. ", + action="store_true", + dest="verbose", + default=False, + ) + + parser.add_argument( + "--16pft", + help="Create and/or modify 16-PFT surface datasets (e.g. for a FATES run) ", + action="store_true", + dest="pft_16", + default=False, + ) + + parser.add_argument( + "-m", + "--mixed", + help="Do not overwrite surface dataset to be just one dominant PFT at 100%", + action="store_true", + dest="mixed", + default=False, + ) + + return parser + + +def execute(command): + """ + Function for running a command on shell. + Args: + command (str): + command that we want to run. + Raises: + Error with the return code from shell. + """ + print("\n", " >> ", *command, "\n") + + try: + subprocess.check_call(command, stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT) + + except subprocess.CalledProcessError as e: + print(e) + + +def main(): + """ + Loop through neon sites and execute subset and modify commands + """ + args = get_parser().parse_args() + + if args.verbose: + logging.basicConfig(level=logging.DEBUG) + + neon_sites = pd.read_csv("neon_sites_dompft.csv") + + for i, row in tqdm.tqdm(neon_sites.iterrows()): + lat = row["Lat"] + lon = row["Lon"] + site = row["Site"] + pft = row["pft"] + clmsite = "1x1_NEON_" + site + print("Now processing site :", site) + + if args.mixed and args.pft_16: + # use surface dataset with 16 pfts, and don't overwrite with 100% 1 dominant PFT + # don't set crop flag + # don't set a dominant pft + subset_command = [ + "./subset_data", + "point", + "--lat", + str(lat), + "--lon", + str(lon), + "--site", + clmsite, + "--create-surface", + "--uniform-snowpack", + "--cap-saturation", + "--verbose", + "--overwrite", + ] + modify_command = [ + "./modify_singlept_site_neon", + "--neon_site", + site, + "--surf_dir", + "subset_data_single_point", + "--16pft", + ] + elif args.pft_16: + # use surface dataset with 16 pfts, but overwrite to 100% 1 dominant PFT + # don't set crop flag + # set dominant pft + subset_command = [ + "./subset_data", + "point", + "--lat", + str(lat), + "--lon", + str(lon), + "--site", + clmsite, + "--dompft", + str(pft), + "--create-surface", + "--uniform-snowpack", + "--cap-saturation", + "--verbose", + "--overwrite", + ] + modify_command = [ + "./modify_singlept_site_neon", + "--neon_site", + site, + "--surf_dir", + "subset_data_single_point", + "--16pft", + ] + elif args.mixed: + # use surface dataset with 78 pfts, and don't overwrite with 100% 1 dominant PFT + # NOTE: FATES will currently not run with a 78-PFT surface dataset + # set crop flag + # don't set dominant pft + subset_command = [ + "./subset_data", + "point", + "--lat", + str(lat), + "--lon", + str(lon), + "--site", + clmsite, + "--crop", + "--create-surface", + "--uniform-snowpack", + "--cap-saturation", + "--verbose", + "--overwrite", + ] + modify_command = [ + "./modify_singlept_site_neon", + "--neon_site", + site, + "--surf_dir", + "subset_data_single_point", + ] + else: + # use surface dataset with 78 pfts, and overwrite to 100% 1 dominant PFT + # NOTE: FATES will currently not run with a 78-PFT surface dataset + # set crop flag + # set dominant pft + subset_command = [ + "./subset_data", + "point", + "--lat", + str(lat), + "--lon", + str(lon), + "--site", + clmsite, + "--crop", + "--dompft", + str(pft), + "--create-surface", + "--uniform-snowpack", + "--cap-saturation", + "--verbose", + "--overwrite", + ] + modify_command = [ + "./modify_singlept_site_neon", + "--neon_site", + site, + "--surf_dir", + "subset_data_single_point", + ] + execute(subset_command) + execute(modify_command) diff --git a/tools/site_and_regional/run_neon.py b/python/ctsm/site_and_regional/run_neon.py similarity index 84% rename from tools/site_and_regional/run_neon.py rename to python/ctsm/site_and_regional/run_neon.py index 84c00715fb..104e325617 100755 --- a/tools/site_and_regional/run_neon.py +++ b/python/ctsm/site_and_regional/run_neon.py @@ -51,37 +51,32 @@ # Import libraries - +import argparse +import datetime +import glob +import logging import os +import re +import shutil import sys import time -import shutil -import logging -import requests -import argparse -import re -import subprocess import pandas as pd -import glob -import datetime -from getpass import getuser # Get the ctsm util tools and then the cime tools. -_CTSM_PYTHON = os.path.abspath( - os.path.join(os.path.dirname(__file__), "..", "..", "python") -) +_CTSM_PYTHON = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "python")) sys.path.insert(1, _CTSM_PYTHON) -from ctsm import add_cime_to_path +from CIME import build +from CIME.case import Case +from CIME.utils import safe_copy, expect, symlink_force + from ctsm.path_utils import path_to_ctsm_root +from ctsm.utils import parse_isoduration from ctsm.download_utils import download_file -import CIME.build as build +from ctsm import add_cime_to_path + from standard_script_setup import * -from CIME.case import Case -from CIME.utils import safe_copy, expect, symlink_force, run_cmd_no_fail -from argparse import RawTextHelpFormatter -from CIME.locked_files import lock_file, unlock_file logger = logging.getLogger(__name__) @@ -194,7 +189,7 @@ def get_parser(args, description, valid_neon_sites): ) parser.add_argument( - "--prism", + "--prism", help=""" Uses the PRISM reanaylsis precipitation data for the site instead of the NEON data (only available over Continental US) @@ -205,9 +200,8 @@ def get_parser(args, description, valid_neon_sites): default=False, ) - parser.add_argument( - "--experiment", + "--experiment", help=""" Appends the case name with string for model experiment """, @@ -274,12 +268,11 @@ def get_parser(args, description, valid_neon_sites): """, action="store", dest="user_version", - required = False, - type = str, - choices= ['v1','v2','v3'], + required=False, + type=str, + choices=["v1", "v2", "v3"], ) - args = CIME.utils.parse_args_and_handle_standard_logging_options(args, parser) if "all" in args.neon_sites: @@ -299,7 +292,8 @@ def get_parser(args, description, valid_neon_sites): elif args.run_type == "postad": run_length = "100Y" else: - # The transient run length is set by cdeps atm buildnml to the last date of the available tower data + # The transient run length is set by cdeps atm buildnml to + # the last date of the available tower data # this value is not used run_length = "4Y" else: @@ -310,7 +304,8 @@ def get_parser(args, description, valid_neon_sites): if args.base_case_root: base_case_root = os.path.abspath(args.base_case_root) - # Reduce output level for this script unless --debug or --verbose is provided on the command line + # Reduce output level for this script unless --debug or + # --verbose is provided on the command line if not args.debug and not args.verbose: root_logger = logging.getLogger() root_logger.setLevel(logging.WARN) @@ -332,31 +327,6 @@ def get_parser(args, description, valid_neon_sites): ) -def get_isosplit(s, split): - if split in s: - n, s = s.split(split) - else: - n = 0 - return n, s - - -def parse_isoduration(s): - """ - simple ISO 8601 duration parser, does not account for leap years and assumes 30 day months - """ - # Remove prefix - s = s.split("P")[-1] - - # Step through letter dividers - years, s = get_isosplit(s, "Y") - months, s = get_isosplit(s, "M") - days, s = get_isosplit(s, "D") - - # Convert all to timedelta - dt = datetime.timedelta(days=int(days) + 365 * int(years) + 30 * int(months)) - return int(dt.total_seconds() / 86400) - - class NeonSite: """ A class for encapsulating neon sites. @@ -381,9 +351,7 @@ def __init__(self, name, start_year, end_year, start_month, end_month, finidat): def __str__(self): return ( - str(self.__class__) - + "\n" - + "\n".join((str(item) + " = " for item in (self.__dict__))) + str(self.__class__) + "\n" + "\n".join((str(item) + " = " for item in (self.__dict__))) ) def build_base_case( @@ -408,9 +376,7 @@ def build_base_case( """ print("---- building a base case -------") self.base_case_root = output_root - user_mods_dirs = [ - os.path.join(cesmroot, "cime_config", "usermods_dirs", "NEON", self.name) - ] + user_mods_dirs = [os.path.join(cesmroot, "cime_config", "usermods_dirs", "NEON", self.name)] if not output_root: output_root = os.getcwd() case_path = os.path.join(output_root, self.name) @@ -449,9 +415,15 @@ def build_base_case( existingcompname = case.get_value("COMPSET") match = re.search("^HIST", existingcompname, flags=re.IGNORECASE) if re.search("^HIST", compset, flags=re.IGNORECASE) is None: - expect( match == None, "Existing base case is a historical type and should not be -- rerun with the --orverwrite option" ) + expect( + match is None, + "Existing base case is a historical type and should not be -- rerun with the --overwrite option", + ) else: - expect( match != None, "Existing base case should be a historical type and is not -- rerun with the --orverwrite option" ) + expect( + match is not None, + "Existing base case should be a historical type and is not -- rerun with the --overwrite option", + ) # reset the case case.case_setup(reset=True) case_path = case.get_value("CASEROOT") @@ -471,25 +443,27 @@ def build_base_case( return case_path def diff_month(self): + """ + Determine difference between two dates in months + """ d1 = datetime.datetime(self.end_year, self.end_month, 1) d2 = datetime.datetime(self.start_year, self.start_month, 1) return (d1.year - d2.year) * 12 + d1.month - d2.month - + def get_batch_query(self, case): """ - Function for querying the batch queue query command for a case, depending on the - user's batch system. + Function for querying the batch queue query command for a case, depending on the + user's batch system. Args: case: case object """ - + if case.get_value("BATCH_SYSTEM") == "none": - return "none" - else: - return case.get_value("batch_query") - + return "none" + return case.get_value("batch_query") + def run_case( self, base_case_root, @@ -503,10 +477,34 @@ def run_case( rerun=False, experiment=False, ): + """ + Run case. + + Args: + self + base_case_root: str, opt + file path of base case + run_type: str, opt + transient, post_ad, or ad case, default transient + prism: bool, opt + if True, use PRISM precipitation, default False + run_length: str, opt + length of run, default '4Y' + user_version: str, opt + default 'latest' + overwrite: bool, opt + default False + setup_only: bool, opt + default False; if True, set up but do not run case + no_batch: bool, opt + default False + rerun: bool, opt + default False + experiment: str, opt + name of experiment, default False + """ user_mods_dirs = [ - os.path.join( - self.cesmroot, "cime_config", "usermods_dirs", "NEON", self.name - ) + os.path.join(self.cesmroot, "cime_config", "usermods_dirs", "NEON", self.name) ] expect( os.path.isdir(base_case_root), @@ -516,15 +514,13 @@ def run_case( if user_version: version = user_version else: - version = 'latest' + version = "latest" - print ("using this version:", version) + print("using this version:", version) - if experiment != None: + if experiment is not None: self.name = self.name + "." + experiment - case_root = os.path.abspath( - os.path.join(base_case_root, "..", self.name + "." + run_type) - ) + case_root = os.path.abspath(os.path.join(base_case_root, "..", self.name + "." + run_type)) rundir = None if os.path.isdir(case_root): @@ -538,15 +534,17 @@ def run_case( existingcompname = case.get_value("COMPSET") match = re.search("^HIST", existingcompname, flags=re.IGNORECASE) if re.search("^HIST", compset, flags=re.IGNORECASE) is None: - expect( match == None, "Existing base case is a historical type and should not be -- rerun with the --orverwrite option" ) + expect( + match is None, + "Existing base case is a historical type and should not be -- rerun with the --overwrite option", + ) else: - expect( match != None, "Existing base case should be a historical type and is not -- rerun with the --orverwrite option" ) - if os.path.isfile(os.path.join(rundir, "ESMF_Profile.summary")): - print( - "Case {} appears to be complete, not rerunning.".format( - case_root - ) + expect( + match is not None, + "Existing base case should be a historical type and is not -- rerun with the --overwrite option", ) + if os.path.isfile(os.path.join(rundir, "ESMF_Profile.summary")): + print("Case {} appears to be complete, not rerunning.".format(case_root)) elif not setup_only: print("Resubmitting case {}".format(case_root)) case.submit(no_batch=no_batch) @@ -557,17 +555,13 @@ def run_case( print(f"Use {batch_query} to check its run status") return else: - logger.warning( - "Case already exists in {}, not overwritting.".format(case_root) - ) + logger.warning("Case already exists in {}, not overwritting.".format(case_root)) return if run_type == "postad": adcase_root = case_root.replace(".postad", ".ad") if not os.path.isdir(adcase_root): - logger.warning( - "postad requested but no ad case found in {}".format(adcase_root) - ) + logger.warning("postad requested but no ad case found in {}".format(adcase_root)) return if not os.path.isdir(case_root): @@ -580,15 +574,14 @@ def run_case( # that the shell_commands file is copied, as well as taking care of the DATM inputs. # See https://github.com/ESCOMP/CTSM/pull/1872#pullrequestreview-1169407493 # - basecase.create_clone( - case_root, keepexe=True, user_mods_dirs=user_mods_dirs - ) + basecase.create_clone(case_root, keepexe=True, user_mods_dirs=user_mods_dirs) with Case(case_root, read_only=False) as case: if run_type != "transient": - # in order to avoid the complication of leap years we always set the run_length in units of days. - case.set_value("STOP_OPTION", "ndays") - case.set_value("REST_OPTION", "end") + # in order to avoid the complication of leap years, + # we always set the run_length in units of days. + case.set_value("STOP_OPTION", "ndays") + case.set_value("REST_OPTION", "end") case.set_value("CONTINUE_RUN", False) case.set_value("NEONVERSION", version) if prism: @@ -639,6 +632,9 @@ def run_case( print(f"Use {batch_query} to check its run status") def set_ref_case(self, case): + """ + Set an existing case as the reference case, eg for use with spinup. + """ rundir = case.get_value("RUNDIR") case_root = case.get_value("CASEROOT") if case_root.endswith(".postad"): @@ -660,9 +656,7 @@ def set_ref_case(self, case): case.set_value("RUN_REFDIR", refrundir) case.set_value("RUN_REFCASE", os.path.basename(ref_case_root)) refdate = None - for reffile in glob.iglob( - refrundir + "/{}{}.clm2.r.*.nc".format(self.name, root) - ): + for reffile in glob.iglob(refrundir + "/{}{}.clm2.r.*.nc".format(self.name, root)): m = re.search("(\d\d\d\d-\d\d-\d\d)-\d\d\d\d\d.nc", reffile) if m: refdate = m.group(1) @@ -677,9 +671,7 @@ def set_ref_case(self, case): if not os.path.isdir(os.path.join(rundir, "inputdata")) and os.path.isdir( os.path.join(refrundir, "inputdata") ): - symlink_force( - os.path.join(refrundir, "inputdata"), os.path.join(rundir, "inputdata") - ) + symlink_force(os.path.join(refrundir, "inputdata"), os.path.join(rundir, "inputdata")) case.set_value("RUN_REFDATE", refdate) if case_root.endswith(".postad"): @@ -688,14 +680,16 @@ def set_ref_case(self, case): return True def modify_user_nl(self, case_root, run_type, rundir): + """ + Modify user namelist. If transient, include finidat in user_nl; + Otherwise, adjust user_nl to include different mfilt, nhtfrq, and variables in hist_fincl1. + """ user_nl_fname = os.path.join(case_root, "user_nl_clm") user_nl_lines = None if run_type == "transient": if self.finidat: user_nl_lines = [ - "finidat = '{}/inputdata/lnd/ctsm/initdata/{}'".format( - rundir, self.finidat - ) + "finidat = '{}/inputdata/lnd/ctsm/initdata/{}'".format(rundir, self.finidat) ] else: user_nl_lines = [ @@ -766,16 +760,15 @@ def parse_neon_listing(listing_file, valid_neon_sites): tmp_df = tmp_df[tmp_df[9].str.contains("\d\d\d\d-\d\d.nc")] # -- find all the data versions - versions = tmp_df[7].unique() - #print ("all versions available for ", site_name,":", *versions) + # versions = tmp_df[7].unique() + # print ("all versions available for ", site_name,":", *versions) latest_version = tmp_df[7].iloc[-1] - #print ("latests version available for ", site_name,":", latest_version) + # print ("latests version available for ", site_name,":", latest_version) tmp_df = tmp_df[tmp_df[7].str.contains(latest_version)] # -- remove .nc from the file names tmp_df[9] = tmp_df[9].str.replace(".nc", "", regex=False) - tmp_df2 = tmp_df[9].str.split("-", expand=True) # ignore any prefix in file name and just get year @@ -800,9 +793,7 @@ def parse_neon_listing(listing_file, valid_neon_sites): if site_name in line: finidat = line.split(",")[0].split("/")[-1] - neon_site = NeonSite( - site_name, start_year, end_year, start_month, end_month, finidat - ) + neon_site = NeonSite(site_name, start_year, end_year, start_month, end_month, finidat) logger.debug(neon_site) available_list.append(neon_site) @@ -810,6 +801,10 @@ def parse_neon_listing(listing_file, valid_neon_sites): def main(description): + """ + Determine valid neon sites. Make an output directory if it does not exist. + Loop through requested sites and run CTSM at that site. + """ cesmroot = path_to_ctsm_root() # Get the list of supported neon sites from usermods valid_neon_sites = glob.glob( @@ -847,9 +842,9 @@ def main(description): res = "CLM_USRDAT" if run_type == "transient": - compset = "IHist1PtClm51Bgc" + compset = "IHist1PtClm51Bgc" else: - compset = "I1PtClm51Bgc" + compset = "I1PtClm51Bgc" # -- Looping over neon sites @@ -875,8 +870,3 @@ def main(description): rerun, experiment, ) - - - -if __name__ == "__main__": - main(__doc__) diff --git a/python/ctsm/test/test_sys_modify_singlept_site_neon.py b/python/ctsm/test/test_sys_modify_singlept_site_neon.py new file mode 100755 index 0000000000..b5ded30399 --- /dev/null +++ b/python/ctsm/test/test_sys_modify_singlept_site_neon.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 + +""" +System tests for modify_singlept_site_neon.py +""" + +import os +import unittest +import tempfile +import shutil +import sys + +from ctsm.path_utils import path_to_ctsm_root +from ctsm import unit_testing +from ctsm.site_and_regional.modify_singlept_site_neon import main, get_parser + +# Allow test names that pylint doesn't like; otherwise hard to make them +# readable +# pylint: disable=invalid-name + + +class TestSysModifySingleptSiteNeon(unittest.TestCase): + """System tests for modify_singlept_site_neon""" + + def setUp(self): + """ + Make /_tempdir for use by these tests. + Check tempdir for history files + """ + self._tempdir = tempfile.mkdtemp() + testinputs_path = os.path.join(path_to_ctsm_root(), "python/ctsm/test/testinputs") + self._cfg_file_path = os.path.join( + testinputs_path, "modify_singlept_site_neon_opt_sections.cfg" + ) + + def tearDown(self): + """ + Remove temporary directory + """ + shutil.rmtree(self._tempdir, ignore_errors=True) + + def test_modify_site(self): + """ + Test modifying a singple point site. + This test currently checks that the run fails due to dir structure + + TODO: The primary items to test here are the following: + 1) Fields are overwritten with site-specific data for neon sites + 2) Downloaded data is used in surface dataset + 3) Check specific fields listed in update_metadata for correct output + 4) Check that a netcdf with correct formatting is created + """ + sys.argv = [ + "modify_singlept_site_neon", + "--neon_site", + path_to_ctsm_root() + "/ctsm/cime_config/usermods_dirs/NEON/ABBY", + ] + # TODO: the above requires a full path instead of site name + # because of how run_neon is configured. + # This needs to be fixed in run_neon. + parser = get_parser() + with self.assertRaises(SystemExit): + print( + """This should currently fail due to directory structure in run_neon + and the directory structure listed in sys.argv""" + ) + main() + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/test_sys_run_neon.py b/python/ctsm/test/test_sys_run_neon.py new file mode 100755 index 0000000000..b6814ee2bc --- /dev/null +++ b/python/ctsm/test/test_sys_run_neon.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 + +"""System tests for run_neon + +""" + +import glob +import os +import unittest +import tempfile +import shutil +import sys + +from ctsm import unit_testing +from ctsm.site_and_regional.run_neon import main, get_parser +from ctsm.path_utils import path_to_ctsm_root + +# Allow test names that pylint doesn't like; otherwise hard to make them +# readable +# pylint: disable=invalid-name + + +class TestSysRunNeon(unittest.TestCase): + """System tests for run_neon""" + + def setUp(self): + """ + Make /_tempdir for use by these tests. + Check tempdir for history files + """ + self._tempdir = tempfile.mkdtemp() + os.chdir(self._tempdir) # cd to tempdir + + def tearDown(self): + """ + Remove temporary directory + """ + shutil.rmtree(self._tempdir, ignore_errors=True) + + def test_one_site(self): + """ + This test specifies a site to run + Run the tool, check that file structure is set up correctly + """ + + # run the run_neon tool + sys.argv = [ + os.path.join(path_to_ctsm_root(), "tools", "site_and_regional", "run_neon"), + "--neon-sites", + "BART", + "--setup-only", + "--output-root", + self._tempdir, + ] + valid_neon_sites = ["ABBY", "OSBS", "BART"] + parser = get_parser(sys.argv, "description_for_parser", valid_neon_sites) + main("") + + # assert that BART directories were created during setup + self.assertTrue("BART" in glob.glob(self._tempdir + "/BART*")[0]) + + # TODO: Would also be useful to test the following items: + # It might be good to ensure the log files are working as expected? + # Test running transient, ad and post ad cases. + # Test use of base case root. + # Test for using prism? + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/test_unit_iso_utils.py b/python/ctsm/test/test_unit_iso_utils.py new file mode 100755 index 0000000000..c58beef52e --- /dev/null +++ b/python/ctsm/test/test_unit_iso_utils.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python3 + +"""Unit tests for the iso functions in utils +""" + +import unittest + +from ctsm import unit_testing +from ctsm.utils import parse_isoduration, get_isosplit + +# Allow names that pylint doesn't like, because otherwise I find it hard +# to make readable unit test names +# pylint: disable=invalid-name + + +class TestIsoUtils(unittest.TestCase): + """Tests of iso functions in utils""" + + def test_iso_split_for_Year(self): + """ + Tests the get_isosplit function for a strings with Years + """ + iso_string = "0Y" + self.assertEqual(get_isosplit(iso_string, "Y"), ("0", "")) + iso_string = "1Y" + self.assertEqual(get_isosplit(iso_string, "Y"), ("1", "")) + iso_string = "4Y" + self.assertEqual(get_isosplit(iso_string, "Y"), ("4", "")) + iso_string = "100Y" + self.assertEqual(get_isosplit(iso_string, "Y"), ("100", "")) + iso_string = "999999Y" + self.assertEqual(get_isosplit(iso_string, "Y"), ("999999", "")) + + def test_parse_isoduration_for_Years(self): + """ + Tests the parse_isoduration function for iso strings with Years + """ + days_in_year = 365 + iso_string = "0Y" + self.assertEqual(parse_isoduration(iso_string), 0) + iso_string = "1Y" + self.assertEqual(parse_isoduration(iso_string), days_in_year) + iso_string = "4Y" + self.assertEqual(parse_isoduration(iso_string), 4 * days_in_year) + iso_string = "100Y" + self.assertEqual(parse_isoduration(iso_string), 100 * days_in_year) + iso_string = "999999Y" + self.assertEqual(parse_isoduration(iso_string), 999999 * days_in_year) + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/test_unit_modify_singlept_site_neon.py b/python/ctsm/test/test_unit_modify_singlept_site_neon.py new file mode 100755 index 0000000000..ecd96357b3 --- /dev/null +++ b/python/ctsm/test/test_unit_modify_singlept_site_neon.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python3 +""" +Unit tests for modify_singlept_site_neon + +You can run this by: + python -m unittest test_unit_modify_singlept_site_neon.py +""" + +import os +import shutil +import sys +import tempfile +import unittest +from datetime import date +import xarray as xr + +# -- add python/ctsm to path (needed if we want to run the test 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.path_utils import path_to_ctsm_root + +# pylint: disable=wrong-import-position +from ctsm import unit_testing +from ctsm.site_and_regional.modify_singlept_site_neon import ( + get_neon, + find_surffile, + update_metadata, + update_time_tag, + check_neon_time, +) + +# pylint: disable=invalid-name + + +class TestModifySingleptSiteNeon(unittest.TestCase): + """ + Basic class for testing modify_singlept_site_neon.py. + """ + + def setUp(self): + """ + Make /_tempdir for use by these tests. + Check tempdir for history files + """ + self._tempdir = tempfile.mkdtemp() + + def tearDown(self): + """ + Remove temporary directory + """ + shutil.rmtree(self._tempdir, ignore_errors=True) + + def test_get_neon(self): + """ + Test to see if neon data for valid site name is found + """ + site_name = "ABBY" + neon_dir = self._tempdir + file = get_neon(neon_dir, site_name) + self.assertEqual(file.split("/")[-1][:4], "ABBY", "CSV file did not download as expected") + + def test_get_neon_false_site(self): + """ + Test to see if neon data for invalid site name is found + """ + site_name = "INVALID_SITE" + neon_dir = self._tempdir + with self.assertRaises(SystemExit): + get_neon(neon_dir, site_name) + + def test_find_surffile(self): + """ + Test that surface file does not exist in tempdir and raises system exit error + """ + surf_dir = self._tempdir + site_name = "BART" + pft_16 = True + with self.assertRaises(SystemExit): + find_surffile(surf_dir, site_name, pft_16) + + def test_find_soil_structure(self): + """ + Test to ensure that correct attributes are found for find_soil_structure. + soil_texture_raw_data_file_name should be found, and test should go through sysexit. + """ + surf_file_name = "surfdata_1x1_mexicocityMEX_hist_16pfts_Irrig_CMIP6_simyr2000_c221206.nc" + surf_file = os.path.join( + path_to_ctsm_root(), + "python/ctsm/test/testinputs/", + surf_file_name, + ) + f1 = xr.open_dataset(surf_file) + self.assertEqual( + f1.attrs["Soil_texture_raw_data_file_name"], + "mksrf_soitex.10level.c010119.nc", + "did not retrieve expected surface soil texture filename from surf file", + ) + + def test_update_metadata(self): + """ + Test to ensure that the file was updated today. + """ + surf_file = "surfdata_1x1_mexicocityMEX_hist_16pfts_Irrig_CMIP6_simyr2000_c221206.nc" + neon_file = "dummy_neon_file.nc" + zb_flag = True + f1 = xr.open_dataset( + os.path.join(path_to_ctsm_root(), "python/ctsm/test/testinputs/") + surf_file + ) + f2 = update_metadata(f1, surf_file, neon_file, zb_flag) + today = date.today() + today_string = today.strftime("%Y-%m-%d") + self.assertEqual(f2.attrs["Updated_on"], today_string, "File was not updated as expected") + + def test_update_time_tag(self): + """ + Test that file ending is updated + """ + self.assertEqual( + update_time_tag("test_YYMMDD.nc")[-9:-3], + date.today().strftime("%y%m%d"), + "File ending not as expected", + ) + + def test_check_neon_time(self): + """ + Test that dictionary containing last modified information is correctly downloaded + """ + previous_dir = os.getcwd() + os.chdir(self._tempdir) # cd to tempdir + last_abby_download = check_neon_time()[ + "https://storage.neonscience.org/neon-ncar/NEON/surf_files/v1/ABBY_surfaceData.csv" + ] + self.assertEqual( + len(last_abby_download), + 19, + "last ABBY download has unexpected date format or does not exist", + ) + # Note: this checks that data is not pulled from before 2021; + # we may want to update this occassionally, + # but in any case it confirms that the oldest data is not found + self.assertGreater( + int(last_abby_download[:4]), 2021, "ABBY download is older than expected" + ) + # change back to previous dir once listing.csv file is created in tempdir and test complete + os.chdir(previous_dir) + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/test_unit_neon_surf_wrapper.py b/python/ctsm/test/test_unit_neon_surf_wrapper.py new file mode 100755 index 0000000000..443af2079b --- /dev/null +++ b/python/ctsm/test/test_unit_neon_surf_wrapper.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 +""" +Unit tests for neon_surf_wrapper + +You can run this by: + python -m unittest test_unit_neon_surf_wrapper.py +""" + +import unittest +import os +import sys + +# -- add python/ctsm to path (needed if we want to run the test stand-alone) +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir) +sys.path.insert(1, _CTSM_PYTHON) + +# pylint: disable=wrong-import-position +from ctsm import unit_testing +from ctsm.site_and_regional.neon_surf_wrapper import get_parser + +# pylint: disable=invalid-name + + +class TestNeonSurfWrapper(unittest.TestCase): + """ + Basic class for testing neon_surf_wrapper.py. + """ + + def test_parser(self): + """ + Test that parser has same defaults as expected + """ + + self.assertEqual(get_parser().argument_default, None, "Parser not working as expected") + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/test_unit_run_neon.py b/python/ctsm/test/test_unit_run_neon.py new file mode 100755 index 0000000000..a35608e249 --- /dev/null +++ b/python/ctsm/test/test_unit_run_neon.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +""" +Unit tests for run_neon + +You can run this by: + python -m unittest test_unit_run_neon.py +""" + +import unittest +import tempfile +import shutil +import os +import sys + +# -- add python/ctsm to path (needed if we want to run the test stand-alone) +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir) +sys.path.insert(1, _CTSM_PYTHON) + +# pylint: disable=wrong-import-position +from ctsm import unit_testing +from ctsm.site_and_regional.run_neon import check_neon_listing + +# pylint: disable=invalid-name + + +class TestRunNeon(unittest.TestCase): + """ + Basic class for testing run_neon.py. + """ + + def setUp(self): + """ + Make /_tempdir for use by these tests. + """ + self._tempdir = tempfile.mkdtemp() + + def tearDown(self): + """ + Remove temporary directory + """ + shutil.rmtree(self._tempdir, ignore_errors=True) + + def test_check_neon_listing(self): + """ + Test that neon listing is available for valid sites + """ + valid_neon_sites = ["ABBY", "BART"] + previous_dir = os.getcwd() + os.chdir(self._tempdir) # cd to tempdir + available_list = check_neon_listing(valid_neon_sites) + self.assertEqual( + available_list[0].name, "ABBY", "available list of actual sites not as expected" + ) + self.assertEqual( + available_list[1].name, "BART", "available list of actual sites not as expected" + ) + # change to previous dir once listing.csv file is created in tempdir and test complete + os.chdir(previous_dir) + + def test_check_neon_listing_misspelled(self): + """ + Test that neon listing is not available for invalid sites + """ + valid_neon_sites = ["INVALID_SITE1", "INVALID_SITE2"] + previous_dir = os.getcwd() + os.chdir(self._tempdir) # cd to tempdir + available_list = check_neon_listing(valid_neon_sites) + self.assertEqual( + available_list, [], "available list of incorrect dummy site not as expected" + ) + # change to previous dir once listing.csv file is created in tempdir and test complete + os.chdir(previous_dir) + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/testinputs/README.md b/python/ctsm/test/testinputs/README.md index 45451b53e1..1775de48b7 100644 --- a/python/ctsm/test/testinputs/README.md +++ b/python/ctsm/test/testinputs/README.md @@ -6,7 +6,7 @@ Installing Git LFS on your machine is a two-step process; step (1) needs to be done once per machine, and step (2) needs to be done once per user: 1. Install the Git LFS tool: Follow the instructions on the [Git LFS page](https://git-lfs.github.com/) for installing Git LFS on your platform. - - On cheyenne, Git LFS is already available as long as you are using a git + - On cheyenne and casper, Git LFS is already available as long as you are using a git module rather than the default system-level git. So just make sure that you are always using git via a git module (`module load git`). - On a Mac using homebrew, this can be done with `brew install git-lfs`. diff --git a/python/ctsm/utils.py b/python/ctsm/utils.py index 42444e32c5..8578ea860c 100644 --- a/python/ctsm/utils.py +++ b/python/ctsm/utils.py @@ -7,7 +7,7 @@ import re import pdb -from datetime import date +from datetime import date, timedelta from getpass import getuser from ctsm.git_utils import get_ctsm_git_short_hash @@ -189,3 +189,33 @@ def write_output(file, file_in, file_out, file_type): file.to_netcdf(path=file_out, mode="w", format="NETCDF3_64BIT") logger.info("Successfully created: %s", file_out) file.close() + + +def get_isosplit(iso_string, split): + """ + Split a string (iso_string) by the character sent in from split + Returns the number for that character split + Only used by parse_isoduration + """ + if split in iso_string: + num, iso_string = iso_string.split(split) + else: + num = 0 + return num, iso_string + + +def parse_isoduration(iso_string): + """ + simple ISO 8601 duration parser, does not account for leap years and assumes 30 day months + """ + # Remove prefix + iso_string = iso_string.split("P")[-1] + + # Step through letter dividers + years, iso_string = get_isosplit(iso_string, "Y") + months, iso_string = get_isosplit(iso_string, "M") + days, iso_string = get_isosplit(iso_string, "D") + + # Convert all to timedelta + delta_t = timedelta(days=int(days) + 365 * int(years) + 30 * int(months)) + return int(delta_t.total_seconds() / 86400) diff --git a/src/biogeochem/CNVegCarbonFluxType.F90 b/src/biogeochem/CNVegCarbonFluxType.F90 index c7aa3469e2..298e7b3053 100644 --- a/src/biogeochem/CNVegCarbonFluxType.F90 +++ b/src/biogeochem/CNVegCarbonFluxType.F90 @@ -3729,8 +3729,8 @@ subroutine RestartBulkOnly ( this, bounds, ncid, flag ) ! BACKWARDS_COMPATIBILITY(wjs/ssr, 2022-06-10) See note in CallRestartvarDimOK() if (CallRestartvarDimOK(ncid, flag, 'mxharvests')) then do k = repr_grain_min, repr_grain_max - data2dptr => this%repr_grainc_to_food_perharv_patch(:,:,k) ! e.g., grainc_to_food_perharv + data2dptr => this%repr_grainc_to_food_perharv_patch(:,:,k) varname = get_repr_rest_fname(k)//'c_to_food_perharv' call restartvar(ncid=ncid, flag=flag, varname=varname, & xtype=ncd_double, & @@ -3742,12 +3742,26 @@ subroutine RestartBulkOnly ( this, bounds, ncid, flag ) readvar=readvar, & scale_by_thickness=.false., & interpinic_flag='interp', data=data2dptr) + + ! e.g., grainc_to_seed_perharv + data2dptr => this%repr_grainc_to_seed_perharv_patch(:,:,k) + varname = get_repr_rest_fname(k)//'c_to_seed_perharv' + call restartvar(ncid=ncid, flag=flag, varname=varname, & + xtype=ncd_double, & + dim1name='pft', & + dim2name='mxharvests', & + switchdim=.true., & + long_name=get_repr_longname(k)//' C to seed per harvest; should only be output annually', & + units='gC/m2', & + readvar=readvar, & + scale_by_thickness=.false., & + interpinic_flag='interp', data=data2dptr) end do end if do k = repr_grain_min, repr_grain_max - data1dptr => this%repr_grainc_to_food_thisyr_patch(:,k) ! e.g., grainc_to_food_thisyr + data1dptr => this%repr_grainc_to_food_thisyr_patch(:,k) varname = get_repr_rest_fname(k)//'c_to_food_thisyr' call restartvar(ncid=ncid, flag=flag, varname=varname, & xtype=ncd_double, & @@ -3755,8 +3769,9 @@ subroutine RestartBulkOnly ( this, bounds, ncid, flag ) long_name=get_repr_longname(k)//' C to food per calendar year; should only be output annually', & units='gC/m2', & interpinic_flag='interp', readvar=readvar, data=data1dptr) - data1dptr => this%repr_grainc_to_seed_thisyr_patch(:,k) + ! e.g., grainc_to_seed_thisyr + data1dptr => this%repr_grainc_to_seed_thisyr_patch(:,k) varname = get_repr_rest_fname(k)//'c_to_seed_thisyr' call restartvar(ncid=ncid, flag=flag, varname=varname, & xtype=ncd_double, & diff --git a/src/biogeochem/CropType.F90 b/src/biogeochem/CropType.F90 index 29c4717ab3..e8de7059c0 100644 --- a/src/biogeochem/CropType.F90 +++ b/src/biogeochem/CropType.F90 @@ -655,6 +655,16 @@ subroutine Restart(this, bounds, ncid, cnveg_state_inst, flag) long_name='crop sowing dates for this patch this year', units='day of year', & scale_by_thickness=.false., & interpinic_flag='interp', readvar=readvar, data=this%sdates_thisyr_patch) + call restartvar(ncid=ncid, flag=flag, varname='swindow_starts_thisyr_patch', xtype=ncd_double, & + dim1name='pft', dim2name='mxsowings', switchdim=.true., & + long_name='sowing window start dates for this patch this year', units='day of year', & + scale_by_thickness=.false., & + interpinic_flag='interp', readvar=readvar, data=this%swindow_starts_thisyr_patch) + call restartvar(ncid=ncid, flag=flag, varname='swindow_ends_thisyr_patch', xtype=ncd_double, & + dim1name='pft', dim2name='mxsowings', switchdim=.true., & + long_name='sowing window end dates for this patch this year', units='day of year', & + scale_by_thickness=.false., & + interpinic_flag='interp', readvar=readvar, data=this%swindow_ends_thisyr_patch) ! Fill variable(s) derived from read-in variable(s) if (flag == 'read' .and. readvar) then do p = bounds%begp,bounds%endp diff --git a/src/biogeophys/SoilTemperatureMod.F90 b/src/biogeophys/SoilTemperatureMod.F90 index b868224a60..0dc8876d24 100644 --- a/src/biogeophys/SoilTemperatureMod.F90 +++ b/src/biogeophys/SoilTemperatureMod.F90 @@ -47,7 +47,7 @@ module SoilTemperatureMod ! o The thermal conductivity of soil is computed from ! the algorithm of Johansen (as reported by Farouki 1981), and the ! conductivity of snow is from the formulation used in - ! SNTHERM (Jordan 1991). + ! Sturm (1997) or Jordan (1991) p. 18 depending on namelist option. ! o Boundary conditions: ! F = Rnet - Hg - LEg (top), F= 0 (base of the soil column). ! o Soil / snow temperature is predicted from heat conduction @@ -100,7 +100,7 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter ! o The thermal conductivity of soil is computed from ! the algorithm of Johansen (as reported by Farouki 1981), and the ! conductivity of snow is from the formulation used in - ! SNTHERM (Jordan 1991). + ! Sturm (1997) or Jordan (1991) p. 18 depending on namelist option. ! o Boundary conditions: ! F = Rnet - Hg - LEg (top), F= 0 (base of the soil column). ! o Soil / snow temperature is predicted from heat conduction @@ -611,18 +611,20 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter ! ! (2) The thermal conductivity of soil is computed from the algorithm of ! Johansen (as reported by Farouki 1981), and of snow is from the - ! formulation used in SNTHERM (Jordan 1991). + ! formulation used in Sturm (1997) or Jordan (1991) p. 18 depending on + ! namelist option. ! The thermal conductivities at the interfaces between two neighboring ! layers (j, j+1) are derived from an assumption that the flux across ! the interface is equal to that from the node j to the interface and the ! flux from the interface to the node j+1. ! ! !USES: + use shr_log_mod , only : errMsg => shr_log_errMsg use clm_varpar , only : nlevsno, nlevgrnd, nlevurb, nlevsoi, nlevmaxurbgrnd use clm_varcon , only : denh2o, denice, tfrz, tkwat, tkice, tkair, cpice, cpliq, thk_bedrock, csol_bedrock use landunit_varcon , only : istice, istwet use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv - use clm_varctl , only : iulog + use clm_varctl , only : iulog, snow_thermal_cond_method ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -647,6 +649,8 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter real(r8) :: fl ! volume fraction of liquid or unfrozen water to total water real(r8) :: satw ! relative total water content of soil. real(r8) :: zh2osfc + + character(len=*),parameter :: subname = 'SoilThermProp' !----------------------------------------------------------------------- call t_startf( 'SoilThermProp' ) @@ -734,11 +738,27 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter endif endif - ! Thermal conductivity of snow, which from Jordan (1991) pp. 18 + ! Thermal conductivity of snow ! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1 if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then bw(c,j) = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/(frac_sno(c)*dz(c,j)) - thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair) + select case (snow_thermal_cond_method) + case ('Jordan1991') + thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair) + case ('Sturm1997') + ! Implemented by Vicky Dutch (VRD), Nick Rutter, and + ! Leanne Wake (LMW) + ! https://tc.copernicus.org/articles/16/4201/2022/ + ! Code provided by Adrien Dams to Will Wieder + if (bw(c,j) <= 156) then !LMW or 0.156 ? + thk(c,j) = 0.023 + 0.234*(bw(c,j)/1000) !LMW - units changed by VRD + else !LMW + thk(c,j) = 0.138 - 1.01*(bw(c,j)/1000) +(3.233*((bw(c,j)/1000)*(bw(c,j)/1000))) ! LMW Sturm I think + end if + case default + write(iulog,*) subname//' ERROR: unknown snow_thermal_cond_method value: ', snow_thermal_cond_method + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select end if end do diff --git a/src/cpl/share_esmf/cropcalStreamMod.F90 b/src/cpl/share_esmf/cropcalStreamMod.F90 index 46696eeba9..0ea63f2c6d 100644 --- a/src/cpl/share_esmf/cropcalStreamMod.F90 +++ b/src/cpl/share_esmf/cropcalStreamMod.F90 @@ -7,7 +7,8 @@ module cropcalStreamMod ! Read crop calendars from streams ! ! !USES: - use ESMF + use ESMF , only : ESMF_LogFoundError, ESMF_LOGERR_PASSTHRU, ESMF_Finalize + use ESMF , only : ESMF_END_ABORT use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_CL, CS => shr_kind_CS use dshr_strdata_mod , only : shr_strdata_type use decompMod , only : bounds_type @@ -322,6 +323,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) integer :: n, g integer :: lsize integer :: rc + integer :: begp, endp real(r8), pointer :: dataptr1d_swindow_start(:) real(r8), pointer :: dataptr1d_swindow_end (:) real(r8), pointer :: dataptr1d_cultivar_gdds(:) @@ -342,6 +344,9 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! Place all data from each type into a temporary 2d array lsize = bounds%endg - bounds%begg + 1 + begp = bounds%begp + endp= bounds%endp + dayspyr = get_curr_days_per_year() ! Read prescribed sowing window start dates from input files @@ -397,17 +402,17 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) ! Ensure that, if mxsowings > 1, sowing windows are ordered such that ENDS are monotonically increasing. This is necessary because of how get_swindow() works. if (mxsowings > 1) then - if (any(ends(:,2:mxsowings) <= ends(:,1:mxsowings-1) .and. & - ends(:,2:mxsowings) >= 1)) then + if (any(ends(begp:endp,2:mxsowings) <= ends(begp:endp,1:mxsowings-1) .and. & + 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 < 1 .or. ends < 1)) then + if (any(starts(begp:endp,:) < 1 .or. ends(begp:endp,:) < 1)) then ! Fail if not allowing fallback to paramfile sowing windows - if ((.not. allow_invalid_swindow_inputs) .and. any(all(starts < 1, dim=2) .and. patch%wtgcell > 0._r8 .and. patch%itype >= npcropmin)) then + 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 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 @@ -422,7 +427,7 @@ subroutine cropcal_interp(bounds, num_pcropp, filter_pcropp, crop_inst) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Fail if a sowing window start date is prescribed without an end date (or vice versa) - else if (any((starts >= 1 .and. ends < 1) .or. (starts < 1 .and. ends >= 1))) then + else if (any((starts(begp:endp,:) >= 1 .and. ends(begp:endp,:) < 1) .or. (starts(begp:endp,:) < 1 .and. 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 diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index acefe4acda..6e89f0952e 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -221,6 +221,8 @@ module clm_varctl ! which snow cover fraction parameterization to use character(len=64), public :: snow_cover_fraction_method + ! which snow thermal conductivity parameterization to use + character(len=25), public :: snow_thermal_cond_method ! atmospheric CO2 molar ratio (by volume) (umol/mol) real(r8), public :: co2_ppmv = 355._r8 ! diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 5937e55b04..deb8c044d8 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -199,7 +199,8 @@ subroutine control_init(dtime) clump_pproc, & create_crop_landunit, nsegspc, co2_ppmv, & albice, soil_layerstruct_predefined, soil_layerstruct_userdefined, & - soil_layerstruct_userdefined_nlevsoi, use_subgrid_fluxes, snow_cover_fraction_method, & + soil_layerstruct_userdefined_nlevsoi, use_subgrid_fluxes, & + snow_thermal_cond_method, snow_cover_fraction_method, & irrigate, run_zero_weight_urban, all_active, & crop_fsat_equals_zero, for_testing_run_ncdiopio_tests, & for_testing_use_second_grain_pool, for_testing_use_repr_structure_pool, & @@ -850,6 +851,7 @@ subroutine control_spmd() ! physics variables call mpi_bcast (nsegspc, 1, MPI_INTEGER, 0, mpicom, ier) call mpi_bcast (use_subgrid_fluxes , 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (snow_thermal_cond_method, len(snow_thermal_cond_method), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (snow_cover_fraction_method , len(snow_cover_fraction_method), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (z0param_method , len(z0param_method), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (use_z0m_snowmelt, 1, MPI_LOGICAL, 0, mpicom, ier) diff --git a/tools/contrib/ssp_anomaly_forcing_smooth b/tools/contrib/ssp_anomaly_forcing_smooth index ae3d189a8b..362e47c67d 100755 --- a/tools/contrib/ssp_anomaly_forcing_smooth +++ b/tools/contrib/ssp_anomaly_forcing_smooth @@ -3,7 +3,7 @@ ssp_anomaly_forcing_smooth -Create anomoly forcing datasets for SSP scenarios that can be used by CESM datm model +Create anomaly forcing datasets for SSP scenarios that can be used by CESM datm model load proper modules first, i.e. @@ -21,6 +21,169 @@ import numpy as np import netCDF4 as netcdf4 +# Adds global attributes, returning hdir and fdir +def add_global_attributes(ds, historydate, histdir, sspdir, num_ens, climo_year, climo_base_nyrs, dpath, dfile, hist_yrstart, hist_yrend, ssp_yrstart, ssp_yrend, timetag): + ds.Created_on = timetag + + ds.title = "anomaly forcing data" + ds.note1 = ( + "Anomaly/scale factors calculated relative to " + + str(climo_year - (climo_base_nyrs - 1) / 2) + + "-" + + str(climo_year + (climo_base_nyrs - 1) / 2) + ) + ds.history = historydate + ": created by " + sys.argv[0] + stdout = os.popen("git describe") + ds.gitdescribe = stdout.read().rstrip() + ds.Source = "CMIP6 CESM simulations" + ds.Conventions = "CF-1.0" + ds.comment = ( + "Monthly scale factors for given SSP scenario compared to a climatology based on" + + " data centered on " + + str(climo_year) + + " over the range given in note1" + ) + ds.number_of_ensemble_members = str(num_ens) + ds.Created_by = getuser() + + for nens in range(num_ens): + hdir = dpath + histdir[nens] + dfile + fdir = dpath + sspdir[nens] + dfile + if nens == 0: + ds.Created_from_historical_dirs = hdir + ds.Created_from_scenario_dirs = fdir + else: + ds.Created_from_historical_dirs += ", " + hdir + ds.Created_from_scenario_dirs += ", " + fdir + + ds.History_years = str(hist_yrstart) + "," + str(hist_yrend) + ds.Scenario_years = str(ssp_yrstart) + "," + str(ssp_yrend) + ds.institution = "National Center for Atmospheric Research" + return hdir,fdir + + +def create_fill_latlon(ds, data, var_name): + + ds.createDimension(var_name, int(data.size)) + wl = ds.createVariable(var_name, np.float64, (var_name,)) + + if var_name == "lat": + wl.units = "degrees_north" + wl.long_name = "Latitude" + elif var_name == "lon": + wl.units = "degrees_east" + wl.long_name = "Longitude" + wl.mode = "time-invariant" + + wl[:] = data + + return ds + + +def create_fill_time(ds, time, ntime, ssp_time_units=None, ssp_time_longname=None, adj_time=False): + if ntime is not None: + ntime = int(ntime) + ds.createDimension("time", ntime) + + wtime = ds.createVariable("time", np.float64, ("time",)) + + if ssp_time_units is not None: + wtime.units = ssp_time_units + if ssp_time_longname is not None: + wtime.long_name = ssp_time_longname + wtime.calendar = "noleap" + + # adjust time to middle of month + if adj_time: + wtime_offset = 15 - time[0] + wtime[:] = time + wtime_offset + else: + wtime[:] = time + + return ds + + +def create_fill_ancillary_vars(ds, landfrac, landmask, area): + + wmask = ds.createVariable("landmask", np.int32, ("lat", "lon")) + warea = ds.createVariable("area", np.float64, ("lat", "lon")) + wfrac = ds.createVariable("landfrac", np.float64, ("lat", "lon")) + + warea.units = "km2" + wfrac.units = "unitless" + wmask.units = "unitless" + + warea.long_name = "Grid cell area" + wfrac.long_name = "Grid cell land fraction" + wmask.long_name = "Grid cell land mask" + + warea.mode = "time-invariant" + wfrac.mode = "time-invariant" + wmask.mode = "time-invariant" + + # write to file -------------------------------------------- + wmask[:, :] = landmask + wfrac[:, :] = landfrac + warea[:, :] = area + + return ds + + +def add_to_dataset(ds, var_name, data, units=None, mode=None, historical_source_files=None, scenario_source_files=None, long_name=None, cell_methods=None): + dims = ("time", "lat", "lon") + data_type = np.float64 + + wvar = ds.createVariable( + var_name, + data_type, + dims, + fill_value=data_type(1.0e36), + ) + + wvar[:, :, :] = data + + if units is not None: + wvar.units = units + if mode is not None: + wvar.mode = mode + if historical_source_files is not None: + wvar.historical_source_files = historical_source_files + if scenario_source_files is not None: + wvar.scenario_source_files = scenario_source_files + if long_name is not None: + wvar.long_name = long_name + if cell_methods is not None: + wvar.cell_methods = cell_methods + + return ds + + +def create_fill_forcing(ds, field_out, units, anomsf, field_out_wind, f, hdir, fdir, histfiles, sspfiles, long_name, anom_fld): + + historical_source_files = "".join(histfiles).replace(hdir, "") + scenario_source_files = "".join(sspfiles).replace(fdir, "") + mode = "time-dependent" + + if field_out[f] == "sfcWind": + long_name = str(long_name) + " U component " + anomsf[f] + var_name = field_out_wind[0] + data = anom_fld / np.sqrt(2) + else: + long_name = str(long_name) + " " + anomsf[f] + var_name = field_out[f] + data = anom_fld + # Was missing cell_methods attribute in original + ds = add_to_dataset(ds, var_name, data, units=units[f], mode=mode, historical_source_files=historical_source_files, scenario_source_files=scenario_source_files, long_name=long_name) + + if field_out[f] == "sfcWind": + long_name = long_name.replace("U component", "V component") + var_name = field_out_wind[1] + # Was missing mode attribute in original + ds = add_to_dataset(ds, var_name, data, units=units[f], historical_source_files=historical_source_files, scenario_source_files=scenario_source_files, long_name=long_name, cell_methods="time: mean") + + return ds + + parser = argparse.ArgumentParser(description="Create anomaly forcing") parser.add_argument( "sspnum", @@ -31,16 +194,25 @@ parser.add_argument( ) parser.add_argument( "--write_climo", + "--write-climo", help="write out climatology files and exit", action="store_true", default=False, ) parser.add_argument( "--print_ssps", + "--print-ssps", help="Just print out directory names and exit", action="store_true", default=False, ) +parser.add_argument( + "--output_dir", "--output-dir", + help="Top-level output directory (default: ./anomaly_forcing/). Sub-directory will be created for the selected scenario.", + type=str, + default=os.path.join(".", "anomaly_forcing"), +) + args = parser.parse_args() if args.sspnum == 0: @@ -48,7 +220,7 @@ if args.sspnum == 0: # ------------------------------------------------------- -print("Create anomoly forcing data that can be used by CTSM in CESM") +print("Create anomaly forcing data that can be used by CTSM in CESM") # Input and output directories make sure they exist datapath = "/glade/campaign/collections/cmip/CMIP6/timeseries-cmip6/" # Path on casper @@ -109,15 +281,13 @@ _v2 is just used for restart files that have been spatially interpolated """ -spath = "./" if os.path.exists(datapath): print("Input data directory:" + datapath) else: sys.exit("Could not find input directory: " + datapath) -if os.path.exists(spath): - print("Output data directory:" + spath) -else: - sys.exit("Could not find output directory: " + spath) +if not os.path.exists(args.output_dir): + os.makedirs(args.output_dir) +print("Output data directory:" + args.output_dir) # Settings to run with today = datetime.date.today() @@ -165,9 +335,6 @@ if args.print_ssps: sspnum = args.sspnum -# hist_case needed? -hist_case = "b.e21.BHIST.f09_g17.CMIP6-historical.010" - if sspnum == 1: # SSP1-26 ssptag = "SSP1-2.6" @@ -196,25 +363,15 @@ if num_ens != len(histdir): print("number of ensemble members not the same") sys.exit("number of members different") -# test w/ 1 ensemble member -num_ens = 3 - # Setup output directory -sspoutdir = "anomaly_forcing/CMIP6-" + ssptag +sspoutdir = "CMIP6-" + ssptag -outdir = spath + sspoutdir +outdir = os.path.join(args.output_dir, sspoutdir) if not os.path.exists(outdir): os.makedirs(outdir) print("Output specific data directory :" + outdir) -# historical files are split by 50 year periods; use last period -hist_suffix = ["200001-201412.nc"] # not standardized?! -# hist_suffix = ['-201412.nc'] -# projections are split 2015/2064 2065/2100 -ssp_suffix = ["201501-206412.nc", "206501-210012.nc"] -# ssp_suffix = ['-206412.nc','-210012.nc'] - climo_year = 2015 # ten years on either side (21 years total) climo_base_nyrs = 21 @@ -244,7 +401,11 @@ field_out_wind = ["uas", "vas"] nfields = len(field_in) +output_format = "NETCDF3_64BIT_DATA" + # -- Loop over forcing fields ------------------------------------ + + for f in range(nfields): # -- Loop over ensemble members ------------------------------ @@ -482,71 +643,31 @@ for f in range(nfields): if write_climo: # Use NetCDF4 format, because using older NetCDF formats are too slow w = netcdf4.Dataset( - outdir + field_out[f] + "_climo" + creationdate + ".nc", + os.path.join(outdir, field_out[f] + "_climo" + creationdate + ".nc"), "w", - format="NETCDF3_64BIT_DATA" + format=output_format, ) - w.createDimension("lat", int(nlat)) - w.createDimension("lon", int(nlon)) - w.createDimension("time", int(nmo)) - - wtime = w.createVariable("time", np.float64, ("time",)) - wlat = w.createVariable("lat", np.float64, ("lat",)) - wlon = w.createVariable("lon", np.float64, ("lon",)) - wvar = w.createVariable( - field_out[f], - np.float64, - ("time", "lat", "lon"), - fill_value=np.float64(1.0e36), - ) - wtime[ - :, - ] = time[0:12] - wlon[ - :, - ] = lon - wlat[ - :, - ] = lat - wvar[:, :, :] = climo + w = create_fill_latlon(w, lat, "lat") + w = create_fill_latlon(w, lon, "lon") + w = create_fill_time(w, time[0:12], nmo) + + add_to_dataset(w, field_out[f], climo) w.close() # Use NetCDF4 format, because using older NetCDF formats are too slow w = netcdf4.Dataset( - outdir + field_out[f] + "_smooth" + creationdate + ".nc", + os.path.join(outdir, field_out[f] + "_smooth" + creationdate + ".nc"), "w", - format="NETCDF3_64BIT_DATA" - ) - w.createDimension("lat", int(nlat)) - w.createDimension("lon", int(nlon)) - w.createDimension("time", int(tm)) - - wtime = w.createVariable("time", np.float64, ("time",)) - wlat = w.createVariable("lat", np.float64, ("lat",)) - wlon = w.createVariable("lon", np.float64, ("lon",)) - wvar = w.createVariable( - field_out[f], - np.float64, - ("time", "lat", "lon"), - fill_value=np.float64(1.0e36), + format=output_format, ) - wvar2 = w.createVariable( - "smooth_" + field_out[f], - np.float64, - ("time", "lat", "lon"), - fill_value=np.float64(1.0e36), - ) - - wtime[:] = time - wlon[ - :, - ] = lon - wlat[ - :, - ] = lat - wvar[:, :, :] = temp_fld - wvar2[:, :, :] = stemp_fld + w = create_fill_latlon(w, lat, "lat") + w = create_fill_latlon(w, lon, "lon") + w = create_fill_time(w, time, tm) + + add_to_dataset(w, field_out[f], temp_fld) + add_to_dataset(w, "smooth_" + field_out[f], stemp_fld) w.close() + print("Exit early after writing out climatology\n\n") sys.exit() @@ -556,9 +677,9 @@ for f in range(nfields): # Use NetCDF4 format, because using older NetCDF formats are too slow # Will need to convert to CDF5 format at the end, as we can't seem to # output in CDF5 format using netCDF4 python interfaces - outfilename = outdir + "/" + "af.allvars" + outfile_suffix + outfilename = os.path.join(outdir, "af.allvars" + outfile_suffix) print("Creating: " + outfilename) - outfile = netcdf4.Dataset(outfilename, "w", format="NETCDF3_64BIT_DATA") + outfile = netcdf4.Dataset(outfilename, "w", format=output_format) # creation date on the file command = 'date "+%Y/%m/%d"' @@ -566,148 +687,21 @@ for f in range(nfields): x = x2.communicate() timetag = x[0].decode("utf-8").strip() - outfile.Created_on = timetag + # Add global attributes and get hdir/fdir + hdir, fdir = add_global_attributes(outfile, historydate, histdir, sspdir, num_ens, climo_year, climo_base_nyrs, dpath, dfile, hist_yrstart, hist_yrend, ssp_yrstart, ssp_yrend, timetag) - outfile.title = "anomaly forcing data" - outfile.note1 = ( - "Anomaly/scale factors calculated relative to " - + str(climo_year - (climo_base_nyrs - 1) / 2) - + "-" - + str(climo_year + (climo_base_nyrs - 1) / 2) - ) - outfile.history = historydate + ": created by " + sys.argv[0] - stdout = os.popen("git describe") - outfile.gitdescribe = stdout.read().rstrip() - outfile.Source = "CMIP6 CESM simulations" - outfile.Conventions = "CF-1.0" - outfile.comment = ( - "Monthly scale factors for given SSP scenario compared to a climatology based on" - + " data centered on " - + str(climo_year) - + " over the range given in note1" - ) - outfile.number_of_ensemble_members = str(num_ens) - outfile.Created_by = getuser() - - for nens in range(num_ens): - hdir = dpath + histdir[nens] + dfile - fdir = dpath + sspdir[nens] + dfile - if nens == 0: - outfile.Created_from_historical_dirs = hdir - outfile.Created_from_scenario_dirs = fdir - else: - outfile.Created_from_historical_dirs += ", " + hdir - outfile.Created_from_scenario_dirs += ", " + fdir - - outfile.History_years = str(hist_yrstart) + "," + str(hist_yrend) - outfile.Scenario_years = str(ssp_yrstart) + "," + str(ssp_yrend) - outfile.institution = "National Center for Atmospheric Research" - - outfile.createDimension("lat", size=int(nlat)) - outfile.createDimension("lon", size=int(nlon)) - outfile.createDimension("time", None) - - wtime = outfile.createVariable("time", np.float64, ("time",)) - wlat = outfile.createVariable("lat", np.float64, ("lat",)) - wlon = outfile.createVariable("lon", np.float64, ("lon",)) - wmask = outfile.createVariable("landmask", np.int32, ("lat", "lon")) - warea = outfile.createVariable("area", np.float64, ("lat", "lon")) - wfrac = outfile.createVariable("landfrac", np.float64, ("lat", "lon")) - wtime.units = ssp_time_units - wlon.units = "degrees_east" - wlat.units = "degrees_north" - warea.units = "km2" - wfrac.units = "unitless" - wmask.units = "unitless" + # Create dimensions + outfile = create_fill_latlon(outfile, lat, "lat") + outfile = create_fill_latlon(outfile, lon, "lon") + outfile = create_fill_time(outfile, ssp_time, None, ssp_time_units=ssp_time_units, ssp_time_longname=ssp_time_longname, adj_time=True) - # wtime.long_name = 'Months since January '+str(fut_yrstart) - wtime.long_name = ssp_time_longname - wlon.long_name = "Longitude" - wlat.long_name = "Latitude" - warea.long_name = "Grid cell area" - wfrac.long_name = "Grid cell land fraction" - wmask.long_name = "Grid cell land mask" - wlon.mode = "time-invariant" - wlat.mode = "time-invariant" - warea.mode = "time-invariant" - wfrac.mode = "time-invariant" - wmask.mode = "time-invariant" - - wtime.calendar = "noleap" - - # write to file -------------------------------------------- - # wtime_offset = 0 - # adjust time to middle of month - # wtime_offset = -15 - wtime_offset = 15 - ssp_time[0] - wtime[:] = ssp_time + wtime_offset - wtime.calendar = "noleap" - wlon[:] = lon - wlat[:] = lat - wmask[:, :] = landmask - wfrac[:, :] = landfrac - warea[:, :] = area + # Create and fill ancillary variables + outfile = create_fill_ancillary_vars(outfile, landfrac, landmask, area) # -- End if on open file - if field_out[f] == "sfcWind": - wvar = outfile.createVariable( - field_out_wind[0], - np.float64, - ("time", "lat", "lon"), - fill_value=np.float64(1.0e36), - ) - else: - wvar = outfile.createVariable( - field_out[f], - np.float64, - ("time", "lat", "lon"), - fill_value=np.float64(1.0e36), - ) - wvar.units = units[f] - wvar.mode = "time-dependent" - - # write to file -------------------------------------------- - if field_out[f] == "sfcWind": - wvar.long_name = str(long_name) + " U component " + anomsf[f] - else: - wvar.long_name = str(long_name) + " " + anomsf[f] - - if field_out[f] == "sfcWind": - wvar[:, :, :] = anom_fld / np.sqrt(2) - else: - wvar[:, :, :] = anom_fld - - # List of source files - wvar.historical_source_files = "".join(histfiles).replace(hdir, "") - wvar.scenario_source_files = "".join(sspfiles).replace(fdir, "") - - # create second wind field for V component - if field_out[f] == "sfcWind": - command = 'date "+%y%m%d"' - x2 = subprocess.Popen(command, stdout=subprocess.PIPE, shell="True") - x = x2.communicate() - timetag = x[0].decode("utf-8").strip() - - wvar = outfile.createVariable( - field_out_wind[1], - np.float64, - ("time", "lat", "lon"), - fill_value=np.float64(1.0e36), - ) - wvar.units = units[f] - wvar.cell_methods = "time: mean" - wvar.long_name = str(long_name) + " V component " + anomsf[f] - - # write to file -------------------------------------------- - wvar[:, :, :] = anom_fld / np.sqrt(2) - - # List of source files - wvar.historical_source_files = "".join(histfiles).replace(hdir, "") - wvar.scenario_source_files = "".join(sspfiles).replace(fdir, "") - - # -- end if statement for write for V field -------- + outfile = create_fill_forcing(outfile, field_out, units, anomsf, field_out_wind, f, hdir, fdir, histfiles, sspfiles, long_name, anom_fld) # -- End Loop over forcing fields ------------------------------------ outfile.close() -print("\n\nSuccessfully made anomoly forcing datasets\n") +print("\n\nSuccessfully made anomaly forcing datasets\n") diff --git a/tools/site_and_regional/modify_singlept_site_neon b/tools/site_and_regional/modify_singlept_site_neon new file mode 100755 index 0000000000..1b790a74ca --- /dev/null +++ b/tools/site_and_regional/modify_singlept_site_neon @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 +""" +This is a just top-level skeleton script that calls +modify_singlept_site_neon.py. +The original code (modify_singlept_site_neon.py) is located under +python/ctsm/site_and_regional folder. + +For full instructions on how to run the code and different options, +please check python/ctsm/site_and_regional/modify_singlept_site_neon.py file. + +This script is for modifying surface dataset at neon sites +using data available from the neon server. + +After creating a single point surface data file from a global +surface data file using subset_data.py, use this script to +overwrite some fields with site-specific data for neon sites. + +This script will do the following: +- Download neon data for the specified site if it does not exist + in the specified directory : (i.e. ../../../neon_surf_files). +- Modify surface dataset with downloaded data. + +---------------------------------------------------------------- +To see all available options for modifying surface datasets at +tower sites: + ./modify_singlept_site_neon --help +---------------------------------------------------------------- +Instructions for running using conda python environments: +../../py_env_create +conda activate ctsm_pylib +""" + +import os +import sys + +# -- add python/ctsm to path +_CTSM_PYTHON = os.path.join( + os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python" +) +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm.site_and_regional.modify_singlept_site_neon import main + +if __name__ == "__main__": + main() diff --git a/tools/site_and_regional/neon_surf_wrapper b/tools/site_and_regional/neon_surf_wrapper new file mode 100755 index 0000000000..306d38a774 --- /dev/null +++ b/tools/site_and_regional/neon_surf_wrapper @@ -0,0 +1,37 @@ +#!/usr/bin/env python3 +""" +This is a just top-level skeleton script that calls +neon_surf_wrapper.py. +The original code (neon_surf_wrapper.py) is located under +python/ctsm/site_and_regional folder. + +For full instructions on how to run the code and different options, +please check python/ctsm/site_and_regional/neon_surf_wrapper.py file. + +This script is a simple wrapper for neon sites that performs the +following: + 1) For neon sites, subset surface dataset from global dataset + (i.e. ./subset_data.py ) + 2) Download neon and update the created surface dataset + based on the downloaded neon data. + (i.e. modify_singlept_site_neon.py) + +---------------------------------------------------------------- +Instructions for running using conda python environments: +../../py_env_create +conda activate ctsm_pylib +""" + +import os +import sys + +# -- add python/ctsm to path +_CTSM_PYTHON = os.path.join( + os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python" +) +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm.site_and_regional.neon_surf_wrapper import main + +if __name__ == "__main__": + main() diff --git a/tools/site_and_regional/neon_surf_wrapper.py b/tools/site_and_regional/neon_surf_wrapper.py deleted file mode 100755 index 3271c72f08..0000000000 --- a/tools/site_and_regional/neon_surf_wrapper.py +++ /dev/null @@ -1,157 +0,0 @@ -#! /usr/bin/env python3 -""" -|------------------------------------------------------------------| -|--------------------- Instructions -----------------------------| -|------------------------------------------------------------------| -This script is a simple wrapper for neon sites that performs the -following: - 1) For neon sites, subset surface dataset from global dataset - (i.e. ./subset_data.py ) - 2) Download neon and update the created surface dataset - based on the downloaded neon data. - (i.e. modify_singlept_site_neon.py) - -Instructions for running using conda python environments: - -../../py_env_create -conda activate ctsm_py - -""" -# TODO -# Automatic downloading of missing files if they are missing -#-[ ] Download neon sites and dom pft file -#-[ ] Make sure verbose works for printing out commands running - -# Import libraries -from __future__ import print_function - -import os -import sys -import tqdm -import logging -import argparse -import subprocess - -import pandas as pd - - - - -def get_parser(): - """ - Get parser object for this script. - """ - parser = argparse.ArgumentParser(description=__doc__, - formatter_class=argparse.RawDescriptionHelpFormatter) - - parser.print_usage = parser.print_help - - parser.add_argument('-v','--verbose', - help='Verbose mode will print more information. ', - action="store_true", - dest="verbose", - default=False) - - parser.add_argument('--16pft', - help='Create and/or modify 16-PFT surface datasets (e.g. for a FATES run) ', - action="store_true", - dest="pft_16", - default=False) - - parser.add_argument('-m', '--mixed', - help='Do not overwrite surface dataset to be just one dominant PFT at 100%', - action="store_true", - dest="mixed", - default=False) - - - return parser - - -def execute(command): - """ - Function for running a command on shell. - Args: - command (str): - command that we want to run. - Raises: - Error with the return code from shell. - """ - print ('\n',' >> ',*command,'\n') - - try: - subprocess.check_call(command, stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT) - - except subprocess.CalledProcessError as e: - #raise RuntimeError("command '{}' return with error (code {}): {}".format(e.cmd, e.returncode, e.output)) - #print (e.ouput) - print (e) - - - - - - -def main(): - - args = get_parser().parse_args() - - if args.verbose: - logging.basicConfig(level=logging.DEBUG) - - - neon_sites = pd.read_csv('neon_sites_dompft.csv') - - - for i, row in tqdm.tqdm(neon_sites.iterrows()): - lat = row['Lat'] - lon = row['Lon'] - site = row['Site'] - pft = row['pft'] - clmsite = "1x1_NEON_"+site - print ("Now processing site :", site) - - if args.mixed and args.pft_16: - # use surface dataset with 16 pfts, and don't overwrite with 100% 1 dominant PFT - # don't set crop flag - # don't set a dominant pft - subset_command = ['./subset_data','point','--lat',str(lat),'--lon',str(lon), - '--site',clmsite, '--create-surface','--uniform-snowpack', - '--cap-saturation','--verbose','--overwrite'] - modify_command = ['./modify_singlept_site_neon.py', '--neon_site', site, '--surf_dir', - 'subset_data_single_point', '--16pft'] - elif args.pft_16: - # use surface dataset with 16 pfts, but overwrite to 100% 1 dominant PFT - # don't set crop flag - # set dominant pft - subset_command = ['./subset_data','point','--lat',str(lat),'--lon',str(lon), - '--site',clmsite,'--dompft',str(pft),'--create-surface', - '--uniform-snowpack','--cap-saturation','--verbose','--overwrite'] - modify_command = ['./modify_singlept_site_neon.py', '--neon_site', site, '--surf_dir', - 'subset_data_single_point', '--16pft'] - elif args.mixed: - # use surface dataset with 78 pfts, and don't overwrite with 100% 1 dominant PFT - # NOTE: FATES will currently not run with a 78-PFT surface dataset - # set crop flag - # don't set dominant pft - subset_command = ['./subset_data','point','--lat',str(lat),'--lon',str(lon), - '--site',clmsite,'--crop','--create-surface', - '--uniform-snowpack','--cap-saturation','--verbose','--overwrite'] - modify_command = ['./modify_singlept_site_neon.py', '--neon_site', site, '--surf_dir', - 'subset_data_single_point'] - else: - # use surface dataset with 78 pfts, and overwrite to 100% 1 dominant PFT - # NOTE: FATES will currently not run with a 78-PFT surface dataset - # set crop flag - # set dominant pft - subset_command = ['./subset_data', 'point', '--lat', str(lat), '--lon', str(lon), - '--site', clmsite,'--crop', '--dompft', str(pft), '--create-surface', - '--uniform-snowpack', '--cap-saturation', '--verbose', '--overwrite'] - modify_command = ['./modify_singlept_site_neon.py', '--neon_site', site, '--surf_dir', - 'subset_data_single_point'] - execute(subset_command) - execute(modify_command) - -if __name__ == "__main__": - main() - diff --git a/tools/site_and_regional/run_neon b/tools/site_and_regional/run_neon new file mode 100755 index 0000000000..ad930f50e3 --- /dev/null +++ b/tools/site_and_regional/run_neon @@ -0,0 +1,47 @@ +#!/usr/bin/env python3 +""" +This is a just top-level skeleton script that calls +run_neon.py. +The original code (run_neon.py) is located under +python/ctsm/site_and_regional folder. + +For full instructions on how to run the code and different options, +please check python/ctsm/site_and_regional/run_neon.py file. + +This script first creates and builds a generic base case. +Next, it will clone the base_case for different neon sites and run +types to reduce the need to build ctsm everytime. + +This script will do the following: + 1) Create a generic base case for cloning. + 2) Make the case for the specific neon site(s). + 3) Make changes to the case, for: + a. AD spinup + b. post-AD spinup + c. transient + #--------------- + d. SASU or Matrix spinup + 4) Build and submit the case. + +---------------------------------------------------------------- +To see all available options for running tower sites: + ./run_neon --help +---------------------------------------------------------------- +Instructions for running using conda python environments: +../../py_env_create +conda activate ctsm_pylib +""" + +import os +import sys + +# -- add python/ctsm to path +_CTSM_PYTHON = os.path.join( + os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python" +) +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm.site_and_regional.run_neon import main + +if __name__ == "__main__": + main(__doc__)