diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 256de592c6..69b3a8c3dc 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1793,6 +1793,11 @@ sub process_namelist_inline_logic { # namelist group: clm_initinterp_inparm # ######################################### setup_logic_initinterp($opts, $nl_flags, $definition, $defaults, $nl); + + ############################### + # namelist group: exice_streams # + ############################### + setup_logic_exice($opts, $nl_flags, $definition, $defaults, $nl); } #------------------------------------------------------------------------------- @@ -2197,6 +2202,7 @@ sub setup_logic_soilstate { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'organic_frac_squared' ); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_bedrock', 'use_fates'=>$nl_flags->{'use_fates'}, 'vichydro'=>$nl_flags->{'vichydro'} ); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_excess_ice'); # excess ice flag should be read before stream vars my $var1 = "soil_layerstruct_predefined"; my $var2 = "soil_layerstruct_userdefined"; @@ -4306,6 +4312,54 @@ sub setup_logic_fates { } } +#------------------------------------------------------------------------------- +sub setup_logic_exice { + # + # excess ice streams + # + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + my $use_exice = $nl->get_value( 'use_excess_ice' ); + my $use_exice_streams = $nl->get_value( 'use_excess_ice_streams' ); + # IF excess ice streams is on + if (defined($use_exice_streams) && value_is_true($use_exice_streams)) { + # Can only be true if excess ice is also on, otherwise fail + if (defined($use_exice) && not value_is_true($use_exice)) { + $log->fatal_error("use_excess_ice_streams can NOT be TRUE when use_excess_ice is FALSE" ); + } + # Otherwise if ice streams are off + } else { + my @list = ( "stream_meshfile_exice", "stream_fldfilename_exice" ); + # fail is excess ice streams files are set + foreach my $var ( @list ) { + if ( defined($nl->get_value($var)) ) { + $log->fatal_error("$var should NOT be set when use_excess_ice_streams=FALSE" ); + } + } + # mapalgo can only be none, if excess ice streams are off + my $map_algo = $nl->get_value("stream_mapalgo_exice"); + if ( defined($map_algo) && ($map_algo ne "none") ) { + $log->fatal_error("stream_mapalgo_exice can ONLY be none when use_excess_ice_streams=FALSE" ); + } + } + # If excess ice is on + if (defined($use_exice) && value_is_true($use_exice)) { + # IF nuopc driver and excess ice streams are on get the stream defaults + if (defined($use_exice_streams) && value_is_true($use_exice_streams)) { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_exice'); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_mapalgo_exice'); + # If excess ice streams on, but NOT the NUOPC driver fail + if ( not $opts->{'driver'} eq "nuopc" ) { + $log->fatal_error("nuopc driver is required when use_excess_ice_streams is set to true" ); + # NUOPC driver needs a mesh file + } else { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_exice'); + } + } + } + + +} # end exice streams + #------------------------------------------------------------------------------- sub setup_logic_misc { @@ -4377,6 +4431,7 @@ sub write_output_files { push @groups, "nitrif_inparm"; push @groups, "lifire_inparm"; push @groups, "ch4finundated"; + push @groups, "exice_streams"; push @groups, "soilbgc_decomp"; push @groups, "clm_canopy_inparm"; if (remove_leading_and_trailing_quotes($nl->get_value('snow_cover_fraction_method')) eq 'SwensonLawrence2012') { diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 3cf9a3ebc0..1d895524e6 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2713,4 +2713,14 @@ use_crop=".true.">lnd/clm2/surfdata_map/ctsm5.1.dev052/landuse.timeseries_mpasa1 general + + + + +.false. + +lnd/clm2/paramdata/exice_init_0.125x0.125_c20220516.nc +lnd/clm2/paramdata/exice_init_0.125x0.125_ESMFmesh_cdf5_c20220802.nc +bilinear + diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index dc693b50ba..45a32cb224 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2848,4 +2848,38 @@ use case.) + + + + +If TRUE turn on the excess ice physics, (Lee et al., 2014; Cai et al., 2020) + + + +If TRUE and use_excess_ice is TRUE, use the excess ice stream to determine the initial values of the excess ice field +if FALSE and use_excess_ice is TRUE, expect excess ice to come from the initial conditions or restart file +Expect to be FALSE is use_excess_ice is FALSE + + + +Filename of input stream data for excess ice data + + + +mesh filename of input stream data for excess ice + + + +Mapping method from excess ice input stream data to the model resolution + bilinear = bilinear interpolation + nn = nearest neighbor + none = no interpolation + + + diff --git a/bld/unit_testers/build-namelist_test.pl b/bld/unit_testers/build-namelist_test.pl index ed6328c166..82a178bf03 100755 --- a/bld/unit_testers/build-namelist_test.pl +++ b/bld/unit_testers/build-namelist_test.pl @@ -163,9 +163,9 @@ sub cat_and_create_namelistinfile { # # Figure out number of tests that will run # -my $ntests = 1975; +my $ntests = 1992; if ( defined($opts{'compare'}) ) { - $ntests += 1344; + $ntests += 1353; } plan( tests=>$ntests ); @@ -323,6 +323,8 @@ sub cat_and_create_namelistinfile { "-res 0.9x1.25 -structure fast", "-res 0.9x1.25 -namelist '&a irrigate=.true./'", "-res 0.9x1.25 -verbose", "-res 0.9x1.25 -ssp_rcp SSP1-2.6", "-res 0.9x1.25 -test", "-res 0.9x1.25 -sim_year 1850", "-res 0.9x1.25 -namelist '&a use_lai_streams=.true.,use_soil_moisture_streams=.true./'", + "-res 0.9x1.25 -namelist '&a use_excess_ice=.true. use_excess_ice_streams=.true./'", + "-res 0.9x1.25 -namelist '&a use_excess_ice=.true. use_excess_ice_streams=.false./'", "-res 0.9x1.25 -use_case 1850_control", "-res 1x1pt_US-UMB -clm_usr_name 1x1pt_US-UMB -namelist '&a fsurdat=\"/dev/null\"/'", "-res 1x1_brazil", @@ -337,6 +339,10 @@ sub cat_and_create_namelistinfile { my $base_options = "-envxml_dir . -driver $driver"; if ( $driver eq "mct" ) { $base_options = "$base_options -lnd_frac $DOMFILE"; + # Skip the MCT test for excess ice streams + if ( $options =~ /use_excess_ice_streams=.true./ ) { + next; + } } else { $base_options = "$base_options -namelist '&a force_send_to_atm = .false./'"; } @@ -524,8 +530,33 @@ sub cat_and_create_namelistinfile { GLC_TWO_WAY_COUPLING=>"FALSE", phys=>"clm4_5", }, - "soilm_stream wo use" =>{ options=>"-res 0.9x1.25 -envxml_dir .", - namelst=>"use_soil_moisture_streams = .false.,stream_fldfilename_soilm='missing_file'", + "soilm_stream off w file" =>{ options=>"-res 0.9x1.25 -envxml_dir .", + namelst=>"use_soil_moisture_streams = .false.,stream_fldfilename_soilm='file_provided_when_off'", + GLC_TWO_WAY_COUPLING=>"FALSE", + phys=>"clm5_0", + }, + "exice_stream off w file" =>{ options=>"-res 0.9x1.25 -envxml_dir .", + namelst=>"use_excess_ice=.true., use_excess_ice_streams = .false.,stream_fldfilename_exice='file_provided_when_off'", + GLC_TWO_WAY_COUPLING=>"FALSE", + phys=>"clm5_0", + }, + "exice_stream off w mesh" =>{ options=>"-res 0.9x1.25 -envxml_dir .", + namelst=>"use_excess_ice=.true., use_excess_ice_streams = .false.,stream_meshfile_exice='file_provided_when_off'", + GLC_TWO_WAY_COUPLING=>"FALSE", + phys=>"clm5_0", + }, + "exice off, but stream on" =>{ options=>"-res 0.9x1.25 -envxml_dir .", + namelst=>"use_excess_ice=.false., use_excess_ice_streams = .true.,stream_fldfilename_exice='file_provided', stream_meshfile_exice='file_provided'", + GLC_TWO_WAY_COUPLING=>"FALSE", + phys=>"clm5_0", + }, + "exice stream off, but setmap"=>{ options=>"-res 0.9x1.25 -envxml_dir .", + namelst=>"use_excess_ice=.true., use_excess_ice_streams = .false.,stream_mapalgo_exice='bilinear'", + GLC_TWO_WAY_COUPLING=>"FALSE", + phys=>"clm5_0", + }, + "exice stream on, but mct" =>{ options=>"--res 0.9x1.25 --envxml_dir . --driver mct --lnd_frac $DOMFILE ", + namelst=>"use_excess_ice=.true., use_excess_ice_streams=.true.", GLC_TWO_WAY_COUPLING=>"FALSE", phys=>"clm5_0", }, diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index 24f7062798..76af28337f 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -29,7 +29,6 @@ - FAIL diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 31c44c307d..aa5dd98697 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -158,6 +158,24 @@ + + + + + + + + + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/clm/ExcessIceStartup_output_sp_exice/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStartup_output_sp_exice/include_user_mods new file mode 100644 index 0000000000..6d8de3732a --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStartup_output_sp_exice/include_user_mods @@ -0,0 +1,2 @@ +../monthly +../../../../usermods_dirs/output_sp_exice diff --git a/cime_config/testdefs/testmods_dirs/clm/ExcessIceStartup_output_sp_exice/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStartup_output_sp_exice/user_nl_clm new file mode 100644 index 0000000000..e479af9449 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStartup_output_sp_exice/user_nl_clm @@ -0,0 +1,3 @@ + use_excess_ice = .true. + finidat = '$DIN_LOC_ROOT/lnd/clm2/initdata_map/clmi.I1850Clm50Sp.0003-01-01.0.9x1.25_gx1v7_exice_simyr1850_c221205.nc' + use_init_interp = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/ExcessIceStreams/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStreams/include_user_mods new file mode 100644 index 0000000000..fe0e18cf88 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStreams/include_user_mods @@ -0,0 +1 @@ +../default diff --git a/cime_config/testdefs/testmods_dirs/clm/ExcessIceStreams/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStreams/user_nl_clm new file mode 100644 index 0000000000..f61ca32a79 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/ExcessIceStreams/user_nl_clm @@ -0,0 +1,2 @@ + use_excess_ice = .true. + use_excess_ice_streams = .true. diff --git a/cime_config/usermods_dirs/output_sp_exice/include_user_mods b/cime_config/usermods_dirs/output_sp_exice/include_user_mods new file mode 100644 index 0000000000..786fa907a9 --- /dev/null +++ b/cime_config/usermods_dirs/output_sp_exice/include_user_mods @@ -0,0 +1,2 @@ +../_includes/output_base +../output_sp diff --git a/cime_config/usermods_dirs/output_sp_exice/user_nl_clm b/cime_config/usermods_dirs/output_sp_exice/user_nl_clm new file mode 100644 index 0000000000..48e680df67 --- /dev/null +++ b/cime_config/usermods_dirs/output_sp_exice/user_nl_clm @@ -0,0 +1,8 @@ +!---------------------------------------------------------------------------------- +! Settings from output_sp_exice +!---------------------------------------------------------------------------------- + +! h3 stream (yearly average, gridcell-level) +! Eyr +hist_fincl4 += 'EXCESS_ICE' + diff --git a/doc/ChangeLog b/doc/ChangeLog index ec7f8303a1..7fa07e58cf 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,123 @@ =============================================================== +Tag name: ctsm5.1.dev132 +Originator(s): mvdebolskiy (NORCE, Bergen, Norway), slevis (Samuel Levis,UCAR/TSS,303-665-1310) +Date: Fri Aug 4 17:52:45 MDT 2023 +One-line Summary: Add parameterization to allow excess ice in soil and subsidence + +Purpose and description of changes +---------------------------------- + +As described in PR #1787: + +Parameterization for excess ice described in Lee et al. (2014): +http://dx.doi.org/10.1088/1748-9326/9/12/124006 + +This code is a modified version of code provided by Lei Cai: +https://github.com/lca041/ctsm/tree/clm5.0.dev92_exice + +Works only for the nuopc driver. + +Files changed: +bld/CLMBuildNamelist.pm, bld/namelist_files/namelist_defaults_ctsm.xml, bld/namelist_files/namelist_definitionss_ctsm.xml -- added namelist options; +src/main/clm_varctl.F90, src/main/controlMod.F90 -- added option to switch excess ice physics and read namelist option; +src/biogeophys/WaterStateType.F90 -- added prognostic excess ice variable and a history field; +src/biogeophys/WaterStateBulkType.F90, src/main/clm_instMod.F90, -- added arguments to soubrutiens for proper initialization; +src/biogeophys/TemperatureType.F90 -- initial soil temperature set to 268.15 K at the cold start (possibly redundant because #1460 is closed) +src/biogeophys/WaterDiagnosticBulkType.F90 -- added two diagnostic excess ice variables and two history fields; +src/biogeophys/SoilTemperatureMod.F90 -- added main excess ice calculations; +src/biogeophys/TotalWaterAndHeatMod.F90 -- added excess ice to total water for balance checks; +src/biogeophys/SoilHydrologyMod.F90 -- added excess ice for ice fraction calculation; + +New files: +src/cpl/share_esmf/ExcessIceStreamType.F90 -- routines to read dataset with initial excess ice concentration + + +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 #): + Fixes #1229 -- excess ice + + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): + Excess ice can EITHER be turned on by using the stream file, OR a restart file that has excess ice on it. + Since, excess ice is expected to melt as time goes on, the use of a restart file is preferred. + But, use of a restart file requires a simulation that was spun up to that point. + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): + New namelist options added: + - use_excess_ice (logical, in clm_inparm) default = .false.; turns on excess ice physics + - stream_meshfile_exice, stream_fldfilename_exice, stream_mapalgo_exice (char, in exice_streams) + meshfile, stream file, spatial interpolation algorithm for initial values of excess ice + defaults - lnd/clm2/paramdata/exice_init_0.125x0.125_ESMFmesh_c20220516.nc, + lnd/clm2/paramdata/exice_init_0.125x0.125_c20220516.nc + and bilinear + Dataset interpolated to 0.125x0.125 degrees grid from Brown et al. (1997) can be found here: + https://drive.google.com/file/d/1mA457Oa52zG_MtLGB7KHuUYQvsS2-P5o/view?usp=sharing + Dataset used only in cold start or hybrid runs (or starting with finidat) that do not have excess ice + + +Notes of particular relevance for developers: +--------------------------------------------- +Changes to tests or testing: + New tests in place for this new code + +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 - PASS + + tools-tests (test/tools) (if tools have been changed): + + cheyenne - + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + cheyenne ---- OK + izumi ------- OK + + 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: No + + +Other details +------------- +Pull Requests that document the changes (include PR ids): + https://github.com/ESCOMP/ctsm/pull/1787 + +=============================================================== +=============================================================== Tag name: ctsm5.1.dev131 Originator(s): samrabin (Sam Rabin,UCAR/TSS) Date: Thu Jul 27 14:24:07 MDT 2023 diff --git a/doc/ChangeSum b/doc/ChangeSum index 2d1812cd13..1702a919d0 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.1.dev132 slevis 08/04/2023 Add parameterization to allow excess ice in soil and subsidence ctsm5.1.dev131 samrabin 07/27/2023 Enable prescribed crop calendars ctsm5.1.dev130 glemieux 07/09/2023 FATES parameter file and test definition update ctsm5.1.dev129 erik 06/22/2023 NEON fixes for TOOL and user-mods, add SP for NEON, some history file updates, black refactor for buildlib/buildnml diff --git a/doc/source/tech_note/References/CLM50_Tech_Note_References.rst b/doc/source/tech_note/References/CLM50_Tech_Note_References.rst index 4bd7a2cf9b..46b835a76d 100644 --- a/doc/source/tech_note/References/CLM50_Tech_Note_References.rst +++ b/doc/source/tech_note/References/CLM50_Tech_Note_References.rst @@ -223,6 +223,13 @@ Geosci. Model Dev., 7, 2193-2222, doi:10.5194/gmd-7-2193-2014. Botta, A et al., 2000. A global prognostic scheme of leaf onset using satellite data. Global Change Biology 6.7, pp. 709-725. +.. _Brownetal1997: + +Brown J., Ferrians O. J. Jr, Heginbottom J. A. and Melnikov E. S. 1997. + Circum-Arctic Map of Permafrost and Ground-Ice Conditions +(Boulder, CO: National Snow and Ice Data Center) version 2, +DOI: 10.3133/cp45 + .. _Brun1989: Brun, E. 1989. Investigation of wet-snow metamorphism in respect of @@ -1085,6 +1092,12 @@ Climate 25:3071-3095. DOI:10.1175/JCLI-D-11-00256.1. Lehner, B. and Döll, P., 2004. Development and validation of a global database of lakes, reservoirs and wetlands, J. Hydrol., 296, 1–22. +.. _Leeetal2014: + +Lee, H., Swenson, S.C., Slater A.G. and Lawrence D.M., 2014. Effects +of excess ground ice on projections of permafrost in a warming climate. +Environmental Research Letters 9:12 124006. DOI: 10.1088/1748-9326/9/12/124006 + .. _Lehneretal2008: Lehner, B., Verdin, K. and Jarvis, A., 2008. New global hydrograhy diff --git a/doc/source/tech_note/Soil_Snow_Temperatures/CLM50_Tech_Note_Soil_Snow_Temperatures.rst b/doc/source/tech_note/Soil_Snow_Temperatures/CLM50_Tech_Note_Soil_Snow_Temperatures.rst index 72040a1514..ab10547339 100644 --- a/doc/source/tech_note/Soil_Snow_Temperatures/CLM50_Tech_Note_Soil_Snow_Temperatures.rst +++ b/doc/source/tech_note/Soil_Snow_Temperatures/CLM50_Tech_Note_Soil_Snow_Temperatures.rst @@ -979,3 +979,50 @@ the top layer is a blend of ice and soil heat capacity c_{1} =c_{1}^{*} +\frac{C_{ice} W_{sno} }{\Delta z_{1} } where :math:`c_{1}^{*}` is calculated from :eq:`6.89` or :eq:`6.92`. + + +.. _Excess Ground Ice: + +Excess Ground Ice +------------------------------------ + +An optional parameterization of excess ground ice melt and respective subsidence based on (:ref:`Lee et al., (2014) `). +Initial excess ground ice concentrations for soil columns are derived from (:ref:`Brown et al., (1997) `). +When the excess ground ice is present in the soil column, soil depth for a given layer (:math:`z_{i}`) +is adjusted by the amount of excess ice in the column: + +.. math:: + :label: 6.94 + + z_{i}^{'}=\Sigma_{j=1}^{i} \ z_{j}^{'}+\frac{w_{exice,\, j}}{\rho_{ice} } + +where :math:`w_{exice,\,j}` is excess ground ice amount (kg m :sup:`-2`) in layer :math:`j` +and :math:`\rho_{ice}` is the density of ice (kg m :sup:`-3`). +After adjustment of layer depths have been made, all of the soil temperature equations (from :eq:`6.80` to :eq:`6.89`) +are calculted based on the adjusted depths. Thermal properties are additionally adjusted (:eq:`6.8` and :eq:`6.8`) in the following way: + +.. math:: + :label: 6.95 + + \begin{array}{lr} + \theta_{sat}^{'} =\frac{\theta _{liq} }{\theta _{liq} +\theta _{ice} +\theta_{exice}}{\theta_{sat}} \\ + \lambda _{sat}^{'} =\lambda _{s}^{1-\theta _{sat}^{'} } \lambda _{liq}^{\frac{\theta _{liq} }{\theta _{liq} +\theta _{ice} +\theta_{exice}} \theta _{sat}^{'} } \lambda _{ice}^{\theta _{sat}^{'} \left(1-\frac{\theta _{liq} }{\theta _{liq} +\theta _{ice} +\theta_{exice}} \right)} \\ + c_{i}^{'} =c_{s,\, i} \left(1-\theta _{sat,\, i}^{'} \right)+\frac{w_{ice,\, i} +w_{exice,\,j}}{\Delta z_{i}^{'} } C_{ice} +\frac{w_{liq,\, i} }{\Delta z_{i}^{'} } C_{liq} + \end{array} + +Soil subsidence at the timestep :math:`n+1` (:math:`z_{exice}^{n+1}`, m) is then calculated as: + +.. math:: + :label: 6.96 + + z_{exice}^{n+1}=\Sigma_{i=1}^{N_{levgrnd}} \ z_{j}^{',\ ,n+1}-z_{j}^{',\ ,n } + +With regards to hydraulic counductivity, excess ground ice is treated the same way normal soil +ice is treated in :numref:`Frozen Soils and Perched Water Table`. +When a soil layer thaws, excess ground ice is only allowed +to melt when no normals soil ice is present in the layer. +When a soil layer refreezes, liquid soil water can only turn into normal soil ice, thus, no new of excess ice can be created but only melted. +The excess liquid soil moisture from excess ice melt is distributed within the soil column according +to :numref:`Lateral Sub-surface Runoff`. + + diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index 5cdffba4ff..4bc6a784de 100644 --- a/src/biogeophys/SoilHydrologyMod.F90 +++ b/src/biogeophys/SoilHydrologyMod.F90 @@ -179,6 +179,7 @@ subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, & integer :: j, fc, c real(r8) :: vol_ice(bounds%begc:bounds%endc,1:nlevsoi) !partial volume of ice lens in layer real(r8) :: icefrac_orig ! original formulation for icefrac + real(r8) :: dz_ext(bounds%begc:bounds%endc,1:nlevsoi) character(len=*), parameter :: subname = 'SetSoilWaterFractions' !----------------------------------------------------------------------- @@ -189,8 +190,9 @@ subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, & watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) eff_porosity => soilstate_inst%eff_porosity_col , & ! Output: [real(r8) (:,:) ] effective porosity = porosity - vol_ice - h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) - h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice water (kg/m2) + h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice water (kg/m2) + excess_ice => waterstatebulk_inst%excess_ice_col , & ! Input: [real(r8) (:,:) ] excess ice (kg/m2) origflag => soilhydrology_inst%origflag , & ! Input: logical icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] @@ -203,7 +205,8 @@ subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, & ! Porosity of soil, partial volume of ice and liquid, fraction of ice in each layer, ! fractional impermeability - vol_ice(c,j) = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) + dz_ext(c,j) = dz(c,j) + excess_ice(c,j)/denice ! extended layer thickness, should be good for all the columns + vol_ice(c,j) = min(watsat(c,j), (h2osoi_ice(c,j) + excess_ice(c,j))/(dz_ext(c,j)*denice)) eff_porosity(c,j) = max(0.01_r8,watsat(c,j)-vol_ice(c,j)) icefrac(c,j) = min(1._r8,vol_ice(c,j)/watsat(c,j)) diff --git a/src/biogeophys/SoilTemperatureMod.F90 b/src/biogeophys/SoilTemperatureMod.F90 index 513413e8a9..50c3cf15e5 100644 --- a/src/biogeophys/SoilTemperatureMod.F90 +++ b/src/biogeophys/SoilTemperatureMod.F90 @@ -114,8 +114,8 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter ! !USES: use clm_time_manager , only : get_step_size_real use clm_varpar , only : nlevsno, nlevgrnd, nlevurb, nlevmaxurbgrnd - use clm_varctl , only : iulog - use clm_varcon , only : cnfac, cpice, cpliq, denh2o + use clm_varctl , only : iulog, use_excess_ice + use clm_varcon , only : cnfac, cpice, cpliq, denh2o, denice use landunit_varcon , only : istsoil, istcrop use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv use BandDiagonalMod , only : BandDiagonal @@ -171,6 +171,10 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter real(r8) :: hs_top_snow(bounds%begc:bounds%endc) ! heat flux on top snow layer [W/m2] real(r8) :: hs_h2osfc(bounds%begc:bounds%endc) ! heat flux on standing water [W/m2] integer :: jbot(bounds%begc:bounds%endc) ! bottom level at each column + real(r8) :: dz_0(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) ! original layer thickness [m] + real(r8) :: z_0(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) ! original layer depth [m] + real(r8) :: zi_0(bounds%begc:bounds%endc,-nlevsno+0:nlevmaxurbgrnd) ! original layer interface level bellow layer "z" [m] + !----------------------------------------------------------------------- associate( & @@ -194,6 +198,7 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter frac_sno_eff => waterdiagnosticbulk_inst%frac_sno_eff_col , & ! Input: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1) snow_depth => waterdiagnosticbulk_inst%snow_depth_col , & ! Input: [real(r8) (:) ] snow height (m) h2osfc => waterstatebulk_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) + excess_ice => waterstatebulk_inst%excess_ice_col , & ! Input: [real(r8) (:,:) ] excess ice (kg/m2) (new) (1:nlevgrnd) frac_h2osfc => waterdiagnosticbulk_inst%frac_h2osfc_col , & ! Input: [real(r8) (:) ] fraction of ground covered by surface water (0 to 1) @@ -271,6 +276,29 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter endif end do + + !-------------------------------------------------------------- + ! Vertical coordinates adjustment for excess ice calculations + !-------------------------------------------------------------- + if ( use_excess_ice ) then + ! Save original soil depth to get put them back in et the end + dz_0(begc:endc,-nlevsno+1:nlevmaxurbgrnd) = dz(begc:endc,-nlevsno+1:nlevmaxurbgrnd) + zi_0(begc:endc,-nlevsno+0:nlevmaxurbgrnd) = zi(begc:endc,-nlevsno+0:nlevmaxurbgrnd) + z_0(begc:endc,-nlevsno+1:nlevmaxurbgrnd) = z(begc:endc,-nlevsno+1:nlevmaxurbgrnd) + ! Adjust column depth for excess ice thickness + do fc = 1, num_nolakec + c = filter_nolakec(fc) + l = col%landunit(c) + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + dz(c,1:nlevmaxurbgrnd) = dz(c,1:nlevmaxurbgrnd) + excess_ice(c,1:nlevmaxurbgrnd) / denice ! add extra layer thickness + do j = 1, nlevmaxurbgrnd ! if excess ice amount dropped to zero there will be no adjustment + zi(c,j) = zi(c,j) + sum(excess_ice(c,1:j)) / denice + z(c,j) = (zi(c,j-1) + zi(c,j)) * 0.5_r8 + end do + end if + end do + end if + !------------------------------------------------------ ! Compute ground surface and soil temperatures !------------------------------------------------------ @@ -488,6 +516,24 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter dhsdT(bounds%begc:bounds%endc), & soilstate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, energyflux_inst, temperature_inst) + !-------------------------------------------------------------- + ! Vertical coordinates adjustment for excess ice calculations + !-------------------------------------------------------------- + ! bringing back the soil depth to the original state + if (use_excess_ice) then + ! Adjust column depth for excess ice thickness + do fc = 1, num_nolakec + c = filter_nolakec(fc) + l = col%landunit(c) + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + dz(c,1:nlevmaxurbgrnd)=dz_0(c,1:nlevmaxurbgrnd) + zi(c,1:nlevmaxurbgrnd)=zi_0(c,1:nlevmaxurbgrnd) + z(c,1:nlevmaxurbgrnd)=z_0(c,1:nlevmaxurbgrnd) + end if + end do + end if + + if ( IsProgBuildTemp() )then call BuildingTemperature(bounds, num_urbanl, filter_urbanl, num_nolakec, filter_nolakec, & tk(bounds%begc:bounds%endc, :), urbanparams_inst, & @@ -628,6 +674,7 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter h2osno_no_layers => waterstatebulk_inst%h2osno_no_layers_col , & ! Input: [real(r8) (:) ] snow not resolved into layers (mm H2O) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) + excess_ice => waterstatebulk_inst%excess_ice_col , & ! Input: [real(r8) (:,:) ] excess ice lenses (kg/m2) (new) (1:nlevgrnd) bw => waterdiagnosticbulk_inst%bw_col , & ! Output: [real(r8) (:,:) ] partial density of water in the snow pack (ice + liquid) [kg/m3] tkmg => soilstate_inst%tkmg_col , & ! Input: [real(r8) (:,:) ] thermal conductivity, soil minerals [W/m-K] @@ -654,7 +701,7 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter col%itype(c) /= icol_roof .and. col%itype(c) /= icol_road_imperv) .or. & (col%itype(c) == icol_road_imperv .and. j > nlev_improad(l))) then - satw = (h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice)/(dz(c,j)*watsat(c,j)) + satw = (h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice +excess_ice(c,j)/denice)/(dz(c,j)*watsat(c,j)) satw = min(1._r8, satw) if (satw > .1e-6_r8) then if (t_soisno(c,j) >= tfrz) then ! Unfrozen soil @@ -663,7 +710,7 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter dke = satw end if fl = (h2osoi_liq(c,j)/(denh2o*dz(c,j))) / (h2osoi_liq(c,j)/(denh2o*dz(c,j)) + & - h2osoi_ice(c,j)/(denice*dz(c,j))) + h2osoi_ice(c,j)/(denice*dz(c,j))+excess_ice(c,j)/(denice*dz(c,j))) dksat = tkmg(c,j)*tkwat**(fl*watsat(c,j))*tkice**((1._r8-fl)*watsat(c,j)) thk(c,j) = dke*dksat + (1._r8-dke)*tkdry(c,j) else @@ -756,7 +803,8 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter .and. col%itype(c) /= icol_sunwall .and. col%itype(c) /= icol_shadewall .and. & col%itype(c) /= icol_roof .and. col%itype(c) /= icol_road_imperv) .or. & (col%itype(c) == icol_road_imperv .and. j > nlev_improad(l))) then - cv(c,j) = csol(c,j)*(1._r8-watsat(c,j))*dz(c,j) + (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq) + cv(c,j) = csol(c,j)*(1._r8-watsat(c,j))*dz(c,j) + (h2osoi_ice(c,j)*cpice + & + h2osoi_liq(c,j)*cpliq) + excess_ice(c,j)*cpice if (j > nbedrock(c)) cv(c,j) = csol_bedrock*dz(c,j) else if (lun%itype(l) == istwet) then cv(c,j) = (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq) @@ -1057,7 +1105,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & use clm_time_manager , only : get_step_size_real use clm_varpar , only : nlevsno, nlevgrnd, nlevurb, nlevmaxurbgrnd use clm_varctl , only : iulog - use clm_varcon , only : tfrz, hfus, grav + use clm_varcon , only : tfrz, hfus, grav, denice use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv use landunit_varcon , only : istsoil, istcrop, istice ! @@ -1081,13 +1129,17 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & real(r8) :: temp1 !temporary variables [kg/m2] real(r8) :: hm(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) !energy residual [W/m2] real(r8) :: xm(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) !melting or freezing within a time step [kg/m2] + real(r8) :: xm2(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) !additional melting or freezing within a time step [kg/m2] (needed for excess ice melt) real(r8) :: wmass0(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd)!initial mass of ice and liquid (kg/m2) real(r8) :: wice0 (bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd)!initial mass of ice (kg/m2) real(r8) :: wliq0 (bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd)!initial mass of liquid (kg/m2) real(r8) :: supercool(bounds%begc:bounds%endc,nlevmaxurbgrnd) !supercooled water in soil (kg/m2) - real(r8) :: propor !proportionality constant (-) + real(r8) :: propor !proportionality constant (-) real(r8) :: tinc(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) !t(n+1)-t(n) [K] - real(r8) :: smp !frozen water potential (mm) + real(r8) :: smp !frozen water potential (mm) + real(r8) :: wexice0(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd) !initial mass of excess_ice at the timestep (kg/m2) + + !----------------------------------------------------------------------- call t_startf( 'PhaseChangebeta' ) @@ -1106,9 +1158,11 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & frac_sno_eff => waterdiagnosticbulk_inst%frac_sno_eff_col , & ! Input: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1) frac_h2osfc => waterdiagnosticbulk_inst%frac_h2osfc_col , & ! Input: [real(r8) (:) ] fraction of ground covered by surface water (0 to 1) snow_depth => waterdiagnosticbulk_inst%snow_depth_col , & ! Input: [real(r8) (:) ] snow height (m) + exice_subs_col => waterdiagnosticbulk_inst%exice_subs_col , & ! Output: [real(r8) (:,:) ] per layer subsidence due to excess ice melt (mm/s) h2osno_no_layers => waterstatebulk_inst%h2osno_no_layers_col , & ! Output: [real(r8) (:) ] snow not resolved into layers (mm H2O) - h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) (new) - h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) (new) + h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) (new) + h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) (new) + excess_ice => waterstatebulk_inst%excess_ice_col , & ! Input: [real(r8) (:,:) ] excess ice (kg/m2) (new) (1:nlevgrnd) qflx_snow_drain => waterfluxbulk_inst%qflx_snow_drain_col , & ! Output: [real(r8) (:) ] drainage from snow pack qflx_snofrz_lyr => waterfluxbulk_inst%qflx_snofrz_lyr_col , & ! Output: [real(r8) (:,:) ] snow freezing rate (positive definite) (col,lyr) [kg m-2 s-1] @@ -1152,9 +1206,14 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & imelt(c,j) = 0 hm(c,j) = 0._r8 xm(c,j) = 0._r8 + xm2(c,j) = 0._r8 wice0(c,j) = h2osoi_ice(c,j) wliq0(c,j) = h2osoi_liq(c,j) - wmass0(c,j) = h2osoi_ice(c,j) + h2osoi_liq(c,j) + wexice0(c,j) = excess_ice(c,j) + wmass0(c,j) = h2osoi_ice(c,j) + h2osoi_liq(c,j) + wexice0(c,j) + if (j >= 1) then + exice_subs_col(c,j) = 0._r8 + endif endif ! end of snow layer if-block if (j <= 0) then @@ -1173,7 +1232,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & ! Melting identification ! If ice exists above melt point, melt some to liquid. - if (h2osoi_ice(c,j) > 0._r8 .AND. t_soisno(c,j) > tfrz) then + if (h2osoi_ice(c,j) > 0._r8 .and. t_soisno(c,j) > tfrz) then imelt(c,j) = 1 ! tinc(c,j) = t_soisno(c,j) - tfrz tinc(c,j) = tfrz - t_soisno(c,j) @@ -1211,6 +1270,13 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & t_soisno(c,j) = tfrz endif + ! melt excess ice after normal ice + if (excess_ice(c,j) > 0._r8 .AND. t_soisno(c,j) > tfrz) then + imelt(c,j) = 1 + tinc(c,j) = tfrz - t_soisno(c,j) + t_soisno(c,j) = tfrz + endif + ! from Zhao (1997) and Koren (1999) supercool(c,j) = 0.0_r8 if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop .or. col%itype(c) == icol_road_perv) then @@ -1325,23 +1391,31 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & endif heatr = 0._r8 - if (xm(c,j) > 0._r8) then + if (xm(c,j) > 0._r8) then !if there is excess heat to melt the ice h2osoi_ice(c,j) = max(0._r8, wice0(c,j)-xm(c,j)) heatr = hm(c,j) - hfus*(wice0(c,j)-h2osoi_ice(c,j))/dtime + xm2(c,j) = xm(c,j) - h2osoi_ice(c,j) !excess ice melting + if (h2osoi_ice(c,j) == 0._r8) then ! this might be redundant + if (excess_ice(c,j) >= 0._r8 .and. xm2(c,j)>0._r8 .and. j>=2) then ! if there is excess ice to melt + excess_ice(c,j) = max(0._r8,wexice0(c,j) - xm2(c,j)) + heatr = hm(c,j) - hfus * (wexice0(c,j)-excess_ice(c,j)+wice0(c,j)-h2osoi_ice(c,j)) / dtime + endif + endif !end of excess ice block else if (xm(c,j) < 0._r8) then if (j <= 0) then h2osoi_ice(c,j) = min(wmass0(c,j), wice0(c,j)-xm(c,j)) ! snow else - if (wmass0(c,j) < supercool(c,j)) then + if (wmass0(c,j) - wexice0(c,j) < supercool(c,j)) then ! even if excess ice is present, it cannot refreeze h2osoi_ice(c,j) = 0._r8 else - h2osoi_ice(c,j) = min(wmass0(c,j) - supercool(c,j),wice0(c,j)-xm(c,j)) + h2osoi_ice(c,j) = min(wmass0(c,j) - wexice0(c,j) - supercool(c,j),wice0(c,j)-xm(c,j)) endif endif heatr = hm(c,j) - hfus*(wice0(c,j)-h2osoi_ice(c,j))/dtime endif - h2osoi_liq(c,j) = max(0._r8,wmass0(c,j)-h2osoi_ice(c,j)) + h2osoi_liq(c,j) = max(0._r8,wmass0(c,j)-h2osoi_ice(c,j)-excess_ice(c,j)) !melted excess ice is added to the respective soil layers + if (abs(heatr) > 0._r8) then if (j == snl(c)+1) then @@ -1373,11 +1447,14 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & endif ! end of heatr > 0 if-block if (j >= 1) then - xmf(c) = xmf(c) + hfus*(wice0(c,j)-h2osoi_ice(c,j))/dtime + xmf(c) = xmf(c) + hfus*(wice0(c,j)-h2osoi_ice(c,j))/dtime + & + hfus*(wexice0(c,j)-excess_ice(c,j))/dtime + ! subsidence calculation + exice_subs_col(c,j) = max(0._r8, (wexice0(c,j)-excess_ice(c,j))/denice) else xmf(c) = xmf(c) + hfus*(wice0(c,j)-h2osoi_ice(c,j))/dtime endif - + if (imelt(c,j) == 1 .AND. j < 1) then qflx_snomelt_lyr(c,j) = max(0._r8,(wice0(c,j)-h2osoi_ice(c,j)))/dtime qflx_snomelt(c) = qflx_snomelt(c) + qflx_snomelt_lyr(c,j) @@ -1398,6 +1475,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, & end do ! end of column-loop enddo ! end of level-loop + ! Needed for history file output do fc = 1,num_nolakec diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index f2d9317f82..d20fc862b1 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -647,11 +647,11 @@ subroutine InitCold(this, bounds, & ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_TKFRZ - use clm_varcon , only : denice, denh2o, sb - use landunit_varcon, only : istwet, istsoil, istdlak, istice + use clm_varcon , only : denice, denh2o + use landunit_varcon, only : istwet, istsoil, istdlak, istice, istcrop use column_varcon , only : icol_road_imperv, icol_roof, icol_sunwall use column_varcon , only : icol_shadewall, icol_road_perv - use clm_varctl , only : iulog, use_vancouver, use_mexicocity + use clm_varctl , only : iulog, use_vancouver, use_mexicocity, use_excess_ice ! ! !ARGUMENTS: class(temperature_type) :: this @@ -742,7 +742,9 @@ subroutine InitCold(this, bounds, & end if else this%t_soisno_col(c,1:nlevgrnd) = 274._r8 - + if (use_excess_ice .and. (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop)) then + this%t_soisno_col(c,1:nlevgrnd) = SHR_CONST_TKFRZ - 5.0_r8 !needs to be below freezing to properly initiate excess ice + end if endif endif end do diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index bfeea81949..885222f33b 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -358,7 +358,8 @@ subroutine AccumulateSoilLiqIceMassNonLake(bounds, num_c, filter_c, & associate( & h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) - h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + excess_ice => waterstate_inst%excess_ice_col & ! Input [real(r8) (:,:) ] excess ice (kg/m2) ) do j = 1, nlevmaxurbgrnd @@ -382,7 +383,7 @@ subroutine AccumulateSoilLiqIceMassNonLake(bounds, num_c, filter_c, & if (has_h2o) then liquid_mass(c) = liquid_mass(c) + h2osoi_liq(c,j) - ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) + ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) + excess_ice(c,j) end if end do end do diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 7804fa3746..9ecf7208a7 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -17,7 +17,7 @@ module WaterDiagnosticBulkType use decompMod , only : bounds_type use abortutils , only : endrun use clm_varctl , only : use_cn, iulog, use_luna - use clm_varpar , only : nlevgrnd, nlevsno, nlevcan + use clm_varpar , only : nlevgrnd, nlevsno, nlevcan, nlevsoi use clm_varcon , only : spval use LandunitType , only : lun use ColumnType , only : col @@ -48,6 +48,10 @@ module WaterDiagnosticBulkType real(r8), pointer :: air_vol_col (:,:) ! col air filled porosity real(r8), pointer :: h2osoi_liqvol_col (:,:) ! col volumetric liquid water content (v/v) real(r8), pointer :: swe_old_col (:,:) ! col initial snow water + real(r8), pointer :: exice_subs_tot_col (:) ! col total subsidence due to excess ice melt (m) + real(r8), pointer :: exice_subs_col (:,:) ! col per layer subsidence due to excess ice melt (m) + real(r8), pointer :: exice_vol_col (:,:) ! col per layer volumetric excess ice content (m3/m3) + real(r8), pointer :: exice_vol_tot_col (:) ! col averaged volumetric excess ice content (m3/m3) real(r8), pointer :: snw_rds_col (:,:) ! col snow grain radius (col,lyr) [m^-6, microns] real(r8), pointer :: snw_rds_top_col (:) ! col snow grain radius (top layer) [m^-6, microns] @@ -166,6 +170,7 @@ subroutine InitBulkAllocate(this, bounds) ! ! !USES: use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + use clm_varpar , only : nlevmaxurbgrnd ! ! !ARGUMENTS: class(waterdiagnosticbulk_type), intent(inout) :: this @@ -194,6 +199,10 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%h2osoi_ice_tot_col (begc:endc)) ; this%h2osoi_ice_tot_col (:) = nan allocate(this%h2osoi_liq_tot_col (begc:endc)) ; this%h2osoi_liq_tot_col (:) = nan allocate(this%swe_old_col (begc:endc,-nlevsno+1:0)) ; this%swe_old_col (:,:) = nan + allocate(this%exice_subs_tot_col (begc:endc)) ; this%exice_subs_tot_col (:) = 0.0_r8 + allocate(this%exice_subs_col (begc:endc, 1:nlevmaxurbgrnd)) ; this%exice_subs_col (:,:) = 0.0_r8 + allocate(this%exice_vol_col (begc:endc, 1:nlevsoi)) ; this%exice_vol_col (:,:) = 0.0_r8 + allocate(this%exice_vol_tot_col (begc:endc)) ; this%exice_vol_tot_col (:) = 0.0_r8 allocate(this%snw_rds_col (begc:endc,-nlevsno+1:0)) ; this%snw_rds_col (:,:) = nan allocate(this%snw_rds_top_col (begc:endc)) ; this%snw_rds_top_col (:) = nan @@ -233,6 +242,7 @@ subroutine InitBulkHistory(this, bounds) ! !USES: use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use histFileMod , only : hist_addfld1d, hist_addfld2d, no_snow_normal, no_snow_zero + use clm_varctl , only : use_excess_ice ! ! !ARGUMENTS: class(waterdiagnosticbulk_type), intent(in) :: this @@ -277,8 +287,28 @@ subroutine InitBulkHistory(this, bounds) fname=this%info%fname('TOTSOILICE'), & units='kg/m2', & avgflag='A', & - long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & + long_name=this%info%lname('vertically summed soil ice (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') + ! excess ice vars + if (use_excess_ice) then + this%exice_vol_tot_col(begc:endc) = 0.0_r8 + call hist_addfld1d ( & + fname=this%info%fname('TOTEXICE_VOL'), & + units='m3/m3', & + avgflag='A', & + l2g_scale_type='veg', & + long_name=this%info%lname('vertically averaged volumetric excess ice concentration (veg landunits only)'), & + ptr_col=this%exice_vol_tot_col) + + this%exice_subs_tot_col(begc:endc) = 0.0_r8 + call hist_addfld1d ( & + fname=this%info%fname('SUBSIDENCE'), & + units='m', & + avgflag='SUM', & + l2g_scale_type='veg', & + long_name=this%info%lname('subsidence due to excess ice melt (veg landunits only)'), & + ptr_col=this%exice_subs_tot_col) + end if this%iwue_ln_patch(begp:endp) = spval call hist_addfld1d ( & @@ -741,9 +771,10 @@ subroutine RestartBulk(this, bounds, ncid, flag, writing_finidat_interp_dest_fil use spmdMod , only : masterproc use clm_varcon , only : pondmx, watmin, spval, nameg use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall - use clm_varctl , only : bound_h2osoi + use clm_varctl , only : bound_h2osoi, use_excess_ice, nsrest, nsrContinue use ncdio_pio , only : file_desc_t, ncd_io, ncd_double use restUtilMod + use ExcessIceStreamType, only : UseExcessIceStreams ! ! !ARGUMENTS: class(waterdiagnosticbulk_type), intent(inout) :: this @@ -755,6 +786,7 @@ subroutine RestartBulk(this, bounds, ncid, flag, writing_finidat_interp_dest_fil ! ! !LOCAL VARIABLES: logical :: readvar + logical :: excess_ice_on_restart !------------------------------------------------------------------------ @@ -875,7 +907,52 @@ subroutine RestartBulk(this, bounds, ncid, flag, writing_finidat_interp_dest_fil interpinic_flag='interp', readvar=readvar, data=this%wf_col) end if - + if (.not. use_excess_ice) then + ! no need to even define the restart vars + this%exice_subs_tot_col(bounds%begc:bounds%endc)=0.0_r8 + this%exice_vol_tot_col(bounds%begc:bounds%endc)=0.0_r8 + this%exice_subs_col(bounds%begc:bounds%endc,1:nlevgrnd)=0.0_r8 + this%exice_vol_col(bounds%begc:bounds%endc,1:nlevsoi)=0.0_r8 + else + ! initialization of these to zero is ok, since they might not be in the restart file + this%exice_subs_col(bounds%begc:bounds%endc,1:nlevgrnd)=0.0_r8 + this%exice_vol_col(bounds%begc:bounds%endc,1:nlevsoi)=0.0_r8 + call RestartExcessIceIssue( & + ncid = ncid, & + flag = flag, & + excess_ice_on_restart = excess_ice_on_restart) + ! have to at least define them + call restartvar(ncid=ncid, flag=flag, varname=this%info%fname('SUBSIDENCE'), & + dim1name='column', xtype=ncd_double, & + long_name=this%info%lname('vertically summed volumetric excess ice concentration (veg landunits only)'), & + units='m', & + interpinic_flag='interp', readvar=readvar, data=this%exice_subs_tot_col) + if (flag == 'read' .and. ((.not. readvar) .or. (.not. excess_ice_on_restart)) ) then ! when reading restart that does not have excess ice in it + if (nsrest == nsrContinue) then + call endrun(msg = "On a continue run, excess ice fields MUST be on the restart file "// & + errMsg(sourcefile, __LINE__)) + else if ( .not. UseExcessIceStreams() )then + call endrun(msg = "This input initial conditions file does NOT include excess ice fields" // & + ", and use_excess_ice_streams is off, one or the other needs to be changed "// & + errMsg(sourcefile, __LINE__)) + end if + this%exice_subs_tot_col(bounds%begc:bounds%endc)=0.0_r8 + this%exice_vol_tot_col(bounds%begc:bounds%endc)=0.0_r8 + this%exice_subs_col(bounds%begc:bounds%endc,1:nlevgrnd)=0.0_r8 + this%exice_vol_col(bounds%begc:bounds%endc,1:nlevsoi)=0.0_r8 + endif + call restartvar(ncid=ncid, flag=flag, varname=this%info%fname('TOTEXICE_VOL'), & + dim1name='column', xtype=ncd_double, & + long_name=this%info%lname('vertically averaged volumetric excess ice concentration (veg landunits only)'), & + units='m3/m3', & + interpinic_flag='interp', readvar=readvar, data=this%exice_vol_tot_col) + if (flag == 'read' .and. ((.not. readvar) .or. (.not. excess_ice_on_restart)) ) then ! when reading restart that does not have excess ice in it + if (nsrest == nsrContinue) then + call endrun(msg = "On a continue run, excess ice fields MUST be on the restart file "// & + errMsg(sourcefile, __LINE__)) + end if + end if + endif end subroutine RestartBulk @@ -981,7 +1058,8 @@ subroutine Summary(this, bounds, & ! Compute end-of-timestep summaries of water diagnostic terms ! ! !USES: - use clm_varpar , only : nlevsoi + use clm_varcon , only : denice + use landunit_varcon, only : istsoil, istcrop ! !ARGUMENTS: class(waterdiagnosticbulk_type) , intent(inout) :: this type(bounds_type) , intent(in) :: bounds @@ -997,6 +1075,8 @@ subroutine Summary(this, bounds, & ! !LOCAL VARIABLES: integer :: fp, p, j, l, fc, c ! Indices real(r8):: fracl ! fraction of soil layer contributing to 10cm total soil water + real(r8):: dz_ext ! extended layer thickness due to excess ice + real(r8):: dz_tot ! total depth with extended thicknesses character(len=*), parameter :: subname = 'Summary' !----------------------------------------------------------------------- @@ -1006,10 +1086,15 @@ subroutine Summary(this, bounds, & h2osoi_ice => waterstate_inst%h2osoi_ice_col, & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_liq => waterstate_inst%h2osoi_liq_col, & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) + excess_ice => waterstate_inst%excess_ice_col, & ! Input: [real(r8) (:,:) ] excess ice lenses (kg/m2) (new) (1:nlevgrnd) + exice_subs_col => this%exice_subs_col , & ! Output: [real(r8) (:,:) ] per layer subsidence due to excess ice melt (m) + exice_vol_col => this%exice_vol_col , & ! Output: [real(r8) (:,:) ] per layer volumetric excess ice content (m3/m3) h2osoi_ice_tot => this%h2osoi_ice_tot_col , & ! Output: [real(r8) (:) ] vertically summed ice lens (kg/m2) h2osoi_liq_tot => this%h2osoi_liq_tot_col , & ! Output: [real(r8) (:) ] vertically summed liquid water (kg/m2) - h2osoi_liqice_10cm => this%h2osoi_liqice_10cm_col & ! Output: [real(r8) (:) ] liquid water + ice lens in top 10cm of soil (kg/m2) + h2osoi_liqice_10cm => this%h2osoi_liqice_10cm_col , & ! Output: [real(r8) (:) ] liquid water + ice lens in top 10cm of soil (kg/m2) + exice_subs_tot_col => this%exice_subs_tot_col , & ! Output [real(r8) (:) ] vertically summed subsidence due to excess ice melt (m) + exice_vol_tot_col => this%exice_vol_tot_col & ! Output [real(r8) (:) ] vertically averaged volumetric excess ice content (m3/m3) ) call this%waterdiagnostic_type%Summary(bounds, & @@ -1042,6 +1127,7 @@ subroutine Summary(this, bounds, & h2osoi_liqice_10cm(c) = 0.0_r8 h2osoi_liq_tot(c) = 0._r8 h2osoi_ice_tot(c) = 0._r8 + exice_subs_tot_col(c) = 0._r8 end if end do do j = 1, nlevsoi @@ -1064,10 +1150,31 @@ subroutine Summary(this, bounds, & end if h2osoi_liq_tot(c) = h2osoi_liq_tot(c) + h2osoi_liq(c,j) h2osoi_ice_tot(c) = h2osoi_ice_tot(c) + h2osoi_ice(c,j) + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + exice_subs_tot_col(c) = exice_subs_tot_col(c) + exice_subs_col(c,j) + endif end if end do end do + do fc = 1, num_nolakec ! extra loop needed since the one above has outer loop with layers + c = filter_nolakec(fc) + l = col%landunit(c) + if (.not. lun%urbpoi(l)) then + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + dz_tot = 0.0_r8 + exice_vol_tot_col(c)=0.0_r8 + do j = 1, nlevsoi + dz_ext = dz(c,j)+excess_ice(c,j)/denice + exice_vol_col(c,j)=excess_ice(c,j)/(denice*dz_ext) + dz_tot=dz_tot+dz_ext + exice_vol_tot_col(c)=exice_vol_tot_col(c)+exice_vol_col(c,j)*dz_ext ! (m) + enddo + exice_vol_tot_col(c)=exice_vol_tot_col(c)/dz_tot ! (m3/m3) + end if + end if + end do + end associate end subroutine Summary diff --git a/src/biogeophys/WaterStateBulkType.F90 b/src/biogeophys/WaterStateBulkType.F90 index 5c0298c8d5..ba3f0513c5 100644 --- a/src/biogeophys/WaterStateBulkType.F90 +++ b/src/biogeophys/WaterStateBulkType.F90 @@ -47,17 +47,17 @@ module WaterStateBulkType !------------------------------------------------------------------------ subroutine InitBulk(this, bounds, info, vars, & - h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer) + h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, NLFilename) class(waterstatebulk_type), intent(inout) :: this type(bounds_type) , intent(in) :: bounds class(water_info_base_type), intent(in), target :: info type(water_tracer_container_type), intent(inout) :: vars real(r8) , intent(in) :: h2osno_input_col(bounds%begc:) - real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) + real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) - logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run - + logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run + character(len=*) , intent(in) :: NLFilename ! Namelist filename call this%Init(bounds = bounds, & info = info, & @@ -65,7 +65,8 @@ subroutine InitBulk(this, bounds, info, vars, & h2osno_input_col = h2osno_input_col, & watsat_col = watsat_col, & t_soisno_col = t_soisno_col, & - use_aquifer_layer = use_aquifer_layer) + use_aquifer_layer = use_aquifer_layer, & + NLFilename = NLFilename) call this%InitBulkAllocate(bounds) @@ -187,7 +188,7 @@ end subroutine InitBulkCold !------------------------------------------------------------------------ subroutine RestartBulk(this, bounds, ncid, flag, & - watsat_col) + watsat_col, t_soisno_col, altmax_lastyear_indx) ! ! !DESCRIPTION: ! Read/Write module information to/from restart file. @@ -199,9 +200,11 @@ subroutine RestartBulk(this, bounds, ncid, flag, & ! !ARGUMENTS: class(waterstatebulk_type), intent(in) :: this type(bounds_type), intent(in) :: bounds - type(file_desc_t), intent(inout) :: ncid ! netcdf id - character(len=*) , intent(in) :: flag ! 'read' or 'write' - real(r8) , intent(in) :: watsat_col (bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) + type(file_desc_t), intent(inout) :: ncid ! netcdf id + character(len=*) , intent(in) :: flag ! 'read' or 'write' + real(r8) , intent(in) :: watsat_col (bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) + real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) + integer , intent(in) :: altmax_lastyear_indx(bounds%begc:) !col active layer index last year ! ! !LOCAL VARIABLES: integer :: c,l,j @@ -211,7 +214,9 @@ subroutine RestartBulk(this, bounds, ncid, flag, & SHR_ASSERT_ALL_FL((ubound(watsat_col) == (/bounds%endc,nlevmaxurbgrnd/)) , sourcefile, __LINE__) call this%restart (bounds, ncid, flag=flag, & - watsat_col=watsat_col(bounds%begc:bounds%endc,:)) + watsat_col=watsat_col(bounds%begc:bounds%endc,:), & + t_soisno_col=t_soisno_col(bounds%begc:, -nlevsno+1:), & + altmax_lastyear_indx=altmax_lastyear_indx(bounds%begc:)) call restartvar(ncid=ncid, flag=flag, & diff --git a/src/biogeophys/WaterStateType.F90 b/src/biogeophys/WaterStateType.F90 index 751f633875..cdbefa2a04 100644 --- a/src/biogeophys/WaterStateType.F90 +++ b/src/biogeophys/WaterStateType.F90 @@ -13,7 +13,8 @@ module WaterStateType use abortutils , only : endrun use decompMod , only : bounds_type use decompMod , only : subgrid_level_patch, subgrid_level_column, subgrid_level_gridcell - use clm_varctl , only : use_bedrock, iulog + use clm_varctl , only : use_bedrock, use_excess_ice, iulog + use spmdMod , only : masterproc use clm_varctl , only : use_fates use clm_varpar , only : nlevgrnd, nlevsoi, nlevurb, nlevmaxurbgrnd, nlevsno use clm_varcon , only : spval @@ -22,6 +23,7 @@ module WaterStateType use WaterInfoBaseType, only : water_info_base_type use WaterTracerContainerType, only : water_tracer_container_type use WaterTracerUtils, only : AllocateVar1d, AllocateVar2d + use ExcessIceStreamType, only : excessicestream_type, UseExcessIceStreams ! implicit none save @@ -46,10 +48,15 @@ module WaterStateType ! For the following dynbal baseline variables: positive values are subtracted to ! avoid counting liquid water content of "virtual" states; negative values are added ! to account for missing states in the model. - real(r8), pointer :: dynbal_baseline_liq_col(:) ! baseline liquid water content subtracted from each column's total liquid water calculation (mm H2O) - real(r8), pointer :: dynbal_baseline_ice_col(:) ! baseline ice content subtracted from each column's total ice calculation (mm H2O) + real(r8), pointer :: dynbal_baseline_liq_col(:) ! baseline liquid water content subtracted from each column's total liquid water calculation (mm H2O) + real(r8), pointer :: dynbal_baseline_ice_col(:) ! baseline ice content subtracted from each column's total ice calculation (mm H2O) - real(r8) :: aquifer_water_baseline ! baseline value for water in the unconfined aquifer (wa_col) for this bulk / tracer (mm) + real(r8) :: aquifer_water_baseline ! baseline value for water in the unconfined aquifer (wa_col) for this bulk / tracer (mm) + + real(r8), pointer :: excess_ice_col (:,:) ! col excess ice (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: exice_bulk_init (:) ! inital value for excess ice (new) (unitless) + + type(excessicestream_type), private :: exicestream ! stream type for excess ice initialization NUOPC only contains @@ -72,7 +79,7 @@ module WaterStateType !------------------------------------------------------------------------ subroutine Init(this, bounds, info, tracer_vars, & - h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer) + h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, NLFilename) class(waterstate_type), intent(inout) :: this type(bounds_type) , intent(in) :: bounds @@ -82,18 +89,19 @@ subroutine Init(this, bounds, info, tracer_vars, & real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run + character(len=*) , intent(in) :: NLFilename ! Namelist filename this%info => info call this%InitAllocate(bounds, tracer_vars) call this%InitHistory(bounds, use_aquifer_layer) - call this%InitCold(bounds = bounds, & - h2osno_input_col = h2osno_input_col, & - watsat_col = watsat_col, & - t_soisno_col = t_soisno_col, & - use_aquifer_layer = use_aquifer_layer) + h2osno_input_col = h2osno_input_col, & + watsat_col = watsat_col, & + t_soisno_col = t_soisno_col, & + use_aquifer_layer = use_aquifer_layer, & + NLFilename = NLFilename) end subroutine Init @@ -150,6 +158,14 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%dynbal_baseline_ice_col, name = 'dynbal_baseline_ice_col', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_column) + !excess ice vars + call AllocateVar2d(var = this%excess_ice_col, name = 'excess_ice_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = subgrid_level_column, & + dim2beg = -nlevsno+1, dim2end = nlevmaxurbgrnd) + call AllocateVar1d(var = this%exice_bulk_init, name = 'exice_bulk_init', & + container = tracer_vars, & + bounds = bounds, subgrid_level = subgrid_level_column) end subroutine InitAllocate @@ -268,6 +284,15 @@ subroutine InitHistory(this, bounds, use_aquifer_layer) ptr_col=this%wa_col, l2g_scale_type='veg') end if + ! Add excess ice fields to history + + if (use_excess_ice) then + data2dptr => this%excess_ice_col(begc:endc,1:nlevsoi) + call hist_addfld2d (fname='EXCESS_ICE', units='kg/m2', type2d='levsoi', & + avgflag='A', long_name='excess soil ice (vegetated landunits only)', & + ptr_col=this%excess_ice_col, l2g_scale_type='veg', default = 'inactive') + end if + ! (rgk 02-02-2017) There is intentionally no entry here for stored plant water ! I think that since the value is zero in all cases except ! for FATES plant hydraulics, it will be confusing for users @@ -281,7 +306,7 @@ end subroutine InitHistory !----------------------------------------------------------------------- subroutine InitCold(this, bounds, & - h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer) + h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, NLFilename) ! ! !DESCRIPTION: ! Initialize time constant variables and cold start conditions @@ -290,20 +315,22 @@ subroutine InitCold(this, bounds, & use shr_const_mod , only : SHR_CONST_TKFRZ use landunit_varcon , only : istwet, istsoil, istcrop, istice use column_varcon , only : icol_road_perv, icol_road_imperv - use clm_varcon , only : denice, denh2o, bdsno + use clm_varcon , only : denice, denh2o, bdsno , zisoi use clm_varcon , only : tfrz, aquifer_water_baseline + use initVerticalMod , only : find_soil_layer_containing_depth ! ! !ARGUMENTS: class(waterstate_type), intent(inout) :: this type(bounds_type) , intent(in) :: bounds real(r8) , intent(in) :: h2osno_input_col(bounds%begc:) - real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) + real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) - logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run + logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run + character(len=*) , intent(in) :: NLFilename ! Namelist filename ! ! !LOCAL VARIABLES: - integer :: c,j,l,nlevs - integer :: nbedrock + integer :: c,j,l,nlevs,g + integer :: nbedrock, n05m ! layer containing 0.5 m real(r8) :: ratio !----------------------------------------------------------------------- @@ -340,7 +367,7 @@ subroutine InitCold(this, bounds, & if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then nlevs = nlevgrnd do j = 1, nlevs - if (use_bedrock) then + if (use_bedrock .and. col%nbedrock(c) <=nlevsoi) then nbedrock = col%nbedrock(c) else nbedrock = nlevsoi @@ -505,39 +532,91 @@ subroutine InitCold(this, bounds, & this%dynbal_baseline_liq_col(bounds%begc:bounds%endc) = 0._r8 this%dynbal_baseline_ice_col(bounds%begc:bounds%endc) = 0._r8 + !Initialize excess ice + if (use_excess_ice .and. NLFilename /= '') then + ! enforce initialization with 0 for everything + this%excess_ice_col(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd)=0.0_r8 + this%exice_bulk_init(bounds%begc:bounds%endc)=0.0_r8 + call this%exicestream%Init(bounds, NLFilename) ! get initial fraction of excess ice per column + ! + ! If excess ice is being read from streams, use the streams to + ! initialize + ! + if ( UseExcessIceStreams() )then + call this%exicestream%CalcExcessIce(bounds, this%exice_bulk_init) + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + l = col%landunit(c) + if (.not. lun%lakpoi(l)) then !not lake + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + if (zisoi(nlevsoi) >= 0.5_r8) then + call find_soil_layer_containing_depth(0.5_r8,n05m) + else + n05m=nlevsoi-1 + endif + if (use_bedrock .and. col%nbedrock(c) <=nlevsoi) then + nbedrock = col%nbedrock(c) + else + nbedrock = nlevsoi + endif + do j = 2, nlevmaxurbgrnd ! ignore first layer + if (n05m= n05m .and. jnlevsoi) then + nbedrock = col%nbedrock(c) + else + nbedrock = nlevsoi + end if + do j = 2, nlevmaxurbgrnd ! ignore first layer + if(altmax_lastyear_indx(c) < nbedrock) then + if (j>altmax_lastyear_indx(c) .and. j shr_kind_r8, CL => shr_kind_CL + use shr_log_mod , only : errMsg => shr_log_errMsg + use spmdMod , only : mpicom, masterproc + use clm_varctl , only : iulog + use abortutils , only : endrun + use decompMod , only : bounds_type + + ! !PUBLIC TYPES: + implicit none + private + + public :: UseExcessIceStreams ! If streams will be used + + type, public :: excessicestream_type + contains + + ! !PUBLIC MEMBER FUNCTIONS: + procedure, public :: Init ! Initialize and read data in + procedure, public :: CalcExcessIce ! Calculate excess ice ammount + + ! !PRIVATE MEMBER FUNCTIONS: + procedure, private :: ReadNML ! Read in namelist + + end type excessicestream_type + ! ! PRIVATE DATA: + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +!============================================================================== +contains +!============================================================================== + + subroutine Init(this, bounds, NLFilename) + ! + ! + ! arguments + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + + ! + ! local variables + + call this%ReadNML( bounds, NLFileName ) + end subroutine Init + + subroutine CalcExcessIce(this,bounds,exice_bulk_init) + + ! only transfers grid values to columns + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + real(r8) , intent(inout) :: exice_bulk_init(bounds%begc:bounds%endc) + ! + ! !LOCAL VARIABLES: + + end subroutine CalcExcessIce + + logical function UseExcessIceStreams() + ! + ! !DESCRIPTION: + ! Return true if + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + ! + ! !LOCAL VARIABLES: + UseExcessIceStreams = .false. +end function UseExcessIceStreams + +subroutine ReadNML(this, bounds, NLFilename) + ! + ! Read the namelist data stream information. + ! + ! Uses: + use shr_nl_mod , only : shr_nl_find_group_name + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_mpi_mod , only : shr_mpi_bcast + ! + ! arguments + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! local variables + integer :: nu_nml ! unit for namelist file + integer :: nml_error ! namelist i/o error flag + logical :: use_excess_ice_streams = .false. ! logical to turn on use of excess ice streams + character(len=CL) :: stream_fldFileName_exice = ' ' + character(len=CL) :: stream_mapalgo_exice = 'none' + character(len=*), parameter :: namelist_name = 'exice_streams' ! MUST agree with name in namelist and read + character(len=*), parameter :: subName = "('exice_streams::ReadNML')" + !----------------------------------------------------------------------- + + namelist /exice_streams/ & ! MUST agree with namelist_name above + stream_mapalgo_exice, stream_fldFileName_exice, use_excess_ice_streams + !----------------------------------------------------------------------- + ! Default values for namelist + + ! Read excess ice namelist + if (masterproc) then + open( newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) + call shr_nl_find_group_name(nu_nml, namelist_name, status=nml_error) + if (nml_error == 0) then + read(nu_nml, nml=exice_streams,iostat=nml_error) ! MUST agree with namelist_name above + if (nml_error /= 0) then + call endrun(msg=' ERROR reading '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + else + call endrun(msg=' ERROR finding '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + close(nu_nml) + endif + + call shr_mpi_bcast(use_excess_ice_streams , mpicom) + + if (masterproc) then + if ( use_excess_ice_streams ) then + call endrun(msg=' ERROR excess ice streams can NOT be on for the MCT driver'//errMsg(sourcefile, __LINE__)) + end if + if ( trim(stream_fldFileName_exice) /= '' ) then + call endrun(msg=' ERROR stream_fldFileName_exice can NOT be set for the MCT driver'//errMsg(sourcefile, __LINE__)) + end if + if ( trim(stream_mapalgo_exice) /= 'none' ) then + call endrun(msg=' ERROR stream_mapalgo_exice can only be none for the MCT driver'//errMsg(sourcefile, __LINE__)) + end if + endif + +end subroutine ReadNML + +end module ExcessIceStreamType diff --git a/src/cpl/share_esmf/ExcessIceStreamType.F90 b/src/cpl/share_esmf/ExcessIceStreamType.F90 new file mode 100644 index 0000000000..92d5632aff --- /dev/null +++ b/src/cpl/share_esmf/ExcessIceStreamType.F90 @@ -0,0 +1,325 @@ +module ExcessIceStreamType + +#include "shr_assert.h" + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Contains methods for reading in excess ice initial bulk values from data stream. + ! Needed in parameterization for excess ice in soil (Lee et al., 2014). + ! Used when use_excess_ice is true for initialization: + ! startup type runs starting from coldstart and initial datasets + ! that do not have required variables + ! or hybrid runs from cases with use_excess_ice was false. + ! Dataset is interpolated to 0.125x0.125 degrees grid from Brown et al., 1997 + ! with values derived from permafrost types. + ! Values represent fraction of excess ice within soil column + ! and are distributed within it later in initialization + ! + ! !USES + use ESMF + use dshr_strdata_mod , only : shr_strdata_type + use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_cl + use shr_log_mod , only : errMsg => shr_log_errMsg + use spmdMod , only : mpicom, masterproc + use clm_varctl , only : iulog + use abortutils , only : endrun + use decompMod , only : bounds_type + + ! !PUBLIC TYPES: + implicit none + private + + public :: UseExcessIceStreams ! If streams will be used + + type, public :: excessicestream_type + real(r8), pointer, private :: exice_bulk (:) ! excess ice bulk value (-) + contains + + ! !PUBLIC MEMBER FUNCTIONS: + procedure, public :: Init ! Initialize and read data in + procedure, public :: CalcExcessIce ! Calculate excess ice ammount + + ! !PRIVATE MEMBER FUNCTIONS: + procedure, private :: InitAllocate ! Allocate data + + end type excessicestream_type + ! ! PRIVATE DATA: + type, private :: streamcontrol_type + character(len=CL) :: stream_fldFileName_exice ! data Filename + character(len=CL) :: stream_meshfile_exice ! mesh Filename + character(len=CL) :: stream_mapalgo_exice ! map algo + contains + procedure, private :: ReadNML ! Read in namelist + end type streamcontrol_type + + logical :: namelist_read = .false. + type(streamcontrol_type), private :: control ! Stream control data + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +!============================================================================== +contains +!============================================================================== + + subroutine Init(this, bounds, NLFilename) + ! + use spmdMod , only : iam + use lnd_comp_shr , only : mesh, model_clock + use dshr_strdata_mod , only : shr_strdata_init_from_inline, shr_strdata_print + use dshr_strdata_mod , only : shr_strdata_advance + use dshr_methods_mod , only : dshr_fldbun_getfldptr + ! + ! arguments + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + + ! + ! local variables + integer :: ig, g, n ! Indices + integer :: year ! year (0, ...) for nstep+1 + integer :: mon ! month (1, ..., 12) for nstep+1 + integer :: day ! day of month (1, ..., 31) for nstep+1 + integer :: sec ! seconds into current date for nstep+1 + integer :: mcdate ! Current model date (yyyymmdd) + type(shr_strdata_type) :: sdat_exice ! input data stream + character(len=16), allocatable :: stream_varnames(:) ! array of stream field names + integer :: rc ! error code + real(r8), pointer :: dataptr1d(:) ! temporary pointer + character(len=*), parameter :: stream_name = 'excess ice' + + call this%InitAllocate( bounds ) + call control%ReadNML( bounds, NLFileName ) + if ( UseExcessIceStreams() )then + allocate(stream_varnames(1)) + stream_varnames = (/"EXICE"/) + + if (masterproc) then + write(iulog,*) ' stream_varnames = ',stream_varnames + write(iulog,*) ' Values will be used if the variable is not on the initial conditions dataset' + end if + + call shr_strdata_init_from_inline(sdat_exice, & + my_task = iam, & + logunit = iulog, & + compname = 'LND', & + model_clock = model_clock, & + model_mesh = mesh, & + stream_meshfile = control%stream_meshfile_exice, & + stream_lev_dimname = 'null', & + stream_mapalgo = control%stream_mapalgo_exice, & + stream_filenames = (/trim(control%stream_fldFileName_exice)/), & + stream_fldlistFile = stream_varnames, & + stream_fldListModel = stream_varnames, & + stream_yearFirst = 1996, & + stream_yearLast = 1996, & + stream_yearAlign = 1, & + stream_offset = 0, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + ! in ch4FinundatedStreamType it is set to linear but we have a single date dataset + stream_tintalgo = 'nearest', & + stream_name = 'excess ice ', & + rc = rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + !TODO + ! Explicitly set current date to a hardcoded constant value. Otherwise + ! using the real date can cause roundoff differences that are + ! detrected as issues with exact restart. EBK M05/20/2017 + ! call get_curr_date(year, mon, day, sec) + year = 1996 + mon = 12 + day = 31 + sec = 0 + mcdate = year*10000 + mon*100 + day + + call shr_strdata_advance(sdat_exice, ymd=mcdate, tod=sec, logunit=iulog, istr='exice', rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Get pointer for stream data that is time and spatially interpolate to model time and grid + do n = 1,size(stream_varnames) + call dshr_fldbun_getFldPtr(sdat_exice%pstrm(1)%fldbun_model, stream_varnames(n), fldptr1=dataptr1d, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + if (trim(stream_varnames(n)) == 'EXICE') then + ig = 0 + do g = bounds%begg,bounds%endg + ig = ig+1 + this%exice_bulk(g) = dataptr1d(ig) + end do + end if + end do + end if + end subroutine Init + + subroutine InitAllocate(this, bounds) + ! + ! !DESCRIPTION: + ! Allocate module variables and data structures + ! + ! !USES: + use shr_infnan_mod, only: nan => shr_infnan_nan, assignment(=) + ! + ! !ARGUMENTS: + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + integer :: begc, endc + integer :: begg, endg + !--------------------------------------------------------------------- + + begc = bounds%begc; endc = bounds%endc + begg = bounds%begg; endg = bounds%endg + + + allocate(this%exice_bulk(begg:endg)) ; this%exice_bulk(:) = nan + + end subroutine InitAllocate + + subroutine CalcExcessIce(this,bounds,exice_bulk_init) + + ! only transfers grid values to columns + use shr_const_mod , only : SHR_CONST_TKFRZ + use landunit_varcon , only : istwet, istsoil, istcrop, istice + use column_varcon , only : icol_road_perv, icol_road_imperv + use clm_varcon , only : denice + use clm_varcon , only : tfrz + use ColumnType , only : col + use LandunitType , only : lun + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + real(r8) , intent(inout) :: exice_bulk_init(bounds%begc:bounds%endc) + ! + ! !LOCAL VARIABLES: + integer :: begc, endc + integer :: begg, endg + integer :: c, l, g !counters + + exice_bulk_init(bounds%begc:bounds%endc)=0.0_r8 + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + l = col%landunit(c) + if ((.not. lun%lakpoi(l)) .and. (.not. lun%urbpoi(l)) .and. (.not. lun%itype(l) == istwet) .and. (.not. lun%itype(l) == istice)) then !not lake + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + exice_bulk_init(c)=this%exice_bulk(g) + else + exice_bulk_init(c) = 0.0_r8 + endif + else + exice_bulk_init(c)=0.0_r8 + endif + enddo + + end subroutine CalcExcessIce + + logical function UseExcessIceStreams() + ! + ! !DESCRIPTION: + ! Return true if + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + ! + ! !LOCAL VARIABLES: + if ( .not. namelist_read ) then + call endrun(msg=' ERROR UseExcessIceStreams being called, but namelist has not been read yet'//errMsg(sourcefile, __LINE__)) + end if + if ( trim(control%stream_fldFileName_exice) == '' )then + UseExcessIceStreams = .false. + else + UseExcessIceStreams = .true. + end if +end function UseExcessIceStreams + +subroutine ReadNML(this, bounds, NLFilename) + ! + ! Read the namelist data stream information. + ! + ! Uses: + use shr_nl_mod , only : shr_nl_find_group_name + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_mpi_mod , only : shr_mpi_bcast + ! + ! arguments + implicit none + class(streamcontrol_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! local variables + integer :: nu_nml ! unit for namelist file + integer :: nml_error ! namelist i/o error flag + logical :: use_excess_ice_streams = .false. ! logical to turn on use of excess ice streams + character(len=CL) :: stream_fldFileName_exice = ' ' + character(len=CL) :: stream_meshfile_exice = ' ' + character(len=CL) :: stream_mapalgo_exice = 'bilinear' + character(len=*), parameter :: namelist_name = 'exice_streams' ! MUST agree with name in namelist and read + character(len=*), parameter :: subName = "('exice_streams::ReadNML')" + !----------------------------------------------------------------------- + + namelist /exice_streams/ & ! MUST agree with namelist_name above + stream_mapalgo_exice, stream_fldFileName_exice, stream_meshfile_exice, use_excess_ice_streams + + ! Default values for namelist + + ! Read excess ice namelist + if (masterproc) then + open( newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) + call shr_nl_find_group_name(nu_nml, namelist_name, status=nml_error) + if (nml_error == 0) then + read(nu_nml, nml=exice_streams,iostat=nml_error) ! MUST agree with namelist_name above + if (nml_error /= 0) then + call endrun(msg=' ERROR reading '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + else + call endrun(msg=' ERROR finding '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + close(nu_nml) + endif + + call shr_mpi_bcast(use_excess_ice_streams , mpicom) + call shr_mpi_bcast(stream_mapalgo_exice , mpicom) + call shr_mpi_bcast(stream_fldFileName_exice , mpicom) + call shr_mpi_bcast(stream_meshfile_exice , mpicom) + + if (masterproc) then + write(iulog,*) ' ' + if ( use_excess_ice_streams ) then + write(iulog,*) 'excess ice streams are enabled: ' + write(iulog,*) namelist_name, ' stream settings:' + write(iulog,*) ' stream_fldFileName_exice = ',stream_fldFileName_exice + write(iulog,*) ' stream_meshfile_exice = ',stream_meshfile_exice + write(iulog,*) ' stream_mapalgo_exice = ',stream_mapalgo_exice + if ( trim(stream_fldFileName_exice) == '' )then + call endrun(msg=' ERROR excess ice streams are on, but stream_fldFileName_exice is NOT set'//errMsg(sourcefile, __LINE__)) + end if + else + write(iulog,*) 'excess ice streams are off' + if ( trim(stream_fldFileName_exice) /= '' )then + call endrun(msg=' ERROR excess ice streams are off, but stream_fldFileName_exice is set'//errMsg(sourcefile, __LINE__)) + end if + end if + endif + this%stream_fldFileName_exice = stream_fldFileName_exice + this%stream_meshfile_exice = stream_meshfile_exice + this%stream_mapalgo_exice = stream_mapalgo_exice + namelist_read = .true. + +end subroutine ReadNML + + +end module ExcessIceStreamType diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index 18447f2843..205bd7e8e2 100644 --- a/src/main/clm_instMod.F90 +++ b/src/main/clm_instMod.F90 @@ -486,6 +486,7 @@ subroutine clm_instRest(bounds, ncid, flag, writing_finidat_interp_dest_file) use ncdio_pio , only : file_desc_t use UrbanParamsType , only : IsSimpleBuildTemp, IsProgBuildTemp use decompMod , only : get_proc_bounds, get_proc_clumps, get_clump_bounds + use clm_varpar , only : nlevsno ! ! !DESCRIPTION: @@ -532,7 +533,9 @@ subroutine clm_instRest(bounds, ncid, flag, writing_finidat_interp_dest_file) call water_inst%restart(bounds, ncid, flag=flag, & writing_finidat_interp_dest_file = writing_finidat_interp_dest_file, & - watsat_col = soilstate_inst%watsat_col(bounds%begc:bounds%endc,:)) + watsat_col = soilstate_inst%watsat_col(bounds%begc:bounds%endc,:), & + t_soisno_col=temperature_inst%t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:), & + altmax_lastyear_indx=active_layer_inst%altmax_lastyear_indx_col(bounds%begc:bounds%endc)) call irrigation_inst%restart (bounds, ncid, flag=flag) diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index f54c0b942c..64aa40d94d 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -325,6 +325,11 @@ module clm_varctl real(r8), public :: soil_layerstruct_userdefined(99) = rundef integer, public :: soil_layerstruct_userdefined_nlevsoi = iundef + !---------------------------------------------------------- + !excess ice physics switch + !---------------------------------------------------------- + logical, public :: use_excess_ice = .false. ! true. => use excess ice physics + !---------------------------------------------------------- ! plant hydraulic stress switch !---------------------------------------------------------- diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 42cd289aba..c1a7959459 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -244,6 +244,8 @@ subroutine control_init(dtime) namelist /clm_inparm/ use_soil_moisture_streams + namelist /clm_inparm/ use_excess_ice + namelist /clm_inparm/ use_lai_streams namelist /clm_inparm/ use_bedrock @@ -733,6 +735,8 @@ subroutine control_spmd() call mpi_bcast (use_soil_moisture_streams, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (use_excess_ice, 1, MPI_LOGICAL, 0, mpicom,ier) + call mpi_bcast (use_lai_streams, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_cropcal_streams, 1, MPI_LOGICAL, 0, mpicom, ier) @@ -869,6 +873,7 @@ subroutine control_print () write(iulog,*) ' use_nitrif_denitrif = ', use_nitrif_denitrif write(iulog,*) ' use_extralakelayers = ', use_extralakelayers write(iulog,*) ' use_vichydro = ', use_vichydro + write(iulog,*) ' use_excess_ice = ', use_excess_ice write(iulog,*) ' use_cn = ', use_cn write(iulog,*) ' use_cndv = ', use_cndv write(iulog,*) ' use_crop = ', use_crop diff --git a/src/main/restFileMod.F90 b/src/main/restFileMod.F90 index d9c23d49f3..6a574406fd 100644 --- a/src/main/restFileMod.F90 +++ b/src/main/restFileMod.F90 @@ -27,6 +27,7 @@ module restFileMod use glcBehaviorMod , only : glc_behavior_type use reweightMod , only : reweight_wrapup use IssueFixedMetadataHandler, only : write_issue_fixed_metadata, read_issue_fixed_metadata + use restUtilMod , only : excess_ice_issue ! ! !PUBLIC TYPES: implicit none @@ -592,6 +593,8 @@ subroutine restFile_write_issues_fixed(ncid, writing_finidat_interp_dest_file) ! !DESCRIPTION: ! Write metadata for issues fixed ! + ! !USES: + use clm_varctl, only : use_excess_ice ! !ARGUMENTS: type(file_desc_t), intent(inout) :: ncid ! local file id logical , intent(in) :: writing_finidat_interp_dest_file ! true if we are writing a finidat_interp_dest file @@ -606,6 +609,15 @@ subroutine restFile_write_issues_fixed(ncid, writing_finidat_interp_dest_file) ncid = ncid, & writing_finidat_interp_dest_file = writing_finidat_interp_dest_file, & issue_num = lake_dynbal_baseline_issue) + ! If running with execess ice then mark the restart file as having excess ice fixed + ! This is a permanent feature, i.e. not expected to be removed from here. + ! It would only be removed if we decided to make use_excess_ice = .true. the default. + if ( use_excess_ice ) then + call write_issue_fixed_metadata( & + ncid = ncid, & + writing_finidat_interp_dest_file = writing_finidat_interp_dest_file, & + issue_num = excess_ice_issue) + end if end subroutine restFile_write_issues_fixed diff --git a/src/unit_test_shr/unittestWaterTypeFactory.F90 b/src/unit_test_shr/unittestWaterTypeFactory.F90 index cc3510e881..f5f417f9ce 100644 --- a/src/unit_test_shr/unittestWaterTypeFactory.F90 +++ b/src/unit_test_shr/unittestWaterTypeFactory.F90 @@ -161,7 +161,8 @@ subroutine create_water_type(this, water_inst, & h2osno_col = col_array(0._r8), & snow_depth_col = col_array(0._r8), & watsat_col = l_watsat_col, & - t_soisno_col = l_t_soisno_col) + t_soisno_col = l_t_soisno_col, & + NLFilename = ' ') end subroutine create_water_type subroutine teardown(this, water_inst) diff --git a/src/unit_test_stubs/CMakeLists.txt b/src/unit_test_stubs/CMakeLists.txt index 38abfb1633..2d7fe23378 100644 --- a/src/unit_test_stubs/CMakeLists.txt +++ b/src/unit_test_stubs/CMakeLists.txt @@ -1,6 +1,7 @@ add_subdirectory(csm_share) add_subdirectory(dyn_subgrid) add_subdirectory(main) +add_subdirectory(share_esmf) add_subdirectory(utils) sourcelist_to_parent(clm_sources) diff --git a/src/unit_test_stubs/share_esmf/CMakeLists.txt b/src/unit_test_stubs/share_esmf/CMakeLists.txt new file mode 100644 index 0000000000..5eb2d42415 --- /dev/null +++ b/src/unit_test_stubs/share_esmf/CMakeLists.txt @@ -0,0 +1,5 @@ +list(APPEND clm_sources + ExcessIceStreamType.F90 + ) + +sourcelist_to_parent(clm_sources) diff --git a/src/unit_test_stubs/share_esmf/ExcessIceStreamType.F90 b/src/unit_test_stubs/share_esmf/ExcessIceStreamType.F90 new file mode 100644 index 0000000000..60bac5ad28 --- /dev/null +++ b/src/unit_test_stubs/share_esmf/ExcessIceStreamType.F90 @@ -0,0 +1,79 @@ +module ExcessIceStreamType + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Stub module for Excess ice streams + ! + ! !USES + use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_cl + use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun + use decompMod , only : bounds_type + + ! !PUBLIC TYPES: + implicit none + private + + public :: UseExcessIceStreams ! If streams will be used + + type, public :: excessicestream_type + contains + + ! !PUBLIC MEMBER FUNCTIONS: + procedure, public :: Init ! Initialize and read data in + procedure, public :: CalcExcessIce ! Calculate excess ice ammount + end type excessicestream_type + ! ! PRIVATE DATA: + + logical :: namelist_read = .false. + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +!============================================================================== +contains +!============================================================================== + + subroutine Init(this, bounds, NLFilename) + ! + ! arguments + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + + namelist_read = .true. + + end subroutine Init + + subroutine CalcExcessIce(this,bounds,exice_bulk_init) + + ! only transfers grid values to columns + implicit none + class(excessicestream_type) :: this + type(bounds_type), intent(in) :: bounds + real(r8) , intent(inout) :: exice_bulk_init(bounds%begc:bounds%endc) + ! + ! !LOCAL VARIABLES: + call endrun(msg=' ERROR CalcExcessIce stub is being called and should NOT be'//errMsg(sourcefile, __LINE__)) + + end subroutine CalcExcessIce + + logical function UseExcessIceStreams() + ! + ! !DESCRIPTION: + ! Return true if + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + ! + ! !LOCAL VARIABLES: + if ( .not. namelist_read ) then + call endrun(msg=' ERROR UseExcessIceStreams being called, but namelist has not been read yet'//errMsg(sourcefile, __LINE__)) + end if + UseExcessIceStreams = .false. +end function UseExcessIceStreams + +end module ExcessIceStreamType diff --git a/src/unit_test_stubs/utils/restUtilMod_stub.F90.in b/src/unit_test_stubs/utils/restUtilMod_stub.F90.in index 2f57303c78..135c977bf6 100644 --- a/src/unit_test_stubs/utils/restUtilMod_stub.F90.in +++ b/src/unit_test_stubs/utils/restUtilMod_stub.F90.in @@ -21,6 +21,8 @@ module restUtilMod public :: restartvar + public :: RestartExcessIceIssue + public :: set_missing_from_template public :: CallRestartvarDimOK @@ -151,6 +153,29 @@ contains end subroutine set_missing_from_template + !----------------------------------------------------------------------- + subroutine RestartExcessIceIssue(ncid, flag, excess_ice_on_restart) + ! + ! !DESCRIPTION: + ! + ! !USES: + ! + ! !ARGUMENTS: + type(file_desc_t), intent(inout) :: ncid ! netcdf id + character(len=*) , intent(in) :: flag ! 'read' or 'write' + logical, intent(out) :: excess_ice_on_restart ! If excess ice is on the restart file + ! + ! !LOCAL VARIABLES: + integer :: attribute_value + + character(len=*), parameter :: subname = 'RestartExcessIceIssue' + !----------------------------------------------------------------------- + + excess_ice_on_restart = .false. + + end subroutine RestartExcessIceIssue + + !----------------------------------------------------------------------- logical function CallRestartvarDimOK (ncid, flag, dimname) ! diff --git a/src/utils/restUtilMod.F90.in b/src/utils/restUtilMod.F90.in index 32686c816c..1bece885e4 100644 --- a/src/utils/restUtilMod.F90.in +++ b/src/utils/restUtilMod.F90.in @@ -19,6 +19,7 @@ module restUtilMod implicit none save private + integer, parameter, public :: excess_ice_issue = 1787 ! save ! !----------------------------------------------------------------------- @@ -88,6 +89,8 @@ module restUtilMod end interface set_missing_vals_to_constant public :: set_missing_vals_to_constant + public :: RestartExcessIceIssue + private :: missing_field_possibly_abort private :: write_interpinic_flag @@ -744,6 +747,47 @@ contains end subroutine write_interpinic_flag + !----------------------------------------------------------------------- + subroutine RestartExcessIceIssue(ncid, flag, excess_ice_on_restart) + ! + ! !DESCRIPTION: + ! Is excess ice on the originating restart file? This is important to have + ! because the init_interp process copies the cold-start values to the + ! interpolated file if they aren't there, and we need to know if good values + ! exist on the originating restart file. + ! + ! !USES: + use ncdio_pio , only : file_desc_t + use IssueFixedMetadataHandler, only : read_issue_fixed_metadata + ! + ! !ARGUMENTS: + type(file_desc_t), intent(inout) :: ncid ! netcdf id + character(len=*) , intent(in) :: flag ! 'read' or 'write' + logical, intent(out) :: excess_ice_on_restart ! If excess ice is on the restart file + ! + ! !LOCAL VARIABLES: + integer :: attribute_value + + + character(len=*), parameter :: subname = 'RestartExcessIceIssue' + !----------------------------------------------------------------------- + + excess_ice_on_restart = .false. + ! The write of the issue metadata is in restFileMod::: restFile_write_issues_fixed + if (flag == 'read' )then + call read_issue_fixed_metadata( & + ncid = ncid, & + issue_num = excess_ice_issue, & + attribute_value = attribute_value) + if (attribute_value == 0) then + excess_ice_on_restart = .false. + else + excess_ice_on_restart = .true. + end if + + end if + + end subroutine RestartExcessIceIssue !----------------------------------------------------------------------- logical function CallRestartvarDimOK (ncid, flag, dimname) !