diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index ae53889fa2..3dea32a937 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2382,6 +2382,7 @@ sub setup_logic_dynamic_subgrid { setup_logic_do_transient_pfts($opts, $nl_flags, $definition, $defaults, $nl); setup_logic_do_transient_crops($opts, $nl_flags, $definition, $defaults, $nl); + setup_logic_do_transient_lakes($opts, $nl_flags, $definition, $defaults, $nl); setup_logic_do_harvest($opts, $nl_flags, $definition, $defaults, $nl); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'reset_dynbal_baselines'); @@ -2547,6 +2548,69 @@ sub setup_logic_do_transient_crops { } } +sub setup_logic_do_transient_lakes { + # + # Set do_transient_lakes default value, and perform error checking on do_transient_lakes + # + # Assumes the following are already set in the namelist (although it's okay + # for them to be unset if that will be their final state): + # - flanduse_timeseries + # + # NOTE(wjs, 2020-08-23) I based this function on setup_logic_do_transient_crops. I'm + # not sure if all of the checks here are truly important for transient lakes (in + # particular, my guess is that collapse_urban could probably be done with transient + # lakes - as well as transient pfts and transient crops for that matter), but some of + # the checks probably are needed, and it seems best to keep transient lakes consistent + # with other transient areas in this respect. + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + + my $var = 'do_transient_lakes'; + + # cannot_be_true will be set to a non-empty string in any case where + # do_transient_lakes should not be true; if it turns out that + # do_transient_lakes IS true in any of these cases, a fatal error will be + # generated + my $cannot_be_true = ""; + + my $n_dom_pfts = $nl->get_value( 'n_dom_pfts' ); + my $n_dom_landunits = $nl->get_value( 'n_dom_landunits' ); + my $toosmall_soil = $nl->get_value( 'toosmall_soil' ); + my $toosmall_crop = $nl->get_value( 'toosmall_crop' ); + my $toosmall_glacier = $nl->get_value( 'toosmall_glacier' ); + my $toosmall_lake = $nl->get_value( 'toosmall_lake' ); + my $toosmall_wetland = $nl->get_value( 'toosmall_wetland' ); + my $toosmall_urban = $nl->get_value( 'toosmall_urban' ); + + if (string_is_undef_or_empty($nl->get_value('flanduse_timeseries'))) { + $cannot_be_true = "$var can only be set to true when running a transient case (flanduse_timeseries non-blank)"; + } + + if (!$cannot_be_true) { + # Note that, if the variable cannot be true, we don't call add_default + # - so that we don't clutter up the namelist with variables that don't + # matter for this case + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var); + } + + # Make sure the value is false when it needs to be false - i.e., that the + # user hasn't tried to set a true value at an inappropriate time. + + if (&value_is_true($nl->get_value($var)) && $cannot_be_true) { + $log->fatal_error($cannot_be_true); + } + + # if do_transient_lakes is .true. and any of these (n_dom_* or toosmall_*) + # are > 0 or collapse_urban = .true., then give fatal error + if (&value_is_true($nl->get_value($var))) { + if (&value_is_true($nl->get_value('collapse_urban'))) { + $log->fatal_error("$var cannot be combined with collapse_urban"); + } + if ($n_dom_pfts > 0 || $n_dom_landunits > 0 || $toosmall_soil > 0 || $toosmall_crop > 0 || $toosmall_glacier > 0 || $toosmall_lake > 0 || $toosmall_wetland > 0 || $toosmall_urban > 0) { + $log->fatal_error("$var cannot be combined with any of the of the following > 0: n_dom_pfts > 0, n_dom_landunit > 0, toosmall_soil > 0._r8, toosmall_crop > 0._r8, toosmall_glacier > 0._r8, toosmall_lake > 0._r8, toosmall_wetland > 0._r8, toosmall_urban > 0._r8"); + } + } +} + sub setup_logic_do_harvest { # # Set do_harvest default value, and perform error checking on do_harvest diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 7f06341036..1ca1015d9e 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -3599,6 +3599,7 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts +.false. .false. diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 70520c926f..8e57bc98e4 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2424,6 +2424,12 @@ If TRUE, apply transient crops from flanduse_timeseries file. (Only valid for transient runs, where there is a flanduse_timeseries file.) + +If TRUE, apply transient lakes from flanduse_timeseries file. +(Only valid for transient runs, where there is a flanduse_timeseries file.) + + If TRUE, apply harvest from flanduse_timeseries file. diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index b51217eb7e..630ca766b2 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1622,6 +1622,16 @@ for ERS test as otherwise it won't work for a sub-day test" + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/include_user_mods new file mode 100644 index 0000000000..399579f425 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/include_user_mods @@ -0,0 +1 @@ +../monthly diff --git a/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm new file mode 100644 index 0000000000..b85d0ba6c1 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm @@ -0,0 +1,12 @@ +do_transient_lakes = .true. + +! This file was created with the following command: +! ncap2 -s 'PCT_LAKE=array(0.0,0.0,PCT_CROP); PCT_LAKE={0.,50.,25.,25.,25.,25.}; HASLAKE=array(1.,1.,AREA); PCT_CROP=array(0.0,0.0,PCT_LAKE); PCT_CROP={0.,25.,12.,12.,12.,12.}' landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_c160127.nc landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_dynLakes_c200928.nc +! Key points are that lake area starts as 0, increases after the first year, then decreases after the second year. +! PCT_CROP is also changed so that PCT_LAKE + PCT_CROP <= 100. (Here, PCT_CROP increases and decreases at the same time as PCT_LAKE in order to exercise the simultaneous increase or decrease of two landunits, but that isn't a critical part of this test.) +! Note that the use of this file means that this testmod can only be used with the 1x1_smallvilleIA grid. +flanduse_timeseries = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_dynLakes_c200928.nc' + +! BUG(wjs, 2020-09-25, ESCOMP/CTSM#43) Dynamic lakes don't work when methane is active, +! so for now disable methane for this test. +use_lch4 = .false. diff --git a/src/biogeochem/ch4Mod.F90 b/src/biogeochem/ch4Mod.F90 index 2432b85cf4..0155dd5019 100644 --- a/src/biogeochem/ch4Mod.F90 +++ b/src/biogeochem/ch4Mod.F90 @@ -1251,7 +1251,7 @@ subroutine DynamicColumnAdjustments(this, bounds, clump_index, column_state_upda character(len=*), parameter :: subname = 'DynamicColumnAdjustments' !----------------------------------------------------------------------- - ! BUG(wjs, 2016-02-16, bugz 2283) Need to do some special handling of finundated for + ! BUG(wjs, 2016-02-16, ESCOMP/CTSM#43) Need to do some special handling of finundated for ! increases in lake area, since lakes are assumed to be 100% inundated. Probably it's ! most appropriate for this special handling to happen elsewhere - i.e., within this ! routine, we do the standard adjustments as they are currently done, but then in the diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index c2d58f711e..69f7cdb167 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -20,6 +20,7 @@ module BalanceCheckMod use SoilHydrologyType , only : soilhydrology_type use SurfaceAlbedoType , only : surfalb_type use WaterStateType , only : waterstate_type + use LakestateType , only : lakestate_type use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type use WaterDiagnosticType, only : waterdiagnostic_type use Wateratm2lndType , only : wateratm2lnd_type @@ -122,7 +123,7 @@ end function GetBalanceCheckSkipSteps !----------------------------------------------------------------------- subroutine BeginWaterBalance(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & - water_inst, soilhydrology_inst, & + water_inst, soilhydrology_inst, lakestate_inst, & use_aquifer_layer) ! ! !DESCRIPTION: @@ -136,6 +137,7 @@ subroutine BeginWaterBalance(bounds, & integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points type(water_type) , intent(inout) :: water_inst + type(lakestate_type) , intent(in) :: lakestate_inst type(soilhydrology_type) , intent(in) :: soilhydrology_inst logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run ! @@ -150,6 +152,7 @@ subroutine BeginWaterBalance(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & soilhydrology_inst, & + lakestate_inst, & water_inst%bulk_and_tracers(i)%waterstate_inst, & water_inst%bulk_and_tracers(i)%waterdiagnostic_inst, & water_inst%bulk_and_tracers(i)%waterbalance_inst, & @@ -161,7 +164,8 @@ end subroutine BeginWaterBalance !----------------------------------------------------------------------- subroutine BeginWaterBalanceSingle(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & - soilhydrology_inst, waterstate_inst, waterdiagnostic_inst, waterbalance_inst, & + soilhydrology_inst, lakestate_inst, waterstate_inst, & + waterdiagnostic_inst, waterbalance_inst, & use_aquifer_layer) ! ! !DESCRIPTION: @@ -175,6 +179,7 @@ subroutine BeginWaterBalanceSingle(bounds, & integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points type(soilhydrology_type) , intent(in) :: soilhydrology_inst + type(lakestate_type) , intent(in) :: lakestate_inst class(waterstate_type) , intent(inout) :: waterstate_inst class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst class(waterbalance_type) , intent(inout) :: waterbalance_inst @@ -210,8 +215,8 @@ subroutine BeginWaterBalanceSingle(bounds, & water_mass = begwb(bounds%begc:bounds%endc)) call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, & - subtract_dynbal_baselines = .false., & + waterstate_inst, lakestate_inst, & + add_lake_water_and_subtract_dynbal_baselines = .false., & water_mass = begwb(bounds%begc:bounds%endc)) call waterstate_inst%CalculateTotalH2osno(bounds, num_nolakec, filter_nolakec, & diff --git a/src/biogeophys/CMakeLists.txt b/src/biogeophys/CMakeLists.txt index a25ba5ed0e..1bb52accc5 100644 --- a/src/biogeophys/CMakeLists.txt +++ b/src/biogeophys/CMakeLists.txt @@ -12,6 +12,7 @@ list(APPEND clm_sources InfiltrationExcessRunoffMod.F90 IrrigationMod.F90 LakeCon.F90 + LakeStateType.F90 QSatMod.F90 RootBiophysMod.F90 SaturatedExcessRunoffMod.F90 diff --git a/src/biogeophys/LakeHydrologyMod.F90 b/src/biogeophys/LakeHydrologyMod.F90 index 44588fafff..f6f83d8956 100644 --- a/src/biogeophys/LakeHydrologyMod.F90 +++ b/src/biogeophys/LakeHydrologyMod.F90 @@ -164,7 +164,7 @@ subroutine LakeHydrology(bounds, & t_soisno => temperature_inst%t_soisno_col , & ! Output: [real(r8) (:,:) ] snow temperature (Kelvin) dTdz_top => temperature_inst%dTdz_top_col , & ! Output: [real(r8) (:) ] temperature gradient in top layer K m-1] !TOD snot_top => temperature_inst%snot_top_col , & ! Output: [real(r8) (:) ] snow temperature in top layer [K] !TODO - t_sno_mul_mss => temperature_inst%t_sno_mul_mss_col , & ! Output: [real(r8) (:) ] col snow temperature multiplied by layer mass, layer sum (K * kg/m2) + t_sno_mul_mss => temperature_inst%t_sno_mul_mss_col , & ! Output: [real(r8) (:) ] col snow temperature multiplied by layer mass, layer sum (K * kg/m2) begwb => b_waterbalance_inst%begwb_col , & ! Input: [real(r8) (:) ] water mass begining of the time step endwb => b_waterbalance_inst%endwb_col , & ! Output: [real(r8) (:) ] water mass end of the time step @@ -649,8 +649,8 @@ subroutine LakeHydrology(bounds, & ! Determine ending water balance and volumetric soil water call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - b_waterstate_inst, & - subtract_dynbal_baselines = .false., & + b_waterstate_inst, lakestate_inst, & + add_lake_water_and_subtract_dynbal_baselines = .false., & water_mass = endwb(bounds%begc:bounds%endc)) do j = 1, nlevgrnd diff --git a/src/biogeophys/LakeTemperatureMod.F90 b/src/biogeophys/LakeTemperatureMod.F90 index f550f294a6..280020b6ef 100644 --- a/src/biogeophys/LakeTemperatureMod.F90 +++ b/src/biogeophys/LakeTemperatureMod.F90 @@ -134,6 +134,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la type(energyflux_type) , intent(inout) :: energyflux_inst type(temperature_type) , intent(inout) :: temperature_inst type(lakestate_type) , intent(inout) :: lakestate_inst + ! ! !LOCAL VARIABLES: real(r8), parameter :: p0 = 1._r8 ! neutral value of turbulent prandtl number @@ -247,7 +248,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la t_grnd => temperature_inst%t_grnd_col , & ! Input: [real(r8) (:) ] ground temperature (Kelvin) t_soisno => temperature_inst%t_soisno_col , & ! Output: [real(r8) (:,:) ] soil (or snow) temperature (Kelvin) t_lake => temperature_inst%t_lake_col , & ! Output: [real(r8) (:,:) ] col lake temperature (Kelvin) - + beta => lakestate_inst%betaprime_col , & ! Output: [real(r8) (:) ] col effective beta: sabg_lyr(p,jtop) for snow layers, beta otherwise lake_icefrac => lakestate_inst%lake_icefrac_col , & ! Output: [real(r8) (:,:) ] col mass fraction of lake layer that is frozen lake_icefracsurf => lakestate_inst%lake_icefracsurf_col , & ! Output: [real(r8) (:,:) ] col mass fraction of surface lake layer that is frozen @@ -999,8 +1000,10 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la ncvts(c) = ncvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) & + cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) fin(c) = fin(c) + phi(c,j) + end do end do + call waterstatebulk_inst%CalculateTotalH2osno(bounds, num_lakec, filter_lakec, & caller = 'LakeTemperature-2', & @@ -1022,6 +1025,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la end do + ! Check energy conservation. do fp = 1, num_lakep p = filter_lakep(fp) diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 5c69733f8a..d9ef9f5cb6 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -116,7 +116,7 @@ module TemperatureType real(r8), pointer :: xmf_h2osfc_col (:) ! latent heat of phase change of surface water real(r8), pointer :: fact_col (:,:) ! used in computing tridiagonal matrix real(r8), pointer :: c_h2osfc_col (:) ! heat capacity of surface water - + contains procedure, public :: Init @@ -620,7 +620,6 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) ptr_patch=this%t_veg10_night_patch, default='inactive') endif - end subroutine InitHistory !----------------------------------------------------------------------- diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 9488eedc70..cd7316af81 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -8,8 +8,8 @@ module TotalWaterAndHeatMod #include "shr_assert.h" use shr_kind_mod , only : r8 => shr_kind_r8 use decompMod , only : bounds_type - use clm_varcon , only : cpice, cpliq, denh2o, tfrz, hfus - use clm_varpar , only : nlevgrnd, nlevsoi, nlevurb + use clm_varcon , only : cpice, cpliq, denh2o, denice, tfrz, hfus + use clm_varpar , only : nlevgrnd, nlevsoi, nlevurb, nlevlak use ColumnType , only : col use LandunitType , only : lun use subgridAveMod , only : p2c @@ -21,6 +21,7 @@ module TotalWaterAndHeatMod use UrbanParamsType , only : urbanparams_type use SoilStateType , only : soilstate_type use TemperatureType , only : temperature_type + use LakeStateType , only : lakestate_type use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall use column_varcon , only : icol_road_perv, icol_road_imperv use landunit_varcon , only : istdlak, istsoil,istcrop,istwet,istice_mec @@ -38,12 +39,14 @@ module TotalWaterAndHeatMod ! routines parallel with the water routines. public :: ComputeWaterMassNonLake ! Compute total water mass of non-lake columns public :: ComputeWaterMassLake ! Compute total water mass of lake columns + public :: AccumulateLiqIceMassLake ! Accumulate lake water mass of lake columns, separated into liquid and ice public :: ComputeLiqIceMassNonLake ! Compute total water mass of non-lake columns, separated into liquid and ice public :: AccumulateSoilLiqIceMassNonLake ! Accumulate soil water mass of non-lake columns, separated into liquid and ice public :: ComputeLiqIceMassLake ! Compute total water mass of lake columns, separated into liquid and ice public :: ComputeHeatNonLake ! Compute heat content of non-lake columns public :: AccumulateSoilHeatNonLake ! Accumulate soil heat content of non-lake columns public :: ComputeHeatLake ! Compute heat content of lake columns + public :: AccumulateHeatLake ! Accumulate heat content of lake water of lake columns public :: AdjustDeltaHeatForDeltaLiq ! Adjusts the change in gridcell heat content due to land cover change to account for the implicit heat flux associated with delta_liq public :: LiquidWaterHeat ! Get the total heat content of some mass of liquid water at a given temperature @@ -104,8 +107,8 @@ subroutine ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, & class(waterstate_type) , intent(in) :: waterstate_inst class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst - ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When we can accept answer changes to methane, - ! remove this argument, always assuming it's true. + ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658 + ! is resolved, remove this argument, always assuming it's true. logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass real(r8) , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2) @@ -139,8 +142,8 @@ end subroutine ComputeWaterMassNonLake !----------------------------------------------------------------------- subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, & - subtract_dynbal_baselines, & + waterstate_inst, lakestate_inst, & + add_lake_water_and_subtract_dynbal_baselines, & water_mass) ! ! !DESCRIPTION: @@ -154,10 +157,14 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points class(waterstate_type) , intent(in) :: waterstate_inst + type(lakestate_type) , intent(in) :: lakestate_inst - ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When we can accept answer changes to methane, - ! remove this argument, always assuming it's true. - logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass + ! Whether to (1) add lake water/ice to total accounting, and (2) subtract + ! dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass + ! + ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658 + ! is resolved, remove this argument, always assuming it's true. + logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines real(r8) , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2) ! @@ -176,7 +183,8 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & num_lakec = num_lakec, & filter_lakec = filter_lakec, & waterstate_inst = waterstate_inst, & - subtract_dynbal_baselines = subtract_dynbal_baselines, & + lakestate_inst = lakestate_inst, & + add_lake_water_and_subtract_dynbal_baselines = add_lake_water_and_subtract_dynbal_baselines, & liquid_mass = liquid_mass(bounds%begc:bounds%endc), & ice_mass = ice_mass(bounds%begc:bounds%endc)) @@ -210,8 +218,8 @@ subroutine ComputeLiqIceMassNonLake(bounds, num_nolakec, filter_nolakec, & class(waterstate_type) , intent(in) :: waterstate_inst class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst - ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When we can accept answer changes to methane, - ! remove this argument, always assuming it's true. + ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658 + ! is resolved, remove this argument, always assuming it's true. logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! computed liquid water mass (kg m-2) @@ -382,8 +390,8 @@ end subroutine AccumulateSoilLiqIceMassNonLake !----------------------------------------------------------------------- subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, & - subtract_dynbal_baselines, & + waterstate_inst, lakestate_inst, & + add_lake_water_and_subtract_dynbal_baselines, & liquid_mass, ice_mass) ! ! !DESCRIPTION: @@ -400,17 +408,21 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points class(waterstate_type), intent(in) :: waterstate_inst + type(lakestate_type) , intent(in) :: lakestate_inst - ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When we can accept answer changes to methane, - ! remove this argument, always assuming it's true. - logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass + ! Whether to (1) add lake water/ice to total accounting, and (2) subtract + ! dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass + ! + ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658 + ! is resolved, remove this argument, always assuming it's true. + logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! computed liquid water mass (kg m-2) real(r8) , intent(inout) :: ice_mass( bounds%begc: ) ! computed ice mass (kg m-2) ! ! !LOCAL VARIABLES: - integer :: c, j, fc ! indices - + integer :: c, j, fc ! indices + character(len=*), parameter :: subname = 'ComputeLiqIceMassLake' !----------------------------------------------------------------------- @@ -419,10 +431,9 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl , & ! Input: [integer (:) ] negative number of snow layers - h2osno_no_layers => waterstate_inst%h2osno_no_layers_col , & ! Input: [real(r8) (:) ] snow water that is not resolved into layers (mm H2O) 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) dynbal_baseline_liq => waterstate_inst%dynbal_baseline_liq_col, & ! Input: [real(r8) (:) ] baseline liquid water content subtracted from each column's total liquid water calculation (mm H2O) dynbal_baseline_ice => waterstate_inst%dynbal_baseline_ice_col & ! Input: [real(r8) (:) ] baseline ice content subtracted from each column's total ice calculation (mm H2O) ) @@ -433,6 +444,34 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & ice_mass(c) = 0._r8 end do + ! ------------------------------------------------------------------------ + ! Start with some large terms that will often cancel (the negative baselines and the + ! lake water content): In order to maintain precision in the other terms, it can help if + ! we deal with these large, often-canceling terms first. (If we accumulated some + ! small terms, then added a big term and then subtracted a big term, we would have + ! lost some precision in the small terms.) + ! ------------------------------------------------------------------------ + + if (add_lake_water_and_subtract_dynbal_baselines) then + ! Subtract baselines set in initialization + do fc = 1, num_lakec + c = filter_lakec(fc) + liquid_mass(c) = liquid_mass(c) - dynbal_baseline_liq(c) + ice_mass(c) = ice_mass(c) - dynbal_baseline_ice(c) + end do + + ! Lake water content + call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, & + lakestate_inst, & + tracer_ratio = waterstate_inst%info%get_ratio(), & + liquid_mass = liquid_mass(bounds%begc:bounds%endc), & + ice_mass = ice_mass(bounds%begc:bounds%endc)) + end if + + ! ------------------------------------------------------------------------ + ! Now add some other terms + ! ------------------------------------------------------------------------ + ! Snow water content do fc = 1, num_lakec c = filter_lakec(fc) @@ -453,19 +492,70 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & end do end do - if (subtract_dynbal_baselines) then - ! Subtract baselines set in initialization - do fc = 1, num_lakec - c = filter_lakec(fc) - liquid_mass(c) = liquid_mass(c) - dynbal_baseline_liq(c) - ice_mass(c) = ice_mass(c) - dynbal_baseline_ice(c) - end do - end if - end associate end subroutine ComputeLiqIceMassLake + !----------------------------------------------------------------------- + subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & + lakestate_inst, tracer_ratio, liquid_mass, ice_mass) + ! + ! !DESCRIPTION: + ! Accumulate lake water mass of lake columns, separated into liquid and ice. + ! + ! Adds to any existing values in liquid_mass and ice_mass. + ! + ! Note: Changes to this routine should generally be accompanied by similar changes to + ! AccumulateHeatLake + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_c ! number of column points in column filter (should not include lake points; can be a subset of the no-lake filter) + integer , intent(in) :: filter_c(:) ! column filter (should not include lake points; can be a subset of the no-lake filter) + type(lakestate_type) , intent(in) :: lakestate_inst + real(r8) , intent(in) :: tracer_ratio ! for water tracers, standard ratio of this tracer to bulk water (1 for bulk water) + real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! accumulated liquid mass (kg m-2) + real(r8) , intent(inout) :: ice_mass( bounds%begc: ) ! accumulated ice mass (kg m-2) + ! + ! !LOCAL VARIABLES: + integer :: c, j, fc ! indices + real(r8) :: h2olak_liq ! liquid water content of lake layer [kg/m²] + real(r8) :: h2olak_ice ! ice water content of lake layer [kg/m²] + + character(len=*), parameter :: subname = 'AccumulateLiqIceMassLake' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL_FL((ubound(liquid_mass) == (/bounds%endc/)), sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(ice_mass) == (/bounds%endc/)), sourcefile, __LINE__) + + associate( & + lake_icefrac => lakestate_inst%lake_icefrac_col, & ! Input: [real(r8) (:,:) ] lake ice fraction + dz_lake => col%dz_lake & ! Input: [real(r8) (:,:) ] lake depth (m) + ) + + ! Lake water content + do j = 1, nlevlak + do fc = 1, num_c + c = filter_c(fc) + ! Lake water volume isn't tracked explicitly, so we calculate it from lake + ! depth. Because it isn't tracked explicitly, we also don't have any water + ! tracer information, so we assume a fixed, standard tracer ratio. + h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) * tracer_ratio + + ! use water density rather than ice density because lake layer depths are not + ! adjusted when the layer freezes + h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) * tracer_ratio + + liquid_mass(c) = liquid_mass(c) + h2olak_liq + ice_mass(c) = ice_mass(c) + h2olak_ice + end do + end do + + end associate + + end subroutine AccumulateLiqIceMassLake + + !----------------------------------------------------------------------- subroutine ComputeHeatNonLake(bounds, num_nolakec, filter_nolakec, & urbanparams_inst, soilstate_inst, & @@ -662,6 +752,11 @@ subroutine ComputeHeatNonLake(bounds, num_nolakec, filter_nolakec, & ! trying to account for the temperature of these baselines at all. This amounts to ! assuming that the baselines that we add / subtract are at the average temperature ! of the real liquid water in the column. + + ! This is different for lakes: there the virtual lake water's temperature is excluded + ! to avoid having it dominating the average temperature calculation + ! see note at top of the AccumulateHeatLake subroutine. + do fc = 1, num_nolakec c = filter_nolakec(fc) heat(c) = heat(c) - dynbal_baseline_heat(c) @@ -799,7 +894,7 @@ end subroutine AccumulateSoilHeatNonLake !----------------------------------------------------------------------- subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & - soilstate_inst, temperature_inst, waterstatebulk_inst, & + soilstate_inst, temperature_inst, waterstatebulk_inst, lakestate_inst, & heat, heat_liquid, cv_liquid) ! ! !DESCRIPTION: @@ -820,10 +915,13 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & type(soilstate_type) , intent(in) :: soilstate_inst type(temperature_type) , intent(in) :: temperature_inst type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst + type(lakestate_type) , intent(in) :: lakestate_inst + real(r8) , intent(inout) :: heat( bounds%begc: ) ! sum of heat content for all columns [J/m^2] real(r8) , intent(inout) :: heat_liquid( bounds%begc: ) ! sum of heat content for all columns: liquid water, excluding latent heat [J/m^2] real(r8) , intent(inout) :: cv_liquid( bounds%begc: ) ! sum of liquid heat capacity for all columns [J/(m^2 K)] + ! ! !LOCAL VARIABLES: integer :: fc @@ -848,7 +946,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & t_soisno => temperature_inst%t_soisno_col, & ! soil temperature (Kelvin) dynbal_baseline_heat => temperature_inst%dynbal_baseline_heat_col, & ! Input: [real(r8) (:) ] baseline heat content subtracted from each column's total heat calculation (J/m2) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col, & ! liquid water (kg/m2) - h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! frozen water (kg/m2) + h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! frozen water (kg/m2) ) do fc = 1, num_lakec @@ -861,6 +959,42 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & latent_heat_liquid(c) = 0._r8 end do + ! ------------------------------------------------------------------------ + ! Start with some large terms that will often cancel (the negative baselines and the + ! lake water content): In order to maintain precision in the other terms, it can help if + ! we deal with these large, often-canceling terms first. (If we accumulated some + ! small terms, then added a big term and then subtracted a big term, we would have + ! lost some precision in the small terms.) + ! ------------------------------------------------------------------------ + + ! Subtract baselines set in initialization + ! + ! NOTE(wjs, 2019-03-01) I haven't given enough thought to how (if at all) we should + ! correct for heat_liquid and cv_liquid, which are used to determine the weighted + ! average liquid water temperature. For example, if we're subtracting out a baseline + ! water amount because a particular water state is fictitious, we probably shouldn't + ! include that particular state when determining the weighted average temperature of + ! liquid water. And conversely, if we're adding a state via these baselines, should + ! we also add some water temperature of that state? The tricky thing here is what to + ! do when we end up subtracting water due to the baselines, so for now I'm simply not + ! trying to account for the temperature of these baselines at all. This amounts to + ! assuming that the baselines that we add / subtract are at the average temperature + ! of the real liquid water in the column. + do fc = 1, num_lakec + c = filter_lakec(fc) + heat(c) = -dynbal_baseline_heat(c) + end do + + ! Lake water heat content + ! Note that we do NOT accumulate heat_liquid and cv_liquid in this call. See the + ! comments at the top of AccumulateHeatLake for rationale. + call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, & + heat) + + ! ------------------------------------------------------------------------ + ! Now add some other terms + ! ------------------------------------------------------------------------ + ! Snow heat content do fc = 1, num_lakec c = filter_lakec(fc) @@ -881,7 +1015,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & end do end do - ! Soil water content of the soil under the lake + ! Soil heat content of the soil under the lake do j = 1,nlevgrnd do fc = 1, num_lakec c = filter_lakec(fc) @@ -898,37 +1032,116 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & TempToHeat(temp = t_soisno(c,j), cv = (h2osoi_ice(c,j)*cpice)) end do end do - - ! TODO(wjs, 2017-03-11) Include heat content of water in lakes, once we include - ! lake water as an explicit water state (https://github.com/ESCOMP/ctsm/issues/200) - + do fc = 1, num_lakec c = filter_lakec(fc) - heat(c) = heat_dry_mass(c) + heat_ice(c) + heat_liquid(c) + latent_heat_liquid(c) + heat(c) = heat(c) + heat_dry_mass(c) + heat_ice(c) + heat_liquid(c) + latent_heat_liquid(c) end do + + end associate - ! Subtract baselines set in initialization + end subroutine ComputeHeatLake + + !----------------------------------------------------------------------- + subroutine AccumulateHeatLake(bounds, num_c, filter_c, & + temperature_inst, lakestate_inst, & + heat) ! - ! NOTE(wjs, 2019-03-01) I haven't given enough thought to how (if at all) we should - ! correct for heat_liquid and cv_liquid, which are used to determine the weighted - ! average liquid water temperature. For example, if we're subtracting out a baseline - ! water amount because a particular water state is fictitious, we probably shouldn't - ! include that particular state when determining the weighted average temperature of - ! liquid water. And conversely, if we're adding a state via these baselines, should - ! we also add some water temperature of that state? The tricky thing here is what to - ! do when we end up subtracting water due to the baselines, so for now I'm simply not - ! trying to account for the temperature of these baselines at all. This amounts to - ! assuming that the baselines that we add / subtract are at the average temperature - ! of the real liquid water in the column. - do fc = 1, num_lakec - c = filter_lakec(fc) - heat(c) = heat(c) - dynbal_baseline_heat(c) + ! !DESCRIPTION: + ! Accumulate heat of lake water in lake columns. + ! + ! This subroutine differs from AccumulateSoilHeatNonLake in the sense that for lake heat + ! the average heat_liquid and cv_liquid are not accumulated. This is because these + ! terms are currently only used to calculate the implicit temperature of the dynbal liquid flux. + ! Because the lake water is virtual (there will never be a change in lake water content, + ! it should not be taken into the average column temperature when adjusting the change in + ! heat content of the grid cell for the change in water content. + ! Now, for lake grid cells, this is only done for the water content of the + ! soil under the lake and the snow on the lake. Since the virtual lake water doesn't generally + ! contribute to the dynbal liquid flux, its temperature shouldn't contribute + ! to the implicit temperature of the dynbal liquid flux. (If we allowed it + ! to contribute, the lake's temperature could dominate the average temperature calculation, + ! since there is so much lake water relative to other water in the grid cell.) + ! + ! We are adopting a different approach in the lake and non-lake columns. + ! For the choices made in a non-lake column, see comment at bottom of ComputeHeatNonLake subroutine + ! + ! Some minor concerns with this approach: + ! In some cases, lake water can have some changes in water content in time, + ! when experiencing phase changes: If a lake was completely liquid in initialization, + ! but then partially froze and then grew / shrank, some dynbal fluxes would be generated: + ! equal and opposite dynbal liquid and ice terms. In this case, it would be appropriate to + ! take the lake temperature along in determining the total heat which is corrected for delta liq. + + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_c ! number of column points in column filter (should not include lake points; can be a subset of the no-lake filter) + integer , intent(in) :: filter_c(:) ! column filter (should not include lake points; can be a subset of the no-lake filter) + type(temperature_type) , intent(in) :: temperature_inst + type(lakestate_type) , intent(in) :: lakestate_inst + real(r8) , intent(inout) :: heat( bounds%begc: ) ! accumulated heat content [J/m^2] + + ! !LOCAL VARIABLES: + integer :: fc + integer :: l, c, j + real(r8) :: h2olak_liq ! liquid water content of lake layer [kg/m²] + real(r8) :: h2olak_ice ! ice water content of lake layer [kg/m²] + real(r8) :: lake_heat_liquid(bounds%begc:bounds%endc) ! sum of heat content: liquid water in lake, excluding latent heat [J/m^2] + real(r8) :: lake_heat_ice(bounds%begc:bounds%endc) ! sum of heat content: ice in lake [J/m^2] + real(r8) :: lake_latent_heat_liquid(bounds%begc:bounds%endc) ! sum of heat content: latent heat of liquid water in lake [J/m^2] + + character(len=*), parameter :: subname = 'AccumulateHeatLake' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL_FL((ubound(heat) == (/bounds%endc/)), sourcefile, __LINE__) + + + associate( & + dz_lake => col%dz_lake, & ! lake layer depth (m) + t_lake => temperature_inst%t_lake_col, & ! lake temperature (K) + lake_icefrac => lakestate_inst%lake_icefrac_col & ! Input: [real(r8) (:,:) ] lake ice fraction + ) + + do fc = 1, num_c + c = filter_c(fc) + + lake_heat_liquid(c) = 0._r8 + lake_heat_ice(c) = 0._r8 + lake_latent_heat_liquid(c) = 0._r8 + end do + + + ! calculate heat content of lake itself + do j = 1, nlevlak + do fc = 1, num_c + c = filter_c(fc) + ! liquid heat + h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + call AccumulateLiquidWaterHeat( & + temp = t_lake(c,j), & + h2o = h2olak_liq, & + heat_liquid = lake_heat_liquid(c), & + latent_heat_liquid = lake_latent_heat_liquid(c)) + ! ice heat + ! use water density rather than ice density because lake layer depths are not + ! adjusted when the layer freezes + h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) + lake_heat_ice(c) = lake_heat_ice(c) + & + TempToHeat(temp=t_lake(c,j), cv = (h2olak_ice * cpice)) + end do + end do + + ! add ice heat and liquid heat together + do fc = 1, num_c + c = filter_c(fc) + heat(c) = heat(c) + (lake_heat_ice(c) + & + lake_heat_liquid(c) + lake_latent_heat_liquid(c)) end do end associate - end subroutine ComputeHeatLake - + end subroutine AccumulateHeatLake + !----------------------------------------------------------------------- subroutine AdjustDeltaHeatForDeltaLiq(bounds, delta_liq, & liquid_water_temp1, liquid_water_temp2, & @@ -955,10 +1168,16 @@ subroutine AdjustDeltaHeatForDeltaLiq(bounds, delta_liq, & ! ice runoff is at heat_base_temp (which is reasonable as long as heat_base_temp = ! tfrz). ! + ! With dynamical lakes, the adjusted delta_heat does not account for the added lake + ! water content due to growing lakes. This is because lake depth is constant, the + ! total lake water content (kg/m^2) does not change. The change in water content of + ! the snow and soil in the lake column are accounted for. + ! ! Eventually, if we begin to explicitly account for the temperature / heat content of ! liquid and ice runoff in CLM, then this routine should be reworked to use the true ! heat contents of both liquid and ice runoff. ! + ! ! Sign convention: delta_liq and delta_heat are positive if the post-landcover change ! value is greater than the pre-landcover change value. ! diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 018cd12176..089715233d 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -20,12 +20,13 @@ module dynConsBiogeophysMod use WaterFluxType , only : waterflux_type use WaterStateBulkType , only : waterstatebulk_type use WaterStateType , only : waterstate_type + use LakestateType , only : lakestate_type use WaterDiagnosticType , only : waterdiagnostic_type use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type use WaterBalanceType , only : waterbalance_type use WaterType , only : water_type - use TotalWaterAndHeatMod , only : AccumulateSoilLiqIceMassNonLake - use TotalWaterAndHeatMod , only : AccumulateSoilHeatNonLake + use TotalWaterAndHeatMod , only : AccumulateSoilLiqIceMassNonLake, AccumulateLiqIceMassLake + use TotalWaterAndHeatMod , only : AccumulateSoilHeatNonLake, AccumulateHeatLake use TotalWaterAndHeatMod , only : ComputeLiqIceMassNonLake, ComputeLiqIceMassLake use TotalWaterAndHeatMod , only : ComputeHeatNonLake, ComputeHeatLake use TotalWaterAndHeatMod , only : AdjustDeltaHeatForDeltaLiq @@ -64,7 +65,9 @@ module dynConsBiogeophysMod !----------------------------------------------------------------------- subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & - urbanparams_inst, soilstate_inst, water_inst, temperature_inst) + num_lakec, filter_lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, & + reset_all_baselines, reset_lake_baselines) ! ! !DESCRIPTION: ! Set start-of-run baseline values for heat and water content in some columns. @@ -77,13 +80,18 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & ! represented). These corrections will typically reduce the fictitious dynbal ! conservation fluxes. ! - ! At a minimum, this should be called during cold start initialization, to initialize - ! the baseline values based on cold start states. In addition, this can be called when - ! starting a new run from existing initial conditions, if the user desires. This - ! optional resetting can further reduce the dynbal fluxes; however, it can break - ! conservation. (So, for example, it can be done when transitioning from an offline - ! spinup to a coupled run, but it should not be done when transitioning from a - ! coupled historical run to a future scenario.) + ! At a minimum, this should be called during cold start initialization with + ! reset_all_baselines set to true, to initialize the baseline values based on cold + ! start states. This should also be called after reading initial conditions, but the + ! various reset_* flags should be set appropriately; setting all of these flags to + ! false will result in no baselines being reset in this call. + + ! Setting reset_all_baselines can be done when starting a new run from existing + ! initial conditions if the user desires. This optional resetting can further reduce + ! the dynbal fluxes; however, it can break conservation. (So, for example, it can be + ! done when transitioning from an offline spinup to a coupled run, but it should not + ! be done when transitioning from a coupled historical run to a future scenario.) + ! Other reset_* flags are described below. ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds @@ -93,11 +101,28 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & ! appropriate baseline value for that point. integer, intent(in) :: num_icemecc ! number of points in filter_icemecc integer, intent(in) :: filter_icemecc(:) ! filter for icemec (i.e., glacier) columns + integer, intent(in) :: num_lakec ! number of points in filter_lakec + integer, intent(in) :: filter_lakec(:) ! filter for lake columns type(urbanparams_type), intent(in) :: urbanparams_inst type(soilstate_type), intent(in) :: soilstate_inst + type(lakestate_type), intent(in) :: lakestate_inst type(water_type), intent(inout) :: water_inst type(temperature_type), intent(inout) :: temperature_inst + + ! Whether to reset baselines for all columns + logical, intent(in) :: reset_all_baselines + + ! Whether to reset baselines for lake columns. Ignored if reset_all_baselines is + ! true. + ! + ! BACKWARDS_COMPATIBILITY(wjs, 2020-09-02) This is needed when reading old initial + ! conditions files created before https://github.com/ESCOMP/CTSM/issues/1140 was + ! resolved via https://github.com/ESCOMP/CTSM/pull/1109: The definition of total + ! column water content has been changed for lakes, so we need to reset baseline values + ! for lakes on older initial conditions files. + logical, intent(in) :: reset_lake_baselines + ! ! !LOCAL VARIABLES: integer :: i @@ -118,22 +143,27 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & associate(bulk_or_tracer => water_inst%bulk_and_tracers(i)) call dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & - num_icemecc, filter_icemecc, & - bulk_or_tracer%waterstate_inst) + num_icemecc, filter_icemecc, num_lakec, filter_lakec, & + bulk_or_tracer%waterstate_inst, lakestate_inst, & + reset_all_baselines = reset_all_baselines, & + reset_lake_baselines = reset_lake_baselines) end associate end do - + call dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & - num_icemecc, filter_icemecc, & - urbanparams_inst, soilstate_inst, water_inst%waterstatebulk_inst, & - temperature_inst) + num_icemecc, filter_icemecc, num_lakec, filter_lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, water_inst%waterstatebulk_inst, & + temperature_inst, & + reset_all_baselines = reset_all_baselines, & + reset_lake_baselines = reset_lake_baselines) end subroutine dyn_hwcontent_set_baselines !----------------------------------------------------------------------- subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & - num_icemecc, filter_icemecc, & - waterstate_inst) + num_icemecc, filter_icemecc, num_lakec, filter_lakec, & + waterstate_inst, lakestate_inst, & + reset_all_baselines, reset_lake_baselines) ! ! !DESCRIPTION: ! Set start-of-run baseline values for water content, for a single water tracer or @@ -146,13 +176,23 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & ! The following filter should include inactive as well as active points integer, intent(in) :: num_icemecc ! number of points in filter_icemecc integer, intent(in) :: filter_icemecc(:) ! filter for icemec (i.e., glacier) columns + integer, intent(in) :: num_lakec ! number of points in filter_lakec + integer, intent(in) :: filter_lakec(:) ! filter for lake columns + type(lakestate_type), intent(in) :: lakestate_inst class(waterstate_type), intent(inout) :: waterstate_inst + + ! See dyn_hwcontent_set_baselines for documentation of these arguments + logical, intent(in) :: reset_all_baselines + logical, intent(in) :: reset_lake_baselines ! ! !LOCAL VARIABLES: + integer :: c, fc ! indices real(r8) :: soil_liquid_mass_col(bounds%begc:bounds%endc) real(r8) :: soil_ice_mass_col(bounds%begc:bounds%endc) - + real(r8) :: lake_liquid_mass_col(bounds%begc:bounds%endc) + real(r8) :: lake_ice_mass_col(bounds%begc:bounds%endc) + character(len=*), parameter :: subname = 'dyn_water_content_set_baselines' !----------------------------------------------------------------------- @@ -161,28 +201,49 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & dynbal_baseline_ice => waterstate_inst%dynbal_baseline_ice_col & ! Output: [real(r8) (:) ] baseline ice content subtracted from each column's total ice calculation (mm H2O) ) - soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8 - soil_ice_mass_col (bounds%begc:bounds%endc) = 0._r8 + if (reset_all_baselines) then + soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8 + soil_ice_mass_col (bounds%begc:bounds%endc) = 0._r8 + + call AccumulateSoilLiqIceMassNonLake(bounds, & + natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, & + waterstate_inst, & + liquid_mass = soil_liquid_mass_col(bounds%begc:bounds%endc), & + ice_mass = soil_ice_mass_col(bounds%begc:bounds%endc)) + + ! For glacier columns, the liquid and ice in the "soil" (i.e., in the glacial ice) is + ! a virtual pool. So we'll subtract this amount when determining the dynbal + ! fluxes. (Note that a positive value in these baseline variables indicates an extra + ! stock that will be subtracted later.) But glacier columns do not represent the soil + ! under the glacial ice. Let's assume that the soil state under each glacier column is + ! the same as the soil state in the natural vegetation landunit on that grid cell. We + ! subtract this from the dynbal baseline variables to indicate a missing stock. + call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, & + vals_col = soil_liquid_mass_col(bounds%begc:bounds%endc), & + baselines_col = dynbal_baseline_liq(bounds%begc:bounds%endc)) + call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, & + vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), & + baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc)) + end if - call AccumulateSoilLiqIceMassNonLake(bounds, & - natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, & - waterstate_inst, & - liquid_mass = soil_liquid_mass_col(bounds%begc:bounds%endc), & - ice_mass = soil_ice_mass_col(bounds%begc:bounds%endc)) - - ! For glacier columns, the liquid and ice in the "soil" (i.e., in the glacial ice) is - ! a virtual pool. So we'll subtract this amount when determining the dynbal - ! fluxes. (Note that a positive value in these baseline variables indicates an extra - ! stock that will be subtracted later.) But glacier columns do not represent the soil - ! under the glacial ice. Let's assume that the soil state under each glacier column is - ! the same as the soil state in the natural vegetation landunit on that grid cell. We - ! subtract this from the dynbal baseline variables to indicate a missing stock. - call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, & - vals_col = soil_liquid_mass_col(bounds%begc:bounds%endc), & - baselines_col = dynbal_baseline_liq(bounds%begc:bounds%endc)) - call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, & - vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), & - baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc)) + if (reset_all_baselines .or. reset_lake_baselines) then + ! set baselines for lake columns + + lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8 + lake_ice_mass_col (bounds%begc:bounds%endc) = 0._r8 + + ! Calculate the total water volume of the lake column + call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, & + tracer_ratio = waterstate_inst%info%get_ratio(), & + liquid_mass = lake_liquid_mass_col(bounds%begc:bounds%endc), & + ice_mass = lake_ice_mass_col(bounds%begc:bounds%endc)) + + do fc = 1, num_lakec + c = filter_lakec(fc) + dynbal_baseline_liq(c) = lake_liquid_mass_col(c) + dynbal_baseline_ice(c) = lake_ice_mass_col(c) + end do + end if end associate @@ -190,9 +251,10 @@ end subroutine dyn_water_content_set_baselines !----------------------------------------------------------------------- subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & - num_icemecc, filter_icemecc, & - urbanparams_inst, soilstate_inst, waterstatebulk_inst, & - temperature_inst) + num_icemecc, filter_icemecc, num_lakec, filter_lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, waterstatebulk_inst, & + temperature_inst, & + reset_all_baselines, reset_lake_baselines) ! ! !DESCRIPTION: ! Set start-of-run baseline values for heat content. @@ -202,19 +264,28 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & type(filter_col_type), intent(in) :: natveg_and_glc_filterc ! filter for natural veg and glacier columns ! The following filter should include inactive as well as active points - integer, intent(in) :: num_icemecc ! number of points in filter_icemecc + integer, intent(in) :: num_icemecc ! number of points in filter_icemecc integer, intent(in) :: filter_icemecc(:) ! filter for icemec (i.e., glacier) columns + integer, intent(in) :: num_lakec ! number of points in filter_lakec + integer, intent(in) :: filter_lakec(:) ! filter for lake columns - type(urbanparams_type), intent(in) :: urbanparams_inst - type(soilstate_type), intent(in) :: soilstate_inst + type(urbanparams_type), intent(in) :: urbanparams_inst + type(soilstate_type) , intent(in) :: soilstate_inst + type(lakestate_type) , intent(in) :: lakestate_inst type(waterstatebulk_type), intent(in) :: waterstatebulk_inst type(temperature_type), intent(inout) :: temperature_inst + + ! See dyn_hwcontent_set_baselines for documentation of these arguments + logical, intent(in) :: reset_all_baselines + logical, intent(in) :: reset_lake_baselines ! ! !LOCAL VARIABLES: - real(r8) :: soil_heat_col(bounds%begc:bounds%endc) ! soil heat content [J/m^2] - real(r8) :: soil_heat_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface - real(r8) :: soil_cv_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface - + integer :: c, fc ! indices + real(r8) :: soil_heat_col(bounds%begc:bounds%endc) ! soil heat content [J/m^2] + real(r8) :: soil_heat_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface + real(r8) :: soil_cv_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface + real(r8) :: lake_heat_col(bounds%begc:bounds%endc) ! lake heat content [J/m^2] + character(len=*), parameter :: subname = 'dyn_heat_content_set_baselines' !----------------------------------------------------------------------- @@ -222,34 +293,51 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & dynbal_baseline_heat => temperature_inst%dynbal_baseline_heat_col & ! Output: [real(r8) (:) ] baseline heat content subtracted from each column's total heat calculation (J/m2) ) - soil_heat_col(bounds%begc:bounds%endc) = 0._r8 - soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8 - soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8 - - call AccumulateSoilHeatNonLake(bounds, & - natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, & - urbanparams_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, & - heat = soil_heat_col(bounds%begc:bounds%endc), & - heat_liquid = soil_heat_liquid_col(bounds%begc:bounds%endc), & - cv_liquid = soil_cv_liquid_col(bounds%begc:bounds%endc)) - - ! See comments in dyn_water_content_set_baselines for rationale for these glacier - ! baselines. Even though the heat in glacier ice can interact with the rest of the - ! system (e.g., giving off heat to the atmosphere or receiving heat from the - ! atmosphere), it is still a virtual state in the sense that the glacier ice column - ! is virtual. And, as for water, we subtract the heat of the soil in the associated - ! natural vegetated landunit to account for the fact that we don't explicitly model - ! the soil under glacial ice. - ! - ! Aside from these considerations of what is virtual vs. real, another rationale - ! driving the use of these baselines is the desire to minimize the dynbal fluxes. For - ! the sake of conservation, it seems okay to pick any fixed baseline when summing the - ! energy (or water) content of a given column (as long as that baseline doesn't - ! change over time). By using the baselines computed here, we reduce the dynbal - ! fluxes to more reasonable values. - call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, & - vals_col = soil_heat_col(bounds%begc:bounds%endc), & - baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc)) + if (reset_all_baselines) then + soil_heat_col(bounds%begc:bounds%endc) = 0._r8 + soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8 + soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8 + + call AccumulateSoilHeatNonLake(bounds, & + natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, & + urbanparams_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, & + heat = soil_heat_col(bounds%begc:bounds%endc), & + heat_liquid = soil_heat_liquid_col(bounds%begc:bounds%endc), & + cv_liquid = soil_cv_liquid_col(bounds%begc:bounds%endc)) + + ! See comments in dyn_water_content_set_baselines for rationale for these glacier + ! baselines. Even though the heat in glacier ice can interact with the rest of the + ! system (e.g., giving off heat to the atmosphere or receiving heat from the + ! atmosphere), it is still a virtual state in the sense that the glacier ice column + ! is virtual. And, as for water, we subtract the heat of the soil in the associated + ! natural vegetated landunit to account for the fact that we don't explicitly model + ! the soil under glacial ice. + ! + ! Aside from these considerations of what is virtual vs. real, another rationale + ! driving the use of these baselines is the desire to minimize the dynbal fluxes. For + ! the sake of conservation, it seems okay to pick any fixed baseline when summing the + ! energy (or water) content of a given column (as long as that baseline doesn't + ! change over time). By using the baselines computed here, we reduce the dynbal + ! fluxes to more reasonable values. + call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, & + vals_col = soil_heat_col(bounds%begc:bounds%endc), & + baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc)) + end if + + if (reset_all_baselines .or. reset_lake_baselines) then + ! Set baselines for lake columns + + lake_heat_col(bounds%begc:bounds%endc) = 0._r8 + + call AccumulateHeatLake(bounds, num_lakec, filter_lakec, & + temperature_inst, lakestate_inst, & + heat = lake_heat_col) + + do fc = 1, num_lakec + c = filter_lakec(fc) + dynbal_baseline_heat(c) = lake_heat_col(c) + end do + end if end associate @@ -334,7 +422,7 @@ subroutine dyn_hwcontent_init(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & - water_inst, temperature_inst) + water_inst, temperature_inst, lakestate_inst) ! ! !DESCRIPTION: ! Compute grid cell-level heat and water content before land cover change @@ -350,7 +438,9 @@ subroutine dyn_hwcontent_init(bounds, & type(urbanparams_type) , intent(in) :: urbanparams_inst type(soilstate_type) , intent(in) :: soilstate_inst type(water_type) , intent(inout) :: water_inst + type(lakestate_type) , intent(in) :: lakestate_inst type(temperature_type) , intent(inout) :: temperature_inst + ! ! !LOCAL VARIABLES: integer :: i @@ -369,6 +459,7 @@ subroutine dyn_hwcontent_init(bounds, & num_lakec, filter_lakec, & bulk_or_tracer%waterstate_inst, & bulk_or_tracer%waterdiagnostic_inst, & + lakestate_inst, & liquid_mass = bulk_or_tracer%waterbalance_inst%liq1_grc(begg:endg), & ice_mass = bulk_or_tracer%waterbalance_inst%ice1_grc(begg:endg)) end associate @@ -380,6 +471,7 @@ subroutine dyn_hwcontent_init(bounds, & urbanparams_inst, soilstate_inst, & temperature_inst, & water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & + lakestate_inst, & heat_grc = temperature_inst%heat1_grc(begg:endg), & liquid_water_temp_grc = temperature_inst%liquid_water_temp1_grc(begg:endg)) @@ -393,7 +485,7 @@ subroutine dyn_hwcontent_final(bounds, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & water_inst, & - temperature_inst, energyflux_inst) + temperature_inst, energyflux_inst, lakestate_inst) ! ! !DESCRIPTION: ! Compute grid cell-level heat and water content and dynbal fluxes after land cover change @@ -409,6 +501,7 @@ subroutine dyn_hwcontent_final(bounds, & type(urbanparams_type) , intent(in) :: urbanparams_inst type(soilstate_type) , intent(in) :: soilstate_inst type(water_type) , intent(inout) :: water_inst + type(lakestate_type) , intent(in) :: lakestate_inst type(temperature_type) , intent(inout) :: temperature_inst type(energyflux_type) , intent(inout) :: energyflux_inst ! @@ -433,6 +526,7 @@ subroutine dyn_hwcontent_final(bounds, & bulk_or_tracer%waterdiagnostic_inst, & bulk_or_tracer%waterbalance_inst, & bulk_or_tracer%waterflux_inst, & + lakestate_inst, & delta_liq = this_delta_liq(begg:endg)) if (i == water_inst%i_bulk) then delta_liq_bulk(begg:endg) = this_delta_liq(begg:endg) @@ -446,6 +540,7 @@ subroutine dyn_hwcontent_final(bounds, & urbanparams_inst, soilstate_inst, & temperature_inst, & water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & + lakestate_inst, & heat_grc = temperature_inst%heat2_grc(begg:endg), & liquid_water_temp_grc = temperature_inst%liquid_water_temp2_grc(begg:endg)) @@ -481,7 +576,7 @@ subroutine dyn_water_content_final(bounds, & num_lakec, filter_lakec, & waterstate_inst, waterdiagnostic_inst, & waterbalance_inst, waterflux_inst, & - delta_liq) + lakestate_inst, delta_liq) ! ! !DESCRIPTION: ! Compute grid cell-level water content and dynbal fluxes after landcover change, for @@ -497,6 +592,7 @@ subroutine dyn_water_content_final(bounds, & class(waterdiagnostic_type) , intent(in) :: waterdiagnostic_inst class(waterbalance_type) , intent(inout) :: waterbalance_inst class(waterflux_type) , intent(inout) :: waterflux_inst + type(lakestate_type) , intent(in) :: lakestate_inst real(r8) , intent(out) :: delta_liq(bounds%begg:) ! change in gridcell h2o liq content ! ! !LOCAL VARIABLES: @@ -515,7 +611,7 @@ subroutine dyn_water_content_final(bounds, & call dyn_water_content(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & - waterstate_inst, waterdiagnostic_inst, & + waterstate_inst, waterdiagnostic_inst, lakestate_inst, & liquid_mass = waterbalance_inst%liq2_grc(bounds%begg:bounds%endg), & ice_mass = waterbalance_inst%ice2_grc(bounds%begg:bounds%endg)) @@ -549,7 +645,7 @@ end subroutine dyn_water_content_final subroutine dyn_water_content(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & - waterstate_inst, waterdiagnostic_inst, & + waterstate_inst, waterdiagnostic_inst, lakestate_inst, & liquid_mass, ice_mass) ! ! !DESCRIPTION: @@ -563,6 +659,7 @@ subroutine dyn_water_content(bounds, & integer , intent(in) :: filter_lakec(:) class(waterstate_type) , intent(in) :: waterstate_inst class(waterdiagnostic_type) , intent(in) :: waterdiagnostic_inst + type(lakestate_type) , intent(in) :: lakestate_inst real(r8) , intent(out) :: liquid_mass( bounds%begg: ) ! kg m-2 real(r8) , intent(out) :: ice_mass( bounds%begg: ) ! kg m-2 ! @@ -584,7 +681,8 @@ subroutine dyn_water_content(bounds, & call ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & waterstate_inst, & - subtract_dynbal_baselines = .true., & + lakestate_inst, & + add_lake_water_and_subtract_dynbal_baselines = .true., & liquid_mass = liquid_mass_col(bounds%begc:bounds%endc), & ice_mass = ice_mass_col(bounds%begc:bounds%endc)) @@ -608,7 +706,7 @@ subroutine dyn_heat_content(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & - temperature_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, & + temperature_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, lakestate_inst, & heat_grc, liquid_water_temp_grc) ! !DESCRIPTION: @@ -631,6 +729,7 @@ subroutine dyn_heat_content(bounds, & type(temperature_type) , intent(in) :: temperature_inst type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst + type(lakestate_type) , intent(in) :: lakestate_inst real(r8) , intent(out) :: heat_grc( bounds%begg: ) ! total heat content for each grid cell [J/m^2] real(r8) , intent(out) :: liquid_water_temp_grc( bounds%begg: ) ! weighted average liquid water temperature for each grid cell (K) @@ -659,7 +758,7 @@ subroutine dyn_heat_content(bounds, & cv_liquid = cv_liquid_col(bounds%begc:bounds%endc)) call ComputeHeatLake(bounds, num_lakec, filter_lakec, & - soilstate_inst, temperature_inst, waterstatebulk_inst, & + soilstate_inst, temperature_inst, waterstatebulk_inst, lakestate_inst, & heat = heat_col(bounds%begc:bounds%endc), & heat_liquid = heat_liquid_col(bounds%begc:bounds%endc), & cv_liquid = cv_liquid_col(bounds%begc:bounds%endc)) diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index 3820d6392a..6d9a703b31 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -24,6 +24,7 @@ module dynSubgridControlMod public :: get_flanduse_timeseries ! return the value of the flanduse_timeseries file name public :: get_do_transient_pfts ! return the value of the do_transient_pfts control flag public :: get_do_transient_crops ! return the value of the do_transient_crops control flag + public :: get_do_transient_lakes ! return the value of the do_transient_lakes control flag public :: run_has_transient_landcover ! returns true if any aspects of prescribed transient landcover are enabled public :: get_do_harvest ! return the value of the do_harvest control flag public :: get_reset_dynbal_baselines ! return the value of the reset_dynbal_baselines control flag @@ -40,6 +41,7 @@ module dynSubgridControlMod character(len=fname_len) :: flanduse_timeseries = ' ' ! transient landuse dataset logical :: do_transient_pfts = .false. ! whether to apply transient natural PFTs from dataset logical :: do_transient_crops = .false. ! whether to apply transient crops from dataset + logical :: do_transient_lakes = .false. ! whether to apply transient lakes from dataset logical :: do_harvest = .false. ! whether to apply harvest from dataset logical :: reset_dynbal_baselines = .false. ! whether to reset baseline values of total column water and energy in the first step of the run @@ -116,6 +118,7 @@ subroutine read_namelist( NLFilename ) character(len=fname_len) :: flanduse_timeseries logical :: do_transient_pfts logical :: do_transient_crops + logical :: do_transient_lakes logical :: do_harvest logical :: reset_dynbal_baselines logical :: for_testing_allow_non_annual_changes @@ -131,6 +134,7 @@ subroutine read_namelist( NLFilename ) flanduse_timeseries, & do_transient_pfts, & do_transient_crops, & + do_transient_lakes, & do_harvest, & reset_dynbal_baselines, & for_testing_allow_non_annual_changes, & @@ -140,6 +144,7 @@ subroutine read_namelist( NLFilename ) flanduse_timeseries = ' ' do_transient_pfts = .false. do_transient_crops = .false. + do_transient_lakes = .false. do_harvest = .false. reset_dynbal_baselines = .false. for_testing_allow_non_annual_changes = .false. @@ -164,6 +169,7 @@ subroutine read_namelist( NLFilename ) call shr_mpi_bcast (flanduse_timeseries, mpicom) call shr_mpi_bcast (do_transient_pfts, mpicom) call shr_mpi_bcast (do_transient_crops, mpicom) + call shr_mpi_bcast (do_transient_lakes, mpicom) call shr_mpi_bcast (do_harvest, mpicom) call shr_mpi_bcast (reset_dynbal_baselines, mpicom) call shr_mpi_bcast (for_testing_allow_non_annual_changes, mpicom) @@ -173,6 +179,7 @@ subroutine read_namelist( NLFilename ) flanduse_timeseries = flanduse_timeseries, & do_transient_pfts = do_transient_pfts, & do_transient_crops = do_transient_crops, & + do_transient_lakes = do_transient_lakes, & do_harvest = do_harvest, & reset_dynbal_baselines = reset_dynbal_baselines, & for_testing_allow_non_annual_changes = for_testing_allow_non_annual_changes, & @@ -218,6 +225,11 @@ subroutine check_namelist_consistency write(iulog,*) 'a flanduse_timeseries file (currently flanduse_timeseries is blank)' call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if (dyn_subgrid_control_inst%do_transient_lakes) then + write(iulog,*) 'ERROR: do_transient_lakes can only be true if you are running with' + write(iulog,*) 'a flanduse_timeseries file (currently flanduse_timeseries is blank)' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if if (dyn_subgrid_control_inst%do_harvest) then write(iulog,*) 'ERROR: do_harvest can only be true if you are running with' write(iulog,*) 'a flanduse_timeseries file (currently flanduse_timeseries is blank)' @@ -236,9 +248,18 @@ subroutine check_namelist_consistency end if end if - if (dyn_subgrid_control_inst%do_transient_pfts .or. dyn_subgrid_control_inst%do_transient_crops) then + ! NOTE(wjs, 2020-08-23) In the following error checks, I'm treating do_transient_lakes + ! similar to do_transient_pfts and do_transient_crops. I'm not sure if all of these + ! checks are truly important for transient lakes (in particular, my guess is that + ! collapse_urban could probably be done with transient lakes - as well as transient + ! pfts and transient crops for that matter), but some of the checks probably are + ! needed, and it seems best to keep transient lakes consistent with other transient + ! areas in this respect. + if (dyn_subgrid_control_inst%do_transient_pfts .or. & + dyn_subgrid_control_inst%do_transient_crops .or. & + dyn_subgrid_control_inst%do_transient_lakes) then if (collapse_urban) then - write(iulog,*) 'ERROR: do_transient_pfts and do_transient_crops are & + write(iulog,*) 'ERROR: do_transient_pfts, do_transient_crops and do_transient_lakes are & incompatible with collapse_urban = .true.' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -246,7 +267,7 @@ subroutine check_namelist_consistency .or. toosmall_soil > 0._r8 .or. toosmall_crop > 0._r8 & .or. toosmall_glacier > 0._r8 .or. toosmall_lake > 0._r8 & .or. toosmall_wetland > 0._r8 .or. toosmall_urban > 0._r8) then - write(iulog,*) 'ERROR: do_transient_pfts and do_transient_crops are & + write(iulog,*) 'ERROR: do_transient_pfts, do_transient_crops and do_transient_lakes are & incompatible with any of the following set to > 0: & n_dom_pfts > 0, n_dom_landunits > 0, & toosmall_soil > 0._r8, toosmall_crop > 0._r8, & @@ -275,8 +296,8 @@ subroutine check_namelist_consistency write(iulog,*) 'ERROR: do_harvest currently does not work with use_fates' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if - + end if + end subroutine check_namelist_consistency !----------------------------------------------------------------------- @@ -316,6 +337,18 @@ logical function get_do_transient_crops() get_do_transient_crops = dyn_subgrid_control_inst%do_transient_crops end function get_do_transient_crops + + !----------------------------------------------------------------------- + logical function get_do_transient_lakes() + ! !DESCRIPTION: + ! Return the value of the do_transient_lakes control flag + !----------------------------------------------------------------------- + + SHR_ASSERT(dyn_subgrid_control_inst%initialized, errMsg(sourcefile, __LINE__)) + + get_do_transient_lakes = dyn_subgrid_control_inst%do_transient_lakes + + end function get_do_transient_lakes !----------------------------------------------------------------------- logical function run_has_transient_landcover() diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index 73c3f9e4b0..676e14da81 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -12,7 +12,7 @@ module dynSubgridDriverMod use decompMod , only : bounds_type, BOUNDS_LEVEL_PROC, BOUNDS_LEVEL_CLUMP use decompMod , only : get_proc_clumps, get_clump_bounds use dynSubgridControlMod , only : get_flanduse_timeseries - use dynSubgridControlMod , only : get_do_transient_pfts, get_do_transient_crops + use dynSubgridControlMod , only : get_do_transient_pfts, get_do_transient_crops, get_do_transient_lakes use dynSubgridControlMod , only : get_do_harvest use dynPriorWeightsMod , only : prior_weights_type use dynPatchStateUpdaterMod , only : patch_state_updater_type @@ -20,6 +20,7 @@ module dynSubgridDriverMod use dynpftFileMod , only : dynpft_init, dynpft_interp use dyncropFileMod , only : dyncrop_init, dyncrop_interp use dynHarvestMod , only : dynHarvest_init, dynHarvest_interp + use dynlakeFileMod , only : dynlake_init, dynlake_interp use dynLandunitAreaMod , only : update_landunit_weights use subgridWeightsMod , only : compute_higher_order_weights, set_subgrid_diagnostic_fields use reweightMod , only : reweight_wrapup @@ -37,6 +38,7 @@ module dynSubgridDriverMod use SoilHydrologyType , only : soilhydrology_type use SoilStateType , only : soilstate_type use WaterType , only : water_type + use LakestateType , only : lakestate_type use TemperatureType , only : temperature_type use CropType , only : crop_type use glc2lndMod , only : glc2lnd_type @@ -120,6 +122,12 @@ subroutine dynSubgrid_init(bounds_proc, glc_behavior, crop_inst) call dynHarvest_init(bounds_proc, harvest_filename=get_flanduse_timeseries()) end if + + ! Initialize stuff for prescribed transient lakes + if (get_do_transient_lakes()) then + call dynlake_init(bounds_proc, dynlake_filename=get_flanduse_timeseries()) + end if + ! ------------------------------------------------------------------------ ! Set initial subgrid weights for aspects that are read from file. This is relevant ! for cold start and use_init_interp-based initialization. @@ -133,6 +141,11 @@ subroutine dynSubgrid_init(bounds_proc, glc_behavior, crop_inst) call dyncrop_interp(bounds_proc, crop_inst) end if + if (get_do_transient_lakes()) then + call dynlake_interp(bounds_proc) + end if + + ! (We don't bother calling dynHarvest_interp, because the harvest information isn't ! needed until the run loop. Harvest has nothing to do with subgrid weights, and in ! some respects doesn't even really belong in this module at all.) @@ -153,7 +166,7 @@ end subroutine dynSubgrid_init !----------------------------------------------------------------------- subroutine dynSubgrid_driver(bounds_proc, & urbanparams_inst, soilstate_inst, water_inst, & - temperature_inst, energyflux_inst, & + temperature_inst, energyflux_inst, lakestate_inst, & canopystate_inst, photosyns_inst, crop_inst, glc2lnd_inst, bgc_vegetation_inst, & soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonstate_inst, & @@ -181,6 +194,7 @@ subroutine dynSubgrid_driver(bounds_proc, type(urbanparams_type) , intent(in) :: urbanparams_inst type(soilstate_type) , intent(in) :: soilstate_inst type(water_type) , intent(inout) :: water_inst + type(lakestate_type) , intent(in) :: lakestate_inst type(temperature_type) , intent(inout) :: temperature_inst type(energyflux_type) , intent(inout) :: energyflux_inst type(canopystate_type) , intent(inout) :: canopystate_inst @@ -221,7 +235,7 @@ subroutine dynSubgrid_driver(bounds_proc, filter(nc)%num_nolakec, filter(nc)%nolakec, & filter(nc)%num_lakec, filter(nc)%lakec, & urbanparams_inst, soilstate_inst, & - water_inst, temperature_inst) + water_inst, temperature_inst, lakestate_inst) call prior_weights%set_prior_weights(bounds_clump) call patch_state_updater%set_old_weights(bounds_clump) @@ -244,7 +258,10 @@ subroutine dynSubgrid_driver(bounds_proc, if (get_do_harvest()) then call dynHarvest_interp(bounds_proc) end if - + + if (get_do_transient_lakes()) then + call dynlake_interp(bounds_proc) + end if ! ========================================================================== ! Do land cover change that does not require I/O ! ========================================================================== @@ -292,7 +309,7 @@ subroutine dynSubgrid_driver(bounds_proc, filter(nc)%num_lakec, filter(nc)%lakec, & urbanparams_inst, soilstate_inst, & water_inst, & - temperature_inst, energyflux_inst) + temperature_inst, energyflux_inst, lakestate_inst) if (use_cn) then call bgc_vegetation_inst%DynamicAreaConservation(bounds_clump, nc, & diff --git a/src/dyn_subgrid/dynlakeFileMod.F90 b/src/dyn_subgrid/dynlakeFileMod.F90 new file mode 100644 index 0000000000..17b21484fd --- /dev/null +++ b/src/dyn_subgrid/dynlakeFileMod.F90 @@ -0,0 +1,133 @@ +module dynlakeFileMod + + !--------------------------------------------------------------------------- + ! !DESCRIPTION: + ! Handle reading of the dataset that specifies transient areas of the lake landunit + ! + ! !USES: + +#include "shr_assert.h" + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_kind_mod , only : r8 => shr_kind_r8 + use decompMod , only : bounds_type, BOUNDS_LEVEL_PROC + use dynFileMod , only : dyn_file_type + use dynVarTimeUninterpMod , only : dyn_var_time_uninterp_type + use clm_varctl , only : iulog + use clm_varcon , only : grlnd + use abortutils , only : endrun + use spmdMod , only : masterproc + + ! !PUBLIC MEMBER FUNCTIONS: + implicit none + private + save + public :: dynlake_init ! initialize information read from landuse.timeseries dataset + public :: dynlake_interp ! get landuse data for the current time step, if needed + ! + ! ! PRIVATE TYPES + type(dyn_file_type), target :: dynlake_file ! information for the file containing transient lake data + type(dyn_var_time_uninterp_type) :: wtlake ! weight of the lake landunit + + ! Names of variables on file + character(len=*), parameter :: lake_varname = 'PCT_LAKE' + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + !--------------------------------------------------------------------------- + +contains + + !----------------------------------------------------------------------- + subroutine dynlake_init(bounds, dynlake_filename) + ! + ! !DESCRIPTION: + ! Initialize dataset containing transient lake info (position it to the right time + ! samples that bound the initial model date) + ! + ! !USES: + use ncdio_pio , only : check_dim + use dynTimeInfoMod , only : YEAR_POSITION_START_OF_TIMESTEP + + + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds ! proc-level bounds + character(len=*) , intent(in) :: dynlake_filename ! name of file containing transient lake information + ! + ! !LOCAL VARIABLES: + integer :: num_points ! number of spatial points + + character(len=*), parameter :: subname = 'dynlake_init' + !----------------------------------------------------------------------- + + SHR_ASSERT(bounds%level == BOUNDS_LEVEL_PROC, subname // ': argument must be PROC-level bounds') + + if (masterproc) then + write(iulog,*) 'Attempting to read lake dynamic landuse data .....' + end if + + ! Get the year from the START of the timestep; this way, we'll update lake areas + ! starting after the year boundary. This is consistent with the timing of glacier + ! updates, and will likely be consistent with the timing of lake updates determined + ! prognostically, if lake areas are ever determined prognostically rather than + ! prescribed ahead of time. + dynlake_file = dyn_file_type(dynlake_filename, YEAR_POSITION_START_OF_TIMESTEP) + + ! read data PCT_LAKE + ! + ! Note: if you want to change transient lakes so that they are interpolated, rather + ! than jumping to each year's value on Jan 1 of that year, simply change wtlake and + ! to be of type dyn_var_time_interp_type (rather than + ! dyn_var_time_uninterp_type), and change the following constructors to construct + ! variables of dyn_var_time_interp_type. That's all you need to do. + num_points = (bounds%endg - bounds%begg + 1) + wtlake = dyn_var_time_uninterp_type( & + dyn_file = dynlake_file, varname=lake_varname, & + dim1name=grlnd, conversion_factor=100._r8, & + do_check_sums_equal_1=.false., data_shape=[num_points]) + + end subroutine dynlake_init + + + !----------------------------------------------------------------------- + subroutine dynlake_interp(bounds) + ! + ! !DESCRIPTION: + ! Get lake cover for model time, when needed. + ! + ! Sets col%wtlunit and lun%wtgcell for lake landunits. + ! + ! Note that lake cover currently jumps to its new value at the start of the year. + ! However, as mentioned above, this behavior can be changed to time interpolation + ! simply by making wtlake and wtcft dyn_var_time_interp_type variables rather than + ! dyn_var_time_uninterp_type. + ! + ! !USES: + use landunit_varcon , only : istdlak + use subgridWeightsMod , only : set_landunit_weight + ! + ! !ARGUMENTS: + type(bounds_type), intent(in) :: bounds ! proc-level bounds + ! + ! !LOCAL VARIABLES: + integer :: g ! indices + real(r8), allocatable :: wtlake_cur(:) ! current weight of the lake landunit + + character(len=*), parameter :: subname = 'dynlake_interp' + !----------------------------------------------------------------------- + + SHR_ASSERT(bounds%level == BOUNDS_LEVEL_PROC, subname // ': argument must be PROC-level bounds') + + call dynlake_file%time_info%set_current_year() + + ! Set new landunit area + allocate(wtlake_cur(bounds%begg:bounds%endg)) + call wtlake%get_current_data(wtlake_cur) + do g = bounds%begg, bounds%endg + call set_landunit_weight(g, istdlak, wtlake_cur(g)) + end do + deallocate(wtlake_cur) + + end subroutine dynlake_interp + +end module dynlakeFileMod diff --git a/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf b/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf index 54f56b333f..47f802e5bf 100644 --- a/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf +++ b/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf @@ -7,7 +7,7 @@ module test_dyn_cons_biogeophys use shr_kind_mod , only : r8 => shr_kind_r8 use unittestSubgridMod use unittestArrayMod, only : col_array - use unittestFilterBuilderMod, only : filter_from_range + use unittestFilterBuilderMod, only : filter_from_range, filter_empty use unittestWaterTypeFactory, only : unittest_water_type_factory_type use clm_varpar, only : nlevgrnd, nlevsno, maxpatch_glcmec use column_varcon, only : icemec_class_to_col_itype @@ -18,6 +18,7 @@ module test_dyn_cons_biogeophys use SoilStateType, only : soilstate_type use TemperatureType, only : temperature_type use UrbanParamsType, only : urbanparams_type + use LakestateType, only : lakestate_type use WaterType, only : water_type use TotalWaterAndHeatMod, only : AccumulateSoilLiqIceMassNonLake use TotalWaterAndHeatMod, only : AccumulateSoilHeatNonLake @@ -30,6 +31,7 @@ module test_dyn_cons_biogeophys type(soilstate_type) :: soilstate_inst type(temperature_type) :: temperature_inst type(urbanparams_type) :: urbanparams_inst + type(lakestate_type) :: lakestate_inst type(water_type) :: water_inst contains procedure :: setUp @@ -99,8 +101,8 @@ contains real(r8), parameter :: wt_veg_col1 = 0.25_r8 real(r8), parameter :: wt_veg_col2 = 0.75_r8 integer :: veg_col1, veg_col2, glc_col - integer :: num_icemecc, num_natveg_and_icemecc - integer, allocatable :: filter_icemecc(:), filter_natveg_and_icemecc(:) + integer :: num_icemecc, num_natveg_and_icemecc, num_lakec + integer, allocatable :: filter_icemecc(:), filter_natveg_and_icemecc(:), filter_lakec(:) integer :: c real(r8), allocatable :: expected_vals_liq_col(:) real(r8), allocatable :: expected_vals_ice_col(:) @@ -160,6 +162,9 @@ contains numf=num_icemecc, filter=filter_icemecc) call filter_from_range(start=veg_col1, end=glc_col, & numf=num_natveg_and_icemecc, filter=filter_natveg_and_icemecc) + ! For now, this test does NOT cover lake columns. So just use an empty lake filter. + ! This also allows us to avoid initializing variables in this%lakestate_inst. + call filter_empty(num_lakec, filter_lakec) ! Initialize necessary variables in soilstate_inst allocate(this%soilstate_inst%csol_col(bounds%begc:bounds%endc, nlevgrnd)) @@ -203,8 +208,10 @@ contains ! ------------------------------------------------------------------------ call dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & - this%urbanparams_inst, this%soilstate_inst, & - this%water_inst, this%temperature_inst) + num_lakec, filter_lakec, & + this%urbanparams_inst, this%soilstate_inst, this%lakestate_inst, & + this%water_inst, this%temperature_inst, & + reset_all_baselines = .true., reset_lake_baselines = .false.) ! ------------------------------------------------------------------------ ! Compute expected values diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index 51647c105e..7e802e0936 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -338,7 +338,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call t_startf('dyn_subgrid') call dynSubgrid_driver(bounds_proc, & urbanparams_inst, soilstate_inst, water_inst, & - temperature_inst, energyflux_inst, & + temperature_inst, energyflux_inst, lakestate_inst, & canopystate_inst, photosyns_inst, crop_inst, glc2lnd_inst, bgc_vegetation_inst, & soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonstate_inst, & @@ -384,7 +384,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call BeginWaterBalance(bounds_clump, & filter(nc)%num_nolakec, filter(nc)%nolakec, & filter(nc)%num_lakec, filter(nc)%lakec, & - water_inst, soilhydrology_inst, & + water_inst, soilhydrology_inst, lakestate_inst, & use_aquifer_layer = use_aquifer_layer()) call t_stopf('begwbal') diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 4616aecc09..912bbb293b 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -14,7 +14,7 @@ module clm_initializeMod use clm_varctl , only : iulog use clm_varctl , only : use_lch4, use_cn, use_cndv, use_c13, use_c14, use_fates use clm_varctl , only : use_soil_moisture_streams - use clm_instur , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, fert_cft, irrig_method, wt_glc_mec, topo_glc_mec + use clm_instur , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, fert_cft, irrig_method, wt_glc_mec, topo_glc_mec, haslake use perf_mod , only : t_startf, t_stopf use readParamsMod , only : readParameters use ncdio_pio , only : file_desc_t @@ -183,7 +183,7 @@ subroutine initialize1(dtime, gindex_ocn) allocate (irrig_method (begg:endg, cft_lb:cft_ub )) allocate (wt_glc_mec (begg:endg, maxpatch_glcmec)) allocate (topo_glc_mec(begg:endg, maxpatch_glcmec)) - + allocate (haslake (begg:endg )) ! Read list of Patches and their corresponding parameter values ! Independent of model resolution, Needs to stay before surfrd_get_data @@ -265,7 +265,7 @@ subroutine initialize1(dtime, gindex_ocn) ! Some things are kept until the end of initialize2; urban_valid is kept through the ! end of the run for error checking. - deallocate (wt_lunit, wt_cft, wt_glc_mec) + deallocate (wt_lunit, wt_cft, wt_glc_mec, haslake) call t_stopf('clm_init1') @@ -278,6 +278,7 @@ subroutine initialize2( ) ! CLM initialization - second phase ! ! !USES: + use shr_orb_mod , only : shr_orb_decl use shr_scam_mod , only : shr_scam_getCloseLatLon use seq_drydep_mod , only : n_drydep, drydep_method, DD_XLND @@ -338,6 +339,8 @@ subroutine initialize2( ) logical :: lexist integer :: closelatidx,closelonidx real(r8) :: closelat,closelon + logical :: reset_dynbal_baselines_all_columns + logical :: reset_dynbal_baselines_lake_columns integer :: begp, endp integer :: begc, endc integer :: begl, endl @@ -487,8 +490,14 @@ subroutine initialize2( ) call dyn_hwcontent_set_baselines(bounds_clump, & filter_inactive_and_active(nc)%num_icemecc, & filter_inactive_and_active(nc)%icemecc, & - urbanparams_inst, soilstate_inst, water_inst, temperature_inst) + filter_inactive_and_active(nc)%num_lakec, & + filter_inactive_and_active(nc)%lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, & + reset_all_baselines = .true., & + ! reset_lake_baselines is irrelevant since reset_all_baselines is true + reset_lake_baselines = .false.) end do + !$OMP END PARALLEL DO ! ------------------------------------------------------------------------ ! Initialize modules (after time-manager initialization in most cases) @@ -540,6 +549,7 @@ subroutine initialize2( ) is_cold_start = .false. is_interpolated_start = .false. + reset_dynbal_baselines_lake_columns = .false. if (nsrest == nsrStartup) then @@ -560,7 +570,8 @@ subroutine initialize2( ) write(iulog,*)'Reading initial conditions from ',trim(finidat) end if call getfil( finidat, fnamer, 0 ) - call restFile_read(bounds_proc, fnamer, glc_behavior) + call restFile_read(bounds_proc, fnamer, glc_behavior, & + reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns) end if else if ((nsrest == nsrContinue) .or. (nsrest == nsrBranch)) then @@ -568,8 +579,8 @@ subroutine initialize2( ) if (masterproc) then write(iulog,*)'Reading restart file ',trim(fnamer) end if - call restFile_read(bounds_proc, fnamer, glc_behavior) - + call restFile_read(bounds_proc, fnamer, glc_behavior, & + reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns) end if ! ------------------------------------------------------------------------ @@ -595,7 +606,8 @@ subroutine initialize2( ) glc_behavior=glc_behavior) ! Read new interpolated conditions file back in - call restFile_read(bounds_proc, finidat_interp_dest, glc_behavior) + call restFile_read(bounds_proc, finidat_interp_dest, glc_behavior, & + reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns) ! Reset finidat to now be finidat_interp_dest ! (to be compatible with routines still using finidat) @@ -610,30 +622,49 @@ subroutine initialize2( ) ! interpolated restart file, if applicable). ! ------------------------------------------------------------------------ - if (get_reset_dynbal_baselines()) then - if (nsrest == nsrStartup) then - if (masterproc) then - write(iulog,*) ' ' - write(iulog,*) 'Resetting dynbal baselines' - write(iulog,*) ' ' - end if - - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) - do nc = 1,nclumps - call get_clump_bounds(nc, bounds_clump) - - call dyn_hwcontent_set_baselines(bounds_clump, & - filter_inactive_and_active(nc)%num_icemecc, & - filter_inactive_and_active(nc)%icemecc, & - urbanparams_inst, soilstate_inst, water_inst, temperature_inst) - end do - else if (nsrest == nsrBranch) then + reset_dynbal_baselines_all_columns = get_reset_dynbal_baselines() + if (nsrest == nsrBranch) then + if (reset_dynbal_baselines_all_columns) then call endrun(msg='ERROR clm_initializeMod: '//& 'Cannot set reset_dynbal_baselines in a branch run') end if - ! nsrContinue not explicitly handled: it's okay for reset_dynbal_baselines to - ! remain set in a continue run, but it has no effect + else if (nsrest == nsrContinue) then + ! It's okay for the reset_dynbal_baselines flag to remain set in a continue + ! run, but we'll ignore it. (This way, the user doesn't have to change their + ! namelist file for the continue run.) + reset_dynbal_baselines_all_columns = .false. end if + ! Note that we will still honor reset_dynbal_baselines_lake_columns even in a branch + ! or continue run: even in these runs, we want to reset those baselines if they are + ! wrong on the restart file. + + if (masterproc) then + if (reset_dynbal_baselines_all_columns) then + write(iulog,*) ' ' + write(iulog,*) 'Resetting dynbal baselines for all columns' + write(iulog,*) ' ' + else if (reset_dynbal_baselines_lake_columns) then + write(iulog,*) ' ' + write(iulog,*) 'Resetting dynbal baselines for lake columns' + write(iulog,*) ' ' + end if + end if + + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + do nc = 1,nclumps + call get_clump_bounds(nc, bounds_clump) + + call dyn_hwcontent_set_baselines(bounds_clump, & + filter_inactive_and_active(nc)%num_icemecc, & + filter_inactive_and_active(nc)%icemecc, & + filter_inactive_and_active(nc)%num_lakec, & + filter_inactive_and_active(nc)%lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, & + water_inst, temperature_inst, & + reset_all_baselines = reset_dynbal_baselines_all_columns, & + reset_lake_baselines = reset_dynbal_baselines_lake_columns) + end do + !$OMP END PARALLEL DO ! ------------------------------------------------------------------------ ! Initialize nitrogen deposition diff --git a/src/main/clm_varsur.F90 b/src/main/clm_varsur.F90 index cd3bacc87a..41740f1e2b 100644 --- a/src/main/clm_varsur.F90 +++ b/src/main/clm_varsur.F90 @@ -45,6 +45,9 @@ module clm_instur ! subgrid glacier_mec sfc elevation real(r8), pointer :: topo_glc_mec(:,:) + + ! whether we have lake to initialise in each grid cell + logical , pointer :: haslake(:) !----------------------------------------------------------------------- end module clm_instur diff --git a/src/main/restFileMod.F90 b/src/main/restFileMod.F90 index e39d1e90cb..8d46e9ab83 100644 --- a/src/main/restFileMod.F90 +++ b/src/main/restFileMod.F90 @@ -25,6 +25,7 @@ module restFileMod use ncdio_pio , only : check_att, ncd_getatt use glcBehaviorMod , only : glc_behavior_type use reweightMod , only : reweight_wrapup + use IssueFixedMetadataHandler, only : write_issue_fixed_metadata, read_issue_fixed_metadata ! ! !PUBLIC TYPES: implicit none @@ -44,6 +45,8 @@ module restFileMod private :: restFile_write_pfile ! Writes restart pointer file private :: restFile_closeRestart ! Close restart file and write restart pointer file private :: restFile_dimset + private :: restFile_write_issues_fixed ! Write metadata for issues fixed + private :: restFile_read_issues_fixed ! Read and process metadata for issues fixed private :: restFile_add_flag_metadata ! Add global metadata for some logical flag private :: restFile_add_ilun_metadata ! Add global metadata defining landunit types private :: restFile_add_icol_metadata ! Add global metadata defining column types @@ -56,6 +59,9 @@ module restFileMod ! ! !PRIVATE TYPES: + ! Issue numbers for issue-fixed metadata + integer, parameter :: lake_dynbal_baseline_issue = 1140 + character(len=*), parameter, private :: sourcefile = & __FILE__ !----------------------------------------------------------------------- @@ -108,6 +114,9 @@ subroutine restFile_write( bounds, file, writing_finidat_interp_dest_file, rdate call hist_restart_ncd (bounds, ncid, flag='define', rdate=rdate ) end if + call restFile_write_issues_fixed(ncid, & + writing_finidat_interp_dest_file = writing_finidat_interp_dest_file) + call restFile_enddef( ncid ) ! Write variables @@ -142,7 +151,7 @@ subroutine restFile_write( bounds, file, writing_finidat_interp_dest_file, rdate end subroutine restFile_write !----------------------------------------------------------------------- - subroutine restFile_read( bounds_proc, file, glc_behavior ) + subroutine restFile_read( bounds_proc, file, glc_behavior, reset_dynbal_baselines_lake_columns ) ! ! !DESCRIPTION: ! Read a CLM restart file. @@ -151,6 +160,13 @@ subroutine restFile_read( bounds_proc, file, glc_behavior ) type(bounds_type) , intent(in) :: bounds_proc ! processor-level bounds character(len=*) , intent(in) :: file ! output netcdf restart file type(glc_behavior_type), intent(in) :: glc_behavior + + ! BACKWARDS_COMPATIBILITY(wjs, 2020-09-02) This is needed when reading old initial + ! conditions files created before https://github.com/ESCOMP/CTSM/issues/1140 was + ! resolved via https://github.com/ESCOMP/CTSM/pull/1109: The definition of total + ! column water content has been changed for lakes, so we need to reset baseline values + ! for lakes on older initial conditions files. + logical, intent(out) :: reset_dynbal_baselines_lake_columns ! ! !LOCAL VARIABLES: type(file_desc_t) :: ncid ! netcdf id @@ -201,6 +217,9 @@ subroutine restFile_read( bounds_proc, file, glc_behavior ) call hist_restart_ncd (bounds_proc, ncid, flag='read' ) + call restFile_read_issues_fixed(ncid, & + reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns) + ! Do error checking on file call restFile_check_consistency(bounds_proc, ncid) @@ -563,6 +582,63 @@ subroutine restFile_dimset( ncid ) end subroutine restFile_dimset + !----------------------------------------------------------------------- + subroutine restFile_write_issues_fixed(ncid, writing_finidat_interp_dest_file) + ! + ! !DESCRIPTION: + ! Write metadata for issues fixed + ! + ! !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 + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'restFile_write_issues_fixed' + !----------------------------------------------------------------------- + + ! See comment associated with reset_dynbal_baselines_lake_columns in restFile_read + call write_issue_fixed_metadata( & + ncid = ncid, & + writing_finidat_interp_dest_file = writing_finidat_interp_dest_file, & + issue_num = lake_dynbal_baseline_issue) + + end subroutine restFile_write_issues_fixed + + !----------------------------------------------------------------------- + subroutine restFile_read_issues_fixed(ncid, reset_dynbal_baselines_lake_columns) + ! + ! !DESCRIPTION: + ! Read and process metadata for issues fixed + ! + ! !ARGUMENTS: + type(file_desc_t), intent(inout) :: ncid ! local file id + logical, intent(out) :: reset_dynbal_baselines_lake_columns ! see comment in restFile_read for details + ! + ! !LOCAL VARIABLES: + integer :: attribute_value + + character(len=*), parameter :: subname = 'restFile_read_issues_fixed' + !----------------------------------------------------------------------- + + ! See comment associated with reset_dynbal_baselines_lake_columns in restFile_read + call read_issue_fixed_metadata( & + ncid = ncid, & + issue_num = lake_dynbal_baseline_issue, & + attribute_value = attribute_value) + if (attribute_value == 0) then + ! Old restart file, from before lake_dynbal_baseline_issue was resolved, so we + ! need to reset dynbal baselines for lake columns later in initialization. + reset_dynbal_baselines_lake_columns = .true. + else + ! Recent restart file, so no need to reset dynbal baselines for lake columns, + ! because they are already up to date. + reset_dynbal_baselines_lake_columns = .false. + end if + + end subroutine restFile_read_issues_fixed + + !----------------------------------------------------------------------- subroutine restFile_add_flag_metadata(ncid, flag, flag_name) ! diff --git a/src/main/subgridMod.F90 b/src/main/subgridMod.F90 index 3f34acda69..c390d9dcaa 100644 --- a/src/main/subgridMod.F90 +++ b/src/main/subgridMod.F90 @@ -38,7 +38,8 @@ module subgridMod public :: subgrid_get_info_glacier_mec public :: subgrid_get_info_crop public :: crop_patch_exists ! returns true if the given crop patch should be created in memory - + public :: lake_landunit_exists ! returns true if the lake landunit should be created in memory + ! !PRIVATE MEMBER FUNCTIONS: private :: subgrid_get_info_urban @@ -392,10 +393,10 @@ subroutine subgrid_get_info_lake(gi, npatches, ncols, nlunits) character(len=*), parameter :: subname = 'subgrid_get_info_lake' !----------------------------------------------------------------------- - ! We currently do NOT allow the lake landunit to expand via dynamic landunits, so we - ! only need to allocate space for it where its weight is currently non-zero. - - if (wt_lunit(gi, istdlak) > 0.0_r8) then + ! We do allow the lake landunit to expand via dynamic landunits, so we + ! need to allocate space for where it is known that the lake unit will grow. + + if (lake_landunit_exists(gi) ) then npatches = 1 ncols = 1 nlunits = 1 @@ -484,7 +485,6 @@ subroutine subgrid_get_info_crop(gi, npatches, ncols, nlunits) !----------------------------------------------------------------------- npatches = 0 - do cft = cft_lb, cft_ub if (crop_patch_exists(gi, cft)) then npatches = npatches + 1 @@ -558,6 +558,45 @@ function crop_patch_exists(gi, cft) result(exists) end function crop_patch_exists +!----------------------------------------------------------------------- + function lake_landunit_exists(gi) result(exists) + ! + ! !DESCRIPTION: + ! Returns true if a land unit for lakes should be created in memory + ! which is defined for gridcells which will grow lake, given by haslake + ! + ! !USES: + use dynSubgridControlMod , only : get_do_transient_lakes + use clm_instur , only : haslake + ! + ! !ARGUMENTS: + logical :: exists ! function result + integer, intent(in) :: gi ! grid cell index + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'lake_landunit_exists' + !----------------------------------------------------------------------- + + if (get_do_transient_lakes()) then + ! To support dynamic landunits, we initialise a lake land unit in each grid cell in which there are lakes. + ! This is defined by the haslake variable + + if (haslake(gi)) then + exists = .true. + else + exists = .false. + end if + + else + ! For a run without transient lakes, only allocate memory for lakes actually present in run) + if (wt_lunit(gi, istdlak) > 0.0_r8) then + exists = .true. + else + exists = .false. + end if + end if + end function lake_landunit_exists end module subgridMod diff --git a/src/main/subgridWeightsMod.F90 b/src/main/subgridWeightsMod.F90 index ddcc9585c9..3b67fd7574 100644 --- a/src/main/subgridWeightsMod.F90 +++ b/src/main/subgridWeightsMod.F90 @@ -301,7 +301,7 @@ logical function is_active_l(l, glc_behavior) ! Determine whether the given landunit is active ! ! !USES: - use landunit_varcon, only : istsoil, istice_mec, isturb_MIN, isturb_MAX + use landunit_varcon, only : istsoil, istice_mec, isturb_MIN, isturb_MAX, istdlak ! ! !ARGUMENTS: implicit none @@ -361,7 +361,15 @@ logical function is_active_l(l, glc_behavior) if (lun%itype(l) == istsoil) then is_active_l = .true. end if - + + ! Set all lake land units to active + ! By doing this, lakes are also run virtually in grid cells which will grow + ! lakes during the transient run. + + if (lun%itype(l) == istdlak) then + is_active_l = .true. + end if + end if end function is_active_l diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index 771165459d..b9d2e1bb3d 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -289,10 +289,13 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) toosmall_soil, toosmall_crop, toosmall_glacier, & toosmall_lake, toosmall_wetland, toosmall_urban, & n_dom_landunits - use fileutils , only : getfil - use domainMod , only : domain_type, domain_init, domain_clean - use clm_instur , only : wt_lunit, topo_glc_mec - use landunit_varcon, only: max_lunit, istsoil, isturb_MIN, isturb_MAX + use fileutils , only : getfil + use domainMod , only : domain_type, domain_init, domain_clean + use clm_instur , only : wt_lunit, topo_glc_mec + use landunit_varcon , only : max_lunit, istsoil, isturb_MIN, isturb_MAX + use dynSubgridControlMod, only : get_flanduse_timeseries + use dynSubgridControlMod, only : get_do_transient_lakes + ! ! !ARGUMENTS: integer, intent(in) :: begg, endg, actual_numcft @@ -312,6 +315,7 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) type(file_desc_t) :: ncid ! netcdf id logical :: istype_domain ! true => input file is of type domain logical :: isgrid2d ! true => intut grid is 2d + character(len=32) :: subname = 'surfrd_get_data' ! subroutine name !----------------------------------------------------------------------- @@ -444,7 +448,12 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) write(iulog,*) 'Successfully read surface boundary data' write(iulog,*) end if - + + ! read the lakemask (necessary for initialisation of dynamical lakes) + if (get_do_transient_lakes()) then + call surfrd_lakemask(begg, endg) + end if + end subroutine surfrd_get_data !----------------------------------------------------------------------- @@ -723,7 +732,7 @@ subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, fert_cft, cftsize, natpft wt_nat_patch(begg:,natpft_lb:natpft_size-1+natpft_lb) = array2D(begg:,:) deallocate( array2D ) - end subroutine surfrd_cftformat + end subroutine surfrd_cftformat !----------------------------------------------------------------------- subroutine surfrd_pftformat( begg, endg, ncid ) @@ -946,4 +955,72 @@ subroutine surfrd_veg_dgvm(begg, endg) end subroutine surfrd_veg_dgvm + !----------------------------------------------------------------------- + subroutine surfrd_lakemask(begg, endg) + ! + ! !DESCRIPTION: + ! Reads the lake mask, indicating where lakes are and will grow + ! of the landuse.timeseries file. + ! Necessary for the initialisation of the lake land units + ! + ! !USES: + use clm_instur , only : haslake + use dynSubgridControlMod , only : get_flanduse_timeseries + use clm_varctl , only : fname_len + use fileutils , only : getfil + ! + ! !ARGUMENTS: + integer, intent(in) :: begg, endg + ! + ! + ! !LOCAL VARIABLES: + type(file_desc_t) :: ncid_dynuse ! netcdf id for landuse timeseries file + character(len=256) :: locfn ! local file name + character(len=fname_len) :: fdynuse ! landuse.timeseries filename + logical :: readvar + real(r8),pointer :: haslake_id(:) + ! + character(len=*), parameter :: subname = 'surfrd_lakemask' + ! + !----------------------------------------------------------------------- + + ! get filename of landuse_timeseries file + fdynuse = get_flanduse_timeseries() + + if (masterproc) then + write(iulog,*) 'Attempting to read landuse.timeseries data .....' + if (fdynuse == ' ') then + write(iulog,*)'fdynuse must be specified' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + call getfil(fdynuse, locfn, 0 ) + + ! open landuse_timeseries file + call ncd_pio_openfile (ncid_dynuse, trim(locfn), 0) + + + allocate(haslake_id(begg:endg)) + + ! read the lakemask + call ncd_io(ncid=ncid_dynuse, varname='HASLAKE' , flag='read', data=haslake_id, & + dim1name=grlnd, readvar=readvar) + if (.not. readvar) call endrun( msg=' ERROR: HASLAKE is not on landuse.timeseries file'//errMsg(sourcefile, __LINE__)) + + where (haslake_id == 1._r8) + haslake = .true. + elsewhere + haslake = .false. + end where + + ! close landuse_timeseries file again + call ncd_pio_closefile(ncid_dynuse) + + deallocate(haslake_id) + + + end subroutine surfrd_lakemask + + end module surfrdMod