From b3e08240430a3e50a67e499299c3677e50adc519 Mon Sep 17 00:00:00 2001 From: Ivanderkelen Date: Wed, 20 Feb 2019 12:31:35 -0700 Subject: [PATCH 01/42] added dynlakeFileMod to account for lake area changes --- src/dyn_subgrid/dynSubgridDriverMod.F90 | 13 ++- src/dyn_subgrid/dynlakeFileMod.F90 | 142 ++++++++++++++++++++++++ 2 files changed, 154 insertions(+), 1 deletion(-) create mode 100644 src/dyn_subgrid/dynlakeFileMod.F90 diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index 73c3f9e4b0..6eb32f5b91 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -20,6 +20,8 @@ module dynSubgridDriverMod use dynpftFileMod , only : dynpft_init, dynpft_interp use dyncropFileMod , only : dyncrop_init, dyncrop_interp use dynHarvestMod , only : dynHarvest_init, dynHarvest_interp + ! use new lake file module + 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 @@ -41,7 +43,10 @@ module dynSubgridDriverMod use CropType , only : crop_type use glc2lndMod , only : glc2lnd_type use filterMod , only : filter, filter_inactive_and_active - ! + + + + ! !PUBLIC MEMBER FUNCTIONS: implicit none private @@ -120,6 +125,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 + ! error if keyword argument like for dynamical land units above: dynlake_filename=get_flanduse_timeseries() + call dynlake_init(bounds_proc, dynlake_filename=get_flanduse_timeseries()) + + ! ------------------------------------------------------------------------ ! Set initial subgrid weights for aspects that are read from file. This is relevant ! for cold start and use_init_interp-based initialization. diff --git a/src/dyn_subgrid/dynlakeFileMod.F90 b/src/dyn_subgrid/dynlakeFileMod.F90 new file mode 100644 index 0000000000..4db8720e07 --- /dev/null +++ b/src/dyn_subgrid/dynlakeFileMod.F90 @@ -0,0 +1,142 @@ +module dynlakeFileMod + + !--------------------------------------------------------------------------- + ! !DESCRIPTION: + ! Handle reading of the dataset that specifies transient areas of the lake landunit + ! + ! !USES: + + ! check this still on unneccesary parts. + +#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, namec + use abortutils , only : endrun + use spmdMod , only : masterproc, mpicom + use LandunitType , only : lun + use ColumnType , only : col + use PatchType , only : patch + + ! !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' + + ! TO DO: account for lake depth? + + 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 crop updates determined + ! prognostically, if crop 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 crops so that they are interpolated, rather + ! than jumping to each year's value on Jan 1 of that year, simply change wtcrop 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 crop landunits. + ! + ! Note that crop 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 wtcrop 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) + +! also account for landunits, patches and columns? + + + end subroutine dynlake_interp + +end module dynlakeFileMod \ No newline at end of file From eb45afd198aaa0b04b097a3e80071a2a0ba6eea0 Mon Sep 17 00:00:00 2001 From: Inne Vanderkelen Date: Fri, 22 Feb 2019 16:55:19 -0700 Subject: [PATCH 02/42] Add computation of lake heat itself in ComputeHeatMod and pass it to dynConsBiogeophysMod.f90 cv is kept constant at water cv for simplicity (altough there is a cv for each lake layer) TEMPORARY change --- src/biogeophys/LakeTemperatureMod.F90 | 4 ++ src/biogeophys/TotalWaterAndHeatMod.F90 | 70 +++++++++++++++++++++++- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 4 +- src/dyn_subgrid/dynSubgridDriverMod.F90 | 4 ++ src/dyn_subgrid/dynlakeFileMod.F90 | 4 +- 5 files changed, 82 insertions(+), 4 deletions(-) diff --git a/src/biogeophys/LakeTemperatureMod.F90 b/src/biogeophys/LakeTemperatureMod.F90 index f550f294a6..25de5f85d4 100644 --- a/src/biogeophys/LakeTemperatureMod.F90 +++ b/src/biogeophys/LakeTemperatureMod.F90 @@ -1001,6 +1001,9 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la fin(c) = fin(c) + phi(c,j) end do end do +! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts + + 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/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 9488eedc70..0fb0ec4303 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -78,6 +78,7 @@ module TotalWaterAndHeatMod ! ! !PRIVATE MEMBER FUNCTIONS: private :: AccumulateLiquidWaterHeat ! For use by ComputeHeat* routines: accumulate quantities that we need to count for liquid water, for a single column + private :: AccumulateLiquidWaterHeatLake ! same as above but able to read in 2D lake temperature private :: TempToHeat ! For use by ComputeHeat* routines: convert temperature to heat content character(len=*), parameter, private :: sourcefile = & @@ -800,7 +801,7 @@ end subroutine AccumulateSoilHeatNonLake !----------------------------------------------------------------------- subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & soilstate_inst, temperature_inst, waterstatebulk_inst, & - heat, heat_liquid, cv_liquid) + heat, heat_liquid, cv_liquid, heat_lake) ! ! !DESCRIPTION: ! Compute total heat content for all lake columns @@ -812,7 +813,13 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & ! ! Note: Changes to this routine - for anything involving liquid or ice - should ! generally be accompanied by similar changes to ComputeLiqIceMassLake - ! + + + + ! REMOVE THIS when done developing + use clm_varctl , only : iulog + + ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_lakec @@ -824,6 +831,8 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & 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)] + real(r8) , intent(inout) :: heat_lake( bounds%begc: ) ! heat content of the lake water in column [J/m^2] + ! ! !LOCAL VARIABLES: integer :: fc @@ -832,6 +841,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & real(r8) :: heat_dry_mass(bounds%begc:bounds%endc) ! sum of heat content: dry mass [J/m^2] real(r8) :: heat_ice(bounds%begc:bounds%endc) ! sum of heat content: ice [J/m^2] real(r8) :: latent_heat_liquid(bounds%begc:bounds%endc) ! sum of latent heat content of liquid water [J/m^2] + real(r8) :: latent_heat_liquid_lake(bounds%begc:bounds%endc) ! sum of latent heat content of liquid water of lake itself [J/m^2] character(len=*), parameter :: subname = 'ComputeHeatLake' !----------------------------------------------------------------------- @@ -843,8 +853,10 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl, & ! number of snow layers dz => col%dz, & ! layer depth (m) + lakedepth => col%lakedepth, & ! lake depth (m) watsat => soilstate_inst%watsat_col, & ! volumetric soil water at saturation (porosity) csol => soilstate_inst%csol_col, & ! heat capacity, soil solids (J/m**3/Kelvin) + t_lake => temperature_inst%t_lake_col, & ! lake temperature (K) 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) @@ -902,6 +914,21 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & ! 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) + ! calculate heat content of lake itself + !set condition on heat_lake when is present (optional argument)? + do fc = 1, num_lakec + c = filter_lakec(fc) + call AccumulateLiquidWaterHeatLake( & + temp = t_lake(c,:), & + h2o = lakedepth(c)*1000, & + cv_liquid = cv_liquid(c), & + heat_liquid = heat_lake(c), & + latent_heat_liquid = latent_heat_liquid_lake(c)) + end do + + write(iulog,*) 'lake heat (J/m^2)', heat_lake(c) + +! Add lake heat here if wanted to incorporate 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) @@ -1087,6 +1114,45 @@ subroutine AccumulateLiquidWaterHeat(temp, h2o, & latent_heat_liquid = latent_heat_liquid + h2o*hfus end subroutine AccumulateLiquidWaterHeat +! Repeating same module, quick fix to account for 2D temperature field of lake temp. (to be removed) + !----------------------------------------------------------------------- + subroutine AccumulateLiquidWaterHeatLake(temp, h2o, & + heat_liquid, latent_heat_liquid, cv_liquid) + ! + ! !DESCRIPTION: + ! In the course of accumulating heat contents: Accumulate quantities that we need to + ! count for liquid water, for a single column + + use clm_varpar ,only : nlevlak +! + ! !ARGUMENTS: + real(r8), intent(in) :: temp(:) ! temperature [K] + real(r8), intent(in) :: h2o ! water mass [kg/m^2] + + real(r8), intent(inout) :: heat_liquid ! accumulated total heat content of liquid water for this column, excluding latent heat [J/m^2] + real(r8), intent(inout) :: latent_heat_liquid ! accumulated total latent heat content of liquid water for this column [J/m^2] + real(r8), intent(inout), optional :: cv_liquid ! accumulated total liquid heat capacity for this column [J/(m^2 K)] + ! + ! !LOCAL VARIABLES: + real(r8) :: cv ! heat capacity [J/(m^2 K)] + integer :: j ! do loop index + character(len=*), parameter :: subname = 'AccumulateLiquidWaterHeatLake' + !----------------------------------------------------------------------- + + cv = h2o*cpliq + if (present(cv_liquid)) then + cv_liquid = cv_liquid + cv + end if + + ! loop over lake levels + do j = 1,nlevlak + heat_liquid = heat_liquid + TempToHeat(temp = temp(j), cv = cv) + end do + latent_heat_liquid = latent_heat_liquid + h2o*hfus + + end subroutine AccumulateLiquidWaterHeatLake + + !----------------------------------------------------------------------- pure function TempToHeat(temp, cv) result(heat) ! diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 018cd12176..9c2e6bdf34 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -642,6 +642,7 @@ subroutine dyn_heat_content(bounds, & real(r8) :: heat_col(bounds%begc:bounds%endc) ! sum of heat content for all columns [J/m^2] real(r8) :: heat_liquid_col(bounds%begc:bounds%endc) ! sum of heat content for all columns: liquid water, excluding latent heat [J/m^2] real(r8) :: cv_liquid_col(bounds%begc:bounds%endc) ! sum of liquid heat capacity for all columns [J/(m^2 K)] + real(r8) :: heat_lake_col(bounds%begc:bounds%endc) ! heat content of lake column [J/m^2] real(r8) :: heat_liquid_grc(bounds%begg:bounds%endg) ! heat_liquid_col averaged to grid cell [J/m^2] real(r8) :: cv_liquid_grc(bounds%begg:bounds%endg) ! cv_liquid_col averaged to grid cell [J/(m^2 K)] @@ -662,7 +663,8 @@ subroutine dyn_heat_content(bounds, & soilstate_inst, temperature_inst, waterstatebulk_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)) + cv_liquid = cv_liquid_col(bounds%begc:bounds%endc), & + heat_lake = heat_lake_col(bounds%begc:bounds%endc)) call c2g(bounds, & carr = heat_col(bounds%begc:bounds%endc), & diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index 6eb32f5b91..b7528c1b24 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -255,6 +255,10 @@ subroutine dynSubgrid_driver(bounds_proc, if (get_do_harvest()) then call dynHarvest_interp(bounds_proc) end if + + ! add lake interp (condition to be added later) + + call dynlake_interp(bounds_proc) ! ========================================================================== ! Do land cover change that does not require I/O diff --git a/src/dyn_subgrid/dynlakeFileMod.F90 b/src/dyn_subgrid/dynlakeFileMod.F90 index 4db8720e07..f3179b1f55 100644 --- a/src/dyn_subgrid/dynlakeFileMod.F90 +++ b/src/dyn_subgrid/dynlakeFileMod.F90 @@ -54,6 +54,7 @@ subroutine dynlake_init(bounds, dynlake_filename) ! !USES: use ncdio_pio , only : check_dim use dynTimeInfoMod , only : YEAR_POSITION_START_OF_TIMESTEP + ! ! !ARGUMENTS: @@ -62,7 +63,7 @@ subroutine dynlake_init(bounds, dynlake_filename) ! ! !LOCAL VARIABLES: integer :: num_points ! number of spatial points - + character(len=*), parameter :: subname = 'dynlake_init' !----------------------------------------------------------------------- @@ -91,6 +92,7 @@ subroutine dynlake_init(bounds, dynlake_filename) 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 From 91ac079a897c59e4ff622e977b7bb6e7cf933ca6 Mon Sep 17 00:00:00 2001 From: Inne Vanderkelen Date: Fri, 1 Mar 2019 05:28:50 -0700 Subject: [PATCH 03/42] Make lake heat output variable of LakeTemperatureMod.F90 and print is as an history field. The history field only contains the heat content of the lake water itself, and is a column output. --- src/biogeophys/LakeTemperatureMod.F90 | 9 ++++++--- src/biogeophys/TemperatureType.F90 | 9 +++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/biogeophys/LakeTemperatureMod.F90 b/src/biogeophys/LakeTemperatureMod.F90 index 25de5f85d4..26f8ee9134 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,8 @@ 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) - + lake_heat => temperature_inst%lake_heat , & ! Output: [real(r8) (:) ] col lake heat (J/m²) + 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,11 +1001,12 @@ 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 -! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts - + write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts + lake_heat(c) = ncvts(c) call waterstatebulk_inst%CalculateTotalH2osno(bounds, num_lakec, filter_lakec, & caller = 'LakeTemperature-2', & diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 5c69733f8a..5dacf3fa55 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -117,6 +117,9 @@ module TemperatureType real(r8), pointer :: fact_col (:,:) ! used in computing tridiagonal matrix real(r8), pointer :: c_h2osfc_col (:) ! heat capacity of surface water + ! lake heat + real(r8), pointer :: lake_heat (:) ! total heat of lake water (J/m²) + contains procedure, public :: Init @@ -279,6 +282,7 @@ subroutine InitAllocate(this, bounds) allocate(this%fact_col (begc:endc, -nlevsno+1:nlevgrnd)) ; this%fact_col (:,:) = nan allocate(this%c_h2osfc_col (begc:endc)) ; this%c_h2osfc_col (:) = nan + allocate(this%lake_heat (begc:endc)) ; this%lake_heat (:) = nan end subroutine InitAllocate !------------------------------------------------------------------------ @@ -620,6 +624,11 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) ptr_patch=this%t_veg10_night_patch, default='inactive') endif + ! add lake heat history field here + this%lake_heat(begc:endc) = spval + call hist_addfld1d (fname='LAKE_HEAT', units='J/m^2', & + avgflag='A', long_name='Heat content of gridcell lake water', & + ptr_col=this%lake_heat, default='active') end subroutine InitHistory From 36b3f1522c077156044f9f6f945ef51c8ebed035 Mon Sep 17 00:00:00 2001 From: Inne Vanderkelen Date: Fri, 8 Mar 2019 03:48:53 -0700 Subject: [PATCH 04/42] Made handle do_transient_lakes to use in namelist for enabling the dynamic lake land unit. Use similar to do_transient_crops --- .../namelist_definition_ctsm.xml | 6 ++++ src/dyn_subgrid/dynSubgridControlMod.F90 | 32 +++++++++++++++++++ src/dyn_subgrid/dynSubgridDriverMod.F90 | 9 +++--- 3 files changed, 42 insertions(+), 5 deletions(-) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 5d142e4626..45c39abeaa 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2426,6 +2426,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/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index 3820d6392a..01ff35ab33 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)' @@ -277,6 +289,14 @@ subroutine check_namelist_consistency end if end if + if (dyn_subgrid_control_inst%do_transient_lakes) then + if (use_fates) then + write(iulog,*) 'ERROR: do_transient_lakes currently does not work with use_fates' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + end subroutine check_namelist_consistency !----------------------------------------------------------------------- @@ -316,6 +336,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 b7528c1b24..c162272e33 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 @@ -256,10 +256,9 @@ subroutine dynSubgrid_driver(bounds_proc, call dynHarvest_interp(bounds_proc) end if - ! add lake interp (condition to be added later) - - call dynlake_interp(bounds_proc) - + if (get_do_transient_lakes()) then + call dynlake_interp(bounds_proc) + end if ! ========================================================================== ! Do land cover change that does not require I/O ! ========================================================================== From 35c16afa909c66874476a4cd4bceffc64cab8a03 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 20 Jun 2019 08:29:38 -0600 Subject: [PATCH 05/42] Minor changes to source code (This commit was split off from 60c0aafc7 using interactive rebase, discarding the tools changes and just keeping the relatively small source code changes in that commit.) --- src/biogeophys/TotalWaterAndHeatMod.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 0fb0ec4303..292bbcb700 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -859,6 +859,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & t_lake => temperature_inst%t_lake_col, & ! lake temperature (K) 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) + lake_heat => temperature_inst%lake_heat, & ! total heat of lake water (J/m²) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col, & ! liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! frozen water (kg/m2) ) @@ -926,7 +927,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & latent_heat_liquid = latent_heat_liquid_lake(c)) end do - write(iulog,*) 'lake heat (J/m^2)', heat_lake(c) + write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c) ! Add lake heat here if wanted to incorporate do fc = 1, num_lakec @@ -1148,7 +1149,10 @@ subroutine AccumulateLiquidWaterHeatLake(temp, h2o, & do j = 1,nlevlak heat_liquid = heat_liquid + TempToHeat(temp = temp(j), cv = cv) end do - latent_heat_liquid = latent_heat_liquid + h2o*hfus + + ! this would assume the whole lake unfrozen? + latent_heat_liquid = latent_heat_liquid + h2o*hfus + end subroutine AccumulateLiquidWaterHeatLake From 587d8f99300f883b7a992907ebac4cef07b3f079 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 1 Aug 2019 01:44:24 -0600 Subject: [PATCH 06/42] Added lake landunit initialisation for gridcells which will grow lake --- src/biogeophys/LakeTemperatureMod.F90 | 5 +- src/biogeophys/TotalWaterAndHeatMod.F90 | 2 +- src/main/clm_initializeMod.F90 | 7 ++- src/main/clm_varsur.F90 | 3 + src/main/subgridMod.F90 | 51 +++++++++++++-- src/main/subgridWeightsMod.F90 | 12 +++- src/main/surfrdMod.F90 | 84 +++++++++++++++++++++++-- 7 files changed, 146 insertions(+), 18 deletions(-) diff --git a/src/biogeophys/LakeTemperatureMod.F90 b/src/biogeophys/LakeTemperatureMod.F90 index 26f8ee9134..6a20856f67 100644 --- a/src/biogeophys/LakeTemperatureMod.F90 +++ b/src/biogeophys/LakeTemperatureMod.F90 @@ -1004,9 +1004,10 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la end do end do - write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts + ! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts - lake_heat(c) = ncvts(c) + ! IV: currently commented out: caused crash. To do: look at this part of the code!!! + ! lake_heat(c) = ncvts(c) call waterstatebulk_inst%CalculateTotalH2osno(bounds, num_lakec, filter_lakec, & caller = 'LakeTemperature-2', & diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 292bbcb700..8f436687cf 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -927,7 +927,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & latent_heat_liquid = latent_heat_liquid_lake(c)) end do - write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c) + ! write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c) ! Add lake heat here if wanted to incorporate do fc = 1, num_lakec diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 442a77c36b..bdb89df7c2 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 @@ -181,7 +181,7 @@ subroutine initialize1(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 @@ -263,7 +263,7 @@ subroutine initialize1(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') @@ -276,6 +276,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 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/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..f73265cd92 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -289,10 +289,12 @@ 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 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 clm_varctl , only : fname_len ! ! !ARGUMENTS: integer, intent(in) :: begg, endg, actual_numcft @@ -310,8 +312,11 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) logical :: readvar ! true => variable is on dataset real(r8) :: rmaxlon,rmaxlat ! local min/max vars type(file_desc_t) :: ncid ! netcdf id + type(file_desc_t) :: ncid_dynuse ! netcdf id for landuse timeseries file logical :: istype_domain ! true => input file is of type domain logical :: isgrid2d ! true => intut grid is 2d + character(len=fname_len) :: fdynuse ! landuse.timeseries filename + character(len=32) :: subname = 'surfrd_get_data' ! subroutine name !----------------------------------------------------------------------- @@ -444,6 +449,33 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) write(iulog,*) 'Successfully read surface boundary data' write(iulog,*) end if + + + ! IV: add here call to subroutine to read in lake mask (necessary for initialisation of dynamical lakes) + ! First open landuse.timeseries file for this. + + if (masterproc) then + write(iulog,*) 'Attempting to read landuse.timeseries data .....' + if (lfsurdat == ' ') then + write(iulog,*)'fdynuse must be specified' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + endif + + + ! open landuse_timeseries file + fdynuse = get_flanduse_timeseries() + + call getfil(fdynuse, locfn, 0 ) + + call ncd_pio_openfile (ncid_dynuse, trim(locfn), 0) + + ! read the lakemask + call surfrd_lakemask(begg, endg, ncid_dynuse) + + ! close landuse_timeseries file again + call ncd_pio_closefile(ncid_dynuse) + end subroutine surfrd_get_data @@ -723,7 +755,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 +978,48 @@ subroutine surfrd_veg_dgvm(begg, endg) end subroutine surfrd_veg_dgvm + !----------------------------------------------------------------------- + subroutine surfrd_lakemask(begg, endg, ncid) + ! + ! !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 + ! !ARGUMENTS: + integer, intent(in) :: begg, endg + type(file_desc_t), intent(inout) :: ncid + ! + ! + ! !LOCAL VARIABLES: + logical :: readvar + integer ,pointer :: haslake_id(:) + ! + ! + character(len=*), parameter :: subname = 'surfrd_lakemask' + ! + !----------------------------------------------------------------------- + + + allocate(haslake_id(begg:endg)) + + call ncd_io(ncid=ncid, 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.) + haslake = .true. + elsewhere + haslake = .false. + end where + + + deallocate(haslake_id) + + + end subroutine surfrd_lakemask + + end module surfrdMod From 63b2452563cff0227482fb69eed6bd826f73d7e2 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Mon, 16 Sep 2019 05:11:52 -0600 Subject: [PATCH 07/42] Add lake water and heat content to gridcell total water and heat computation --- src/biogeophys/BalanceCheckMod.F90 | 6 +- src/biogeophys/LakeHydrologyMod.F90 | 4 +- src/biogeophys/TotalWaterAndHeatMod.F90 | 130 ++++++++++------------- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 24 +++-- src/dyn_subgrid/dynSubgridDriverMod.F90 | 8 +- src/main/clm_driver.F90 | 4 +- 6 files changed, 87 insertions(+), 89 deletions(-) diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index c2d58f711e..65844b7299 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 ! @@ -210,7 +212,7 @@ subroutine BeginWaterBalanceSingle(bounds, & water_mass = begwb(bounds%begc:bounds%endc)) call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, & + waterstate_inst, lakestate_inst, & subtract_dynbal_baselines = .false., & water_mass = begwb(bounds%begc:bounds%endc)) diff --git a/src/biogeophys/LakeHydrologyMod.F90 b/src/biogeophys/LakeHydrologyMod.F90 index 44588fafff..c078e958e1 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,7 +649,7 @@ subroutine LakeHydrology(bounds, & ! Determine ending water balance and volumetric soil water call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - b_waterstate_inst, & + b_waterstate_inst, lakestate_inst, & subtract_dynbal_baselines = .false., & water_mass = endwb(bounds%begc:bounds%endc)) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 8f436687cf..e176555a14 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 @@ -78,7 +79,6 @@ module TotalWaterAndHeatMod ! ! !PRIVATE MEMBER FUNCTIONS: private :: AccumulateLiquidWaterHeat ! For use by ComputeHeat* routines: accumulate quantities that we need to count for liquid water, for a single column - private :: AccumulateLiquidWaterHeatLake ! same as above but able to read in 2D lake temperature private :: TempToHeat ! For use by ComputeHeat* routines: convert temperature to heat content character(len=*), parameter, private :: sourcefile = & @@ -140,7 +140,7 @@ end subroutine ComputeWaterMassNonLake !----------------------------------------------------------------------- subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, & + waterstate_inst, lakestate_inst, & subtract_dynbal_baselines, & water_mass) ! @@ -155,6 +155,7 @@ 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. @@ -177,6 +178,7 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & num_lakec = num_lakec, & filter_lakec = filter_lakec, & waterstate_inst = waterstate_inst, & + lakestate_inst = lakestate_inst, & subtract_dynbal_baselines = subtract_dynbal_baselines, & liquid_mass = liquid_mass(bounds%begc:bounds%endc), & ice_mass = ice_mass(bounds%begc:bounds%endc)) @@ -383,7 +385,7 @@ end subroutine AccumulateSoilLiqIceMassNonLake !----------------------------------------------------------------------- subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, & + waterstate_inst, lakestate_inst,& subtract_dynbal_baselines, & liquid_mass, ice_mass) ! @@ -401,6 +403,8 @@ 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. @@ -410,8 +414,10 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & 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 + real(r8) :: h2olak_liq(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] + real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] + character(len=*), parameter :: subname = 'ComputeLiqIceMassLake' !----------------------------------------------------------------------- @@ -420,10 +426,11 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl , & ! Input: [integer (:) ] negative number of snow layers - + dz_lake => col%dz_lake , & ! Input: [real(r8) (:,:) ] lake depth (m) 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) + lake_icefrac => lakestate_inst%lake_icefrac_col, & ! Input: [real(r8) (:,:) ] lake ice fraction 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) ) @@ -445,6 +452,19 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & end do end do + ! Lake water content + do j = 1, nlevlak + do fc = 1, num_lakec + c = filter_lakec(fc) + ! calculate lake liq and ice content per lake layer first + h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) + + liquid_mass(c) = liquid_mass(c) + h2olak_liq(c,j) + ice_mass(c) = ice_mass(c) + h2olak_ice(c,j) + end do + end do + ! Soil water content of the soil under the lake do j = 1, nlevgrnd do fc = 1, num_lakec @@ -800,8 +820,8 @@ end subroutine AccumulateSoilHeatNonLake !----------------------------------------------------------------------- subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & - soilstate_inst, temperature_inst, waterstatebulk_inst, & - heat, heat_liquid, cv_liquid, heat_lake) + soilstate_inst, temperature_inst, waterstatebulk_inst, lakestate_inst, & + heat, heat_liquid, cv_liquid) ! ! !DESCRIPTION: ! Compute total heat content for all lake columns @@ -816,10 +836,6 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & - ! REMOVE THIS when done developing - use clm_varctl , only : iulog - - ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_lakec @@ -827,11 +843,12 @@ 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)] - real(r8) , intent(inout) :: heat_lake( bounds%begc: ) ! heat content of the lake water in column [J/m^2] ! ! !LOCAL VARIABLES: @@ -841,8 +858,9 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & real(r8) :: heat_dry_mass(bounds%begc:bounds%endc) ! sum of heat content: dry mass [J/m^2] real(r8) :: heat_ice(bounds%begc:bounds%endc) ! sum of heat content: ice [J/m^2] real(r8) :: latent_heat_liquid(bounds%begc:bounds%endc) ! sum of latent heat content of liquid water [J/m^2] - real(r8) :: latent_heat_liquid_lake(bounds%begc:bounds%endc) ! sum of latent heat content of liquid water of lake itself [J/m^2] - + real(r8) :: h2olak_liq(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] + real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] + character(len=*), parameter :: subname = 'ComputeHeatLake' !----------------------------------------------------------------------- @@ -853,7 +871,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl, & ! number of snow layers dz => col%dz, & ! layer depth (m) - lakedepth => col%lakedepth, & ! lake depth (m) + dz_lake => col%dz_lake, & ! lake layer depth (m) watsat => soilstate_inst%watsat_col, & ! volumetric soil water at saturation (porosity) csol => soilstate_inst%csol_col, & ! heat capacity, soil solids (J/m**3/Kelvin) t_lake => temperature_inst%t_lake_col, & ! lake temperature (K) @@ -861,7 +879,8 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & 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) lake_heat => temperature_inst%lake_heat, & ! total heat of lake water (J/m²) 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) + lake_icefrac => lakestate_inst%lake_icefrac_col & ! Input: [real(r8) (:,:) ] lake ice fraction ) do fc = 1, num_lakec @@ -915,21 +934,27 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & ! 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) - ! calculate heat content of lake itself - !set condition on heat_lake when is present (optional argument)? - do fc = 1, num_lakec - c = filter_lakec(fc) - call AccumulateLiquidWaterHeatLake( & - temp = t_lake(c,:), & - h2o = lakedepth(c)*1000, & - cv_liquid = cv_liquid(c), & - heat_liquid = heat_lake(c), & - latent_heat_liquid = latent_heat_liquid_lake(c)) - end do + ! calculate heat content of lake itself + do j = 1, nlevlak + do fc = 1, num_lakec + c = filter_lakec(fc) + ! liquid heat + h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + call AccumulateLiquidWaterHeat( & + temp = t_lake(c,j), & + h2o = h2olak_liq(c,j), & + cv_liquid = cv_liquid(c), & + heat_liquid = heat_liquid(c), & + latent_heat_liquid = latent_heat_liquid(c)) + ! ice heat + h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) + heat_ice(c) = heat_ice(c) + & + TempToHeat(temp=t_lake(c,j), cv = (h2olak_ice(c,j) * cpice)) + end do + end do ! write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c) -! Add lake heat here if wanted to incorporate 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) @@ -1115,47 +1140,6 @@ subroutine AccumulateLiquidWaterHeat(temp, h2o, & latent_heat_liquid = latent_heat_liquid + h2o*hfus end subroutine AccumulateLiquidWaterHeat -! Repeating same module, quick fix to account for 2D temperature field of lake temp. (to be removed) - !----------------------------------------------------------------------- - subroutine AccumulateLiquidWaterHeatLake(temp, h2o, & - heat_liquid, latent_heat_liquid, cv_liquid) - ! - ! !DESCRIPTION: - ! In the course of accumulating heat contents: Accumulate quantities that we need to - ! count for liquid water, for a single column - - use clm_varpar ,only : nlevlak -! - ! !ARGUMENTS: - real(r8), intent(in) :: temp(:) ! temperature [K] - real(r8), intent(in) :: h2o ! water mass [kg/m^2] - - real(r8), intent(inout) :: heat_liquid ! accumulated total heat content of liquid water for this column, excluding latent heat [J/m^2] - real(r8), intent(inout) :: latent_heat_liquid ! accumulated total latent heat content of liquid water for this column [J/m^2] - real(r8), intent(inout), optional :: cv_liquid ! accumulated total liquid heat capacity for this column [J/(m^2 K)] - ! - ! !LOCAL VARIABLES: - real(r8) :: cv ! heat capacity [J/(m^2 K)] - integer :: j ! do loop index - character(len=*), parameter :: subname = 'AccumulateLiquidWaterHeatLake' - !----------------------------------------------------------------------- - - cv = h2o*cpliq - if (present(cv_liquid)) then - cv_liquid = cv_liquid + cv - end if - - ! loop over lake levels - do j = 1,nlevlak - heat_liquid = heat_liquid + TempToHeat(temp = temp(j), cv = cv) - end do - - ! this would assume the whole lake unfrozen? - latent_heat_liquid = latent_heat_liquid + h2o*hfus - - - end subroutine AccumulateLiquidWaterHeatLake - !----------------------------------------------------------------------- pure function TempToHeat(temp, cv) result(heat) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 9c2e6bdf34..95deb10dd8 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -20,6 +20,7 @@ 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 @@ -334,7 +335,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 +351,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 +372,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 +384,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 +398,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 +414,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 +439,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 +453,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)) @@ -549,7 +557,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 +571,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,6 +593,7 @@ subroutine dyn_water_content(bounds, & call ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & waterstate_inst, & + lakestate_inst, & subtract_dynbal_baselines = .true., & liquid_mass = liquid_mass_col(bounds%begc:bounds%endc), & ice_mass = ice_mass_col(bounds%begc:bounds%endc)) @@ -608,7 +618,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 +641,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) @@ -660,11 +671,10 @@ 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), & - heat_lake = heat_lake_col(bounds%begc:bounds%endc)) + cv_liquid = cv_liquid_col(bounds%begc:bounds%endc)) call c2g(bounds, & carr = heat_col(bounds%begc:bounds%endc), & diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index c162272e33..f5ace35b48 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -39,6 +39,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 @@ -164,7 +165,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, & @@ -192,6 +193,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 @@ -232,7 +234,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) @@ -306,7 +308,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/main/clm_driver.F90 b/src/main/clm_driver.F90 index 92c2372a52..60f64ab03d 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -332,7 +332,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, & @@ -378,7 +378,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') From 48cddcc7830f51357ec164349f183b59fc5af457 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 19 Sep 2019 03:29:21 -0600 Subject: [PATCH 08/42] Set initial subgrid weights for aspects that are read from file: add lakes --- src/dyn_subgrid/dynSubgridDriverMod.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index f5ace35b48..c82598d444 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -145,6 +145,9 @@ 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.) From e9a2709fee992647c5da3197b5acd306e842d833 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 19 Sep 2019 03:31:18 -0600 Subject: [PATCH 09/42] Fix arguments to account for total lake heat and water content --- src/main/clm_initializeMod.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index bdb89df7c2..7f596c14f0 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -486,7 +486,9 @@ 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) end do ! ------------------------------------------------------------------------ @@ -624,7 +626,10 @@ 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) end do else if (nsrest == nsrBranch) then call endrun(msg='ERROR clm_initializeMod: '//& From 50bf1975c14e6f599075f45446feed818158c749 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 19 Sep 2019 09:09:23 -0600 Subject: [PATCH 10/42] Add baseline approach for lake dynamical land units --- src/biogeophys/TotalWaterAndHeatMod.F90 | 208 +++++++++++++++++------ src/dyn_subgrid/dynConsBiogeophysMod.F90 | 102 ++++++++--- src/dyn_subgrid/dynSubgridDriverMod.F90 | 2 + 3 files changed, 237 insertions(+), 75 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index e176555a14..60af2a901e 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -39,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 @@ -155,7 +157,7 @@ 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 + 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. @@ -415,9 +417,7 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & ! ! !LOCAL VARIABLES: integer :: c, j, fc ! indices - real(r8) :: h2olak_liq(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] - real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] - + character(len=*), parameter :: subname = 'ComputeLiqIceMassLake' !----------------------------------------------------------------------- @@ -426,11 +426,9 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl , & ! Input: [integer (:) ] negative number of snow layers - dz_lake => col%dz_lake , & ! Input: [real(r8) (:,:) ] lake depth (m) 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) - lake_icefrac => lakestate_inst%lake_icefrac_col, & ! Input: [real(r8) (:,:) ] lake ice fraction 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) ) @@ -451,19 +449,10 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) end do end do - - ! Lake water content - do j = 1, nlevlak - do fc = 1, num_lakec - c = filter_lakec(fc) - ! calculate lake liq and ice content per lake layer first - h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) - h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) - - liquid_mass(c) = liquid_mass(c) + h2olak_liq(c,j) - ice_mass(c) = ice_mass(c) + h2olak_ice(c,j) - end do - end do + + ! Lake water content + call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, & + lakestate_inst, liquid_mass, ice_mass) ! Soil water content of the soil under the lake do j = 1, nlevgrnd @@ -487,6 +476,60 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & end subroutine ComputeLiqIceMassLake + !----------------------------------------------------------------------- + subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & + lakestate_inst, 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 + ! AccumulateSoilHeatLake + ! + ! !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(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(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] + real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] + + character(len=*), parameter :: subname = 'AccumulateLiqIceMassLake' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL((ubound(liquid_mass) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(ice_mass) == [bounds%endc]), errMsg(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) + ! calculate lake liq and ice content per lake layer first + h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) + + liquid_mass(c) = liquid_mass(c) + h2olak_liq(c,j) + ice_mass(c) = ice_mass(c) + h2olak_ice(c,j) + end do + end do + + end associate + + end subroutine AccumulateLiqIceMassLake + + !----------------------------------------------------------------------- subroutine ComputeHeatNonLake(bounds, num_nolakec, filter_nolakec, & urbanparams_inst, soilstate_inst, & @@ -858,9 +901,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & real(r8) :: heat_dry_mass(bounds%begc:bounds%endc) ! sum of heat content: dry mass [J/m^2] real(r8) :: heat_ice(bounds%begc:bounds%endc) ! sum of heat content: ice [J/m^2] real(r8) :: latent_heat_liquid(bounds%begc:bounds%endc) ! sum of latent heat content of liquid water [J/m^2] - real(r8) :: h2olak_liq(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] - real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] - + character(len=*), parameter :: subname = 'ComputeHeatLake' !----------------------------------------------------------------------- @@ -871,16 +912,13 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl, & ! number of snow layers dz => col%dz, & ! layer depth (m) - dz_lake => col%dz_lake, & ! lake layer depth (m) watsat => soilstate_inst%watsat_col, & ! volumetric soil water at saturation (porosity) csol => soilstate_inst%csol_col, & ! heat capacity, soil solids (J/m**3/Kelvin) - t_lake => temperature_inst%t_lake_col, & ! lake temperature (K) 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) lake_heat => temperature_inst%lake_heat, & ! total heat of lake water (J/m²) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col, & ! liquid water (kg/m2) - h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col, & ! frozen water (kg/m2) - lake_icefrac => lakestate_inst%lake_icefrac_col & ! Input: [real(r8) (:,:) ] lake ice fraction + h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! frozen water (kg/m2) ) do fc = 1, num_lakec @@ -913,7 +951,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) @@ -930,35 +968,15 @@ 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) - - ! calculate heat content of lake itself - do j = 1, nlevlak + do fc = 1, num_lakec - c = filter_lakec(fc) - ! liquid heat - h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) - call AccumulateLiquidWaterHeat( & - temp = t_lake(c,j), & - h2o = h2olak_liq(c,j), & - cv_liquid = cv_liquid(c), & - heat_liquid = heat_liquid(c), & - latent_heat_liquid = latent_heat_liquid(c)) - ! ice heat - h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) - heat_ice(c) = heat_ice(c) + & - TempToHeat(temp=t_lake(c,j), cv = (h2olak_ice(c,j) * cpice)) - end do - end do - - ! write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c) - - 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) end do + + ! Lake water heat content + call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, & + heat, heat_liquid, cv_liquid) ! Subtract baselines set in initialization ! @@ -982,6 +1000,92 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & end subroutine ComputeHeatLake + !----------------------------------------------------------------------- + subroutine AccumulateHeatLake(bounds, num_c, filter_c, & + temperature_inst, lakestate_inst, & + heat, heat_liquid, cv_liquid) + ! + ! !DESCRIPTION: + ! Accumulate heat of lake water in lake columns. + ! + ! Adds to any existing values in heat, heat_liquid and cv_liquid. + ! + ! Note: Changes to this routine should generally be accompanied by similar changes to + ! AccumulateLiqIceMassLake + ! + ! !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] + real(r8) , intent(inout) :: heat_liquid( bounds%begc: ) ! accumulated heat content: liquid water, excluding latent heat [J/m^2] + real(r8) , intent(inout) :: cv_liquid( bounds%begc: ) ! accumulated liquid water heat capacity [J/(m^2 K)] + ! + ! !LOCAL VARIABLES: + integer :: fc + integer :: l, c, j + real(r8) :: h2olak_liq(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] + real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! 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((ubound(heat) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(heat_liquid) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(cv_liquid) == [bounds%endc]), errMsg(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(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + call AccumulateLiquidWaterHeat( & + temp = t_lake(c,j), & + h2o = h2olak_liq(c,j), & + cv_liquid = cv_liquid(c), & + heat_liquid = lake_heat_liquid(c), & + latent_heat_liquid = lake_latent_heat_liquid(c)) + ! ice heat + h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) + lake_heat_ice(c) = lake_heat_ice(c) + & + TempToHeat(temp=t_lake(c,j), cv = (h2olak_ice(c,j) * cpice)) + end do + end do + + ! add ice heat and liquid heat together (look at this part) + do fc = 1, num_c + c = filter_c(fc) + heat_liquid(c) = heat_liquid(c) + lake_heat_liquid(c) + heat(c) = heat(c) + lake_heat_ice(c) + & + lake_heat_liquid(c) + lake_latent_heat_liquid(c) + end do + + end associate + + end subroutine AccumulateHeatLake + !----------------------------------------------------------------------- subroutine AdjustDeltaHeatForDeltaLiq(bounds, delta_liq, & liquid_water_temp1, liquid_water_temp2, & diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 95deb10dd8..b30969689f 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -20,13 +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 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 @@ -65,7 +65,8 @@ 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) ! ! !DESCRIPTION: ! Set start-of-run baseline values for heat and water content in some columns. @@ -94,9 +95,13 @@ 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 ! @@ -119,22 +124,23 @@ 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) 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, & + num_icemecc, filter_icemecc, num_lakec, filter_lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, water_inst%waterstatebulk_inst, & temperature_inst) + 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) ! ! !DESCRIPTION: ! Set start-of-run baseline values for water content, for a single water tracer or @@ -147,13 +153,19 @@ 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 ! ! !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' !----------------------------------------------------------------------- @@ -164,7 +176,10 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8 soil_ice_mass_col (bounds%begc:bounds%endc) = 0._r8 - + lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8 + lake_ice_mass_col (bounds%begc:bounds%endc) = 0._r8 + + call AccumulateSoilLiqIceMassNonLake(bounds, & natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, & waterstate_inst, & @@ -184,15 +199,31 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & 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)) - + + + ! set baselines for lake columns + + ! Calculate the total water volume of the lake column + call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, & + 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 associate 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, & + num_icemecc, filter_icemecc, num_lakec, filter_lakec, & + urbanparams_inst, soilstate_inst, lakestate_inst, waterstatebulk_inst, & temperature_inst) ! ! !DESCRIPTION: @@ -203,18 +234,25 @@ 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 ! ! !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] + real(r8) :: lake_heat_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface + real(r8) :: lake_cv_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface character(len=*), parameter :: subname = 'dyn_heat_content_set_baselines' !----------------------------------------------------------------------- @@ -227,12 +265,17 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8 soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8 + lake_heat_col(bounds%begc:bounds%endc) = 0._r8 + lake_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8 + lake_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 @@ -251,7 +294,20 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & 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)) - + + + ! Set baselines for lake columns + call AccumulateHeatLake(bounds, num_lakec, filter_lakec, & + temperature_inst, lakestate_inst, & + heat = lake_heat_col, heat_liquid = lake_heat_liquid_col, cv_liquid = lake_cv_liquid_col) + + do fc = 1, num_lakec + c = filter_lakec(fc) + + dynbal_baseline_heat(c) = lake_heat_col(c) + + end do + end associate end subroutine dyn_heat_content_set_baselines diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index c82598d444..69bbb388a8 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -148,6 +148,8 @@ subroutine dynSubgrid_init(bounds_proc, glc_behavior, crop_inst) 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.) From 079e5c47f1d68ba97fca8298ba496cc57b9750cf Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Fri, 20 Sep 2019 08:01:45 -0600 Subject: [PATCH 11/42] Removed heat_liquid and cv_liquid calculation for lakes --still to add comments to code! --- src/biogeophys/TotalWaterAndHeatMod.F90 | 17 ++++++++--------- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 10 +++------- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 60af2a901e..5b5a7318ee 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -976,7 +976,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & ! Lake water heat content call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, & - heat, heat_liquid, cv_liquid) + heat) ! Subtract baselines set in initialization ! @@ -1003,12 +1003,14 @@ end subroutine ComputeHeatLake !----------------------------------------------------------------------- subroutine AccumulateHeatLake(bounds, num_c, filter_c, & temperature_inst, lakestate_inst, & - heat, heat_liquid, cv_liquid) + heat) ! ! !DESCRIPTION: ! Accumulate heat of lake water in lake columns. ! ! Adds to any existing values in heat, heat_liquid and cv_liquid. + ! HERE EXPLANATIONS about not adding heat_liquid and cv_liquid should come. + ! where this routine differs from AccumulateLiqIceMassLake !! ! ! Note: Changes to this routine should generally be accompanied by similar changes to ! AccumulateLiqIceMassLake @@ -1020,9 +1022,7 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & 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] - real(r8) , intent(inout) :: heat_liquid( bounds%begc: ) ! accumulated heat content: liquid water, excluding latent heat [J/m^2] - real(r8) , intent(inout) :: cv_liquid( bounds%begc: ) ! accumulated liquid water heat capacity [J/(m^2 K)] - ! + ! !LOCAL VARIABLES: integer :: fc integer :: l, c, j @@ -1031,14 +1031,12 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & 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] - + real(r8) :: cv_liquid( bounds%begc:bounds%endc ) ! dummy argument to use AccumulateLiquidWaterHeat - not used character(len=*), parameter :: subname = 'AccumulateHeatLake' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(heat) == [bounds%endc]), errMsg(sourcefile, __LINE__)) - SHR_ASSERT_ALL((ubound(heat_liquid) == [bounds%endc]), errMsg(sourcefile, __LINE__)) - SHR_ASSERT_ALL((ubound(cv_liquid) == [bounds%endc]), errMsg(sourcefile, __LINE__)) associate( & dz_lake => col%dz_lake, & ! lake layer depth (m) @@ -1077,7 +1075,6 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & ! add ice heat and liquid heat together (look at this part) do fc = 1, num_c c = filter_c(fc) - heat_liquid(c) = heat_liquid(c) + lake_heat_liquid(c) heat(c) = heat(c) + lake_heat_ice(c) + & lake_heat_liquid(c) + lake_latent_heat_liquid(c) end do @@ -1116,6 +1113,8 @@ subroutine AdjustDeltaHeatForDeltaLiq(bounds, delta_liq, & ! liquid and ice runoff in CLM, then this routine should be reworked to use the true ! heat contents of both liquid and ice runoff. ! + ! ADD HERE NOTES ABOUT LAKES!!! lake water not accounted. because baselines + ! ! 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 b30969689f..e14dcc1d17 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -251,9 +251,7 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & 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] - real(r8) :: lake_heat_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface - real(r8) :: lake_cv_liquid_col(bounds%begc:bounds%endc) ! unused; just needed for AccumulateSoilHeatNonLake interface - + character(len=*), parameter :: subname = 'dyn_heat_content_set_baselines' !----------------------------------------------------------------------- @@ -266,9 +264,7 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8 lake_heat_col(bounds%begc:bounds%endc) = 0._r8 - lake_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8 - lake_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, & @@ -299,7 +295,7 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & ! Set baselines for lake columns call AccumulateHeatLake(bounds, num_lakec, filter_lakec, & temperature_inst, lakestate_inst, & - heat = lake_heat_col, heat_liquid = lake_heat_liquid_col, cv_liquid = lake_cv_liquid_col) + heat = lake_heat_col) do fc = 1, num_lakec c = filter_lakec(fc) From 5dd892502b0ed3e24ae19568c140dbfb5897c220 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Fri, 20 Sep 2019 12:13:07 -0600 Subject: [PATCH 12/42] Added description for AccumulateLakeHeat --- src/biogeophys/TotalWaterAndHeatMod.F90 | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 5b5a7318ee..47b0a561d2 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -980,6 +980,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & ! 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 @@ -1008,13 +1009,21 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & ! !DESCRIPTION: ! Accumulate heat of lake water in lake columns. ! - ! Adds to any existing values in heat, heat_liquid and cv_liquid. - ! HERE EXPLANATIONS about not adding heat_liquid and cv_liquid should come. - ! where this routine differs from AccumulateLiqIceMassLake !! - ! - ! Note: Changes to this routine should generally be accompanied by similar changes to - ! AccumulateLiqIceMassLake - ! + ! This module differs from AccumulateSoilHeatNonLake in the sense that for lake heat + ! the average heat_liquid and cv_liquid are not calculated. This because the + ! lake water is virtual (there will never by a change in lake water content) + ! and thus 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. + ! + ! Some minor concerns with this approach: + ! In a rare case, 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) From 27d982ed1695182697eec43aacddd9abcf4e3c48 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 3 Oct 2019 03:49:20 -0600 Subject: [PATCH 13/42] Updated comments to apply to lakes and removed redundant file --- src/dyn_subgrid/dynlakeFileMod.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/dyn_subgrid/dynlakeFileMod.F90 b/src/dyn_subgrid/dynlakeFileMod.F90 index f3179b1f55..88f2176510 100644 --- a/src/dyn_subgrid/dynlakeFileMod.F90 +++ b/src/dyn_subgrid/dynlakeFileMod.F90 @@ -82,8 +82,8 @@ subroutine dynlake_init(bounds, dynlake_filename) ! read data PCT_LAKE ! - ! Note: if you want to change transient crops so that they are interpolated, rather - ! than jumping to each year's value on Jan 1 of that year, simply change wtcrop and + ! 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. @@ -103,11 +103,11 @@ subroutine dynlake_interp(bounds) ! !DESCRIPTION: ! Get lake cover for model time, when needed. ! - ! Sets col%wtlunit and lun%wtgcell for crop landunits. + ! Sets col%wtlunit and lun%wtgcell for lake landunits. ! - ! Note that crop cover currently jumps to its new value at the start of the year. + ! 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 wtcrop and wtcft dyn_var_time_interp_type variables rather than + ! simply by making wtlake and wtcft dyn_var_time_interp_type variables rather than ! dyn_var_time_uninterp_type. ! ! !USES: @@ -141,4 +141,4 @@ subroutine dynlake_interp(bounds) end subroutine dynlake_interp -end module dynlakeFileMod \ No newline at end of file +end module dynlakeFileMod From 69a655360ccafbbbbf23d680aad6a90cb6d889a0 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 3 Oct 2019 06:06:21 -0600 Subject: [PATCH 14/42] Updated comments in AccumulateHeatLake and some little fixes --- src/biogeophys/TotalWaterAndHeatMod.F90 | 56 ++++++++++++++---------- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 1 - src/dyn_subgrid/dynSubgridControlMod.F90 | 10 +---- src/dyn_subgrid/dynlakeFileMod.F90 | 2 +- 4 files changed, 36 insertions(+), 33 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 47b0a561d2..27d3d96a03 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -486,7 +486,7 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & ! Adds to any existing values in liquid_mass and ice_mass. ! ! Note: Changes to this routine should generally be accompanied by similar changes to - ! AccumulateSoilHeatLake + ! AccumulateHeatLake ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -497,9 +497,9 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & 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(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] - real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] + 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' !----------------------------------------------------------------------- @@ -517,11 +517,11 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & do fc = 1, num_c c = filter_c(fc) ! calculate lake liq and ice content per lake layer first - h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) - h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) + h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + h2olak_ice = dz_lake(c,j) * denice * lake_icefrac(c,j) - liquid_mass(c) = liquid_mass(c) + h2olak_liq(c,j) - ice_mass(c) = ice_mass(c) + h2olak_ice(c,j) + liquid_mass(c) = liquid_mass(c) + h2olak_liq + ice_mass(c) = ice_mass(c) + h2olak_ice end do end do @@ -726,6 +726,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) @@ -975,6 +980,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & 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) @@ -1009,16 +1015,24 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & ! !DESCRIPTION: ! Accumulate heat of lake water in lake columns. ! - ! This module differs from AccumulateSoilHeatNonLake in the sense that for lake heat - ! the average heat_liquid and cv_liquid are not calculated. This because the - ! lake water is virtual (there will never by a change in lake water content) - ! and thus it should not be taken into the average column temperature when adjusting the change in + ! 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. + ! 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 a rare case, lake water can have some changes in water content in time, + ! 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 @@ -1035,12 +1049,11 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & ! !LOCAL VARIABLES: integer :: fc integer :: l, c, j - real(r8) :: h2olak_liq(bounds%begc:bounds%endc,1:nlevlak) ! liquid water content of lake layer [kg/m²] - real(r8) :: h2olak_ice(bounds%begc:bounds%endc,1:nlevlak) ! ice water content of lake layer [kg/m²] + 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] - real(r8) :: cv_liquid( bounds%begc:bounds%endc ) ! dummy argument to use AccumulateLiquidWaterHeat - not used character(len=*), parameter :: subname = 'AccumulateHeatLake' !----------------------------------------------------------------------- @@ -1067,17 +1080,16 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & do fc = 1, num_c c = filter_c(fc) ! liquid heat - h2olak_liq(c,j) = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) call AccumulateLiquidWaterHeat( & temp = t_lake(c,j), & - h2o = h2olak_liq(c,j), & - cv_liquid = cv_liquid(c), & + h2o = h2olak_liq, & heat_liquid = lake_heat_liquid(c), & latent_heat_liquid = lake_latent_heat_liquid(c)) ! ice heat - h2olak_ice(c,j) = dz_lake(c,j) * denice * lake_icefrac(c,j) + h2olak_ice = dz_lake(c,j) * denice * lake_icefrac(c,j) lake_heat_ice(c) = lake_heat_ice(c) + & - TempToHeat(temp=t_lake(c,j), cv = (h2olak_ice(c,j) * cpice)) + TempToHeat(temp=t_lake(c,j), cv = (h2olak_ice * cpice)) end do end do diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index e14dcc1d17..954092b3e9 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -705,7 +705,6 @@ subroutine dyn_heat_content(bounds, & real(r8) :: heat_col(bounds%begc:bounds%endc) ! sum of heat content for all columns [J/m^2] real(r8) :: heat_liquid_col(bounds%begc:bounds%endc) ! sum of heat content for all columns: liquid water, excluding latent heat [J/m^2] real(r8) :: cv_liquid_col(bounds%begc:bounds%endc) ! sum of liquid heat capacity for all columns [J/(m^2 K)] - real(r8) :: heat_lake_col(bounds%begc:bounds%endc) ! heat content of lake column [J/m^2] real(r8) :: heat_liquid_grc(bounds%begg:bounds%endg) ! heat_liquid_col averaged to grid cell [J/m^2] real(r8) :: cv_liquid_grc(bounds%begg:bounds%endg) ! cv_liquid_col averaged to grid cell [J/(m^2 K)] diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index 01ff35ab33..4cd78c7c1f 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -287,15 +287,7 @@ 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 - - if (dyn_subgrid_control_inst%do_transient_lakes) then - if (use_fates) then - write(iulog,*) 'ERROR: do_transient_lakes currently does not work with use_fates' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end if - + end if end subroutine check_namelist_consistency diff --git a/src/dyn_subgrid/dynlakeFileMod.F90 b/src/dyn_subgrid/dynlakeFileMod.F90 index 88f2176510..304c191a02 100644 --- a/src/dyn_subgrid/dynlakeFileMod.F90 +++ b/src/dyn_subgrid/dynlakeFileMod.F90 @@ -36,7 +36,7 @@ module dynlakeFileMod ! Names of variables on file character(len=*), parameter :: lake_varname = 'PCT_LAKE' - ! TO DO: account for lake depth? + character(len=*), parameter, private :: sourcefile = & __FILE__ From 7f8ef4191d0c34e806e0a5f8f0a419b684fbc4ff Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Mon, 16 Mar 2020 04:08:21 -0600 Subject: [PATCH 15/42] Finetune option to run without dynamical lakes --- src/dyn_subgrid/dynSubgridDriverMod.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index 69bbb388a8..7edcec8f3b 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -129,8 +129,9 @@ subroutine dynSubgrid_init(bounds_proc, glc_behavior, crop_inst) ! Initialize stuff for prescribed transient lakes ! error if keyword argument like for dynamical land units above: dynlake_filename=get_flanduse_timeseries() - call dynlake_init(bounds_proc, dynlake_filename=get_flanduse_timeseries()) - + 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 From ddfba1e3d4d9bd6b9edf09ed40bd793d05c311b7 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Wed, 8 Jul 2020 06:59:21 -0600 Subject: [PATCH 16/42] Update argument lists of subroutines --- src/biogeophys/BalanceCheckMod.F90 | 6 ++++-- src/biogeophys/TotalWaterAndHeatMod.F90 | 7 ++++--- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 7 ++++--- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index 65844b7299..5375880f5c 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -152,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, & @@ -163,8 +164,8 @@ end subroutine BeginWaterBalance !----------------------------------------------------------------------- subroutine BeginWaterBalanceSingle(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & - soilhydrology_inst, waterstate_inst, waterdiagnostic_inst, waterbalance_inst, & - use_aquifer_layer) + soilhydrology_inst, lakestate_inst, waterstate_inst, & + waterdiagnostic_inst, waterbalance_inst, & use_aquifer_layer) ! ! !DESCRIPTION: ! Initialize column-level water balance at beginning of time step, for bulk or a @@ -177,6 +178,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 diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 27d3d96a03..e817cb6511 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -504,8 +504,8 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & character(len=*), parameter :: subname = 'AccumulateLiqIceMassLake' !----------------------------------------------------------------------- - SHR_ASSERT_ALL((ubound(liquid_mass) == [bounds%endc]), errMsg(sourcefile, __LINE__)) - SHR_ASSERT_ALL((ubound(ice_mass) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + 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 @@ -1058,7 +1058,8 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & character(len=*), parameter :: subname = 'AccumulateHeatLake' !----------------------------------------------------------------------- - SHR_ASSERT_ALL((ubound(heat) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL_FL((ubound(heat) == (/bounds%endc/)), sourcefile, __LINE__) + associate( & dz_lake => col%dz_lake, & ! lake layer depth (m) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 954092b3e9..9947203a3f 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -541,7 +541,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 @@ -557,6 +557,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: @@ -575,7 +576,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)) @@ -623,7 +624,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 + 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 ! From 83484a9e7d15a8b4cf02a734775f523d709b0f78 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 16 Jul 2020 19:54:01 -0600 Subject: [PATCH 17/42] Add line break to fix compilation error --- src/biogeophys/BalanceCheckMod.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index 5375880f5c..65f102a548 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -165,7 +165,8 @@ end subroutine BeginWaterBalance subroutine BeginWaterBalanceSingle(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & soilhydrology_inst, lakestate_inst, waterstate_inst, & - waterdiagnostic_inst, waterbalance_inst, & use_aquifer_layer) + waterdiagnostic_inst, waterbalance_inst, & + use_aquifer_layer) ! ! !DESCRIPTION: ! Initialize column-level water balance at beginning of time step, for bulk or a From 124c55262f72515b272a07df4ed3556b9078391e Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 22 Jul 2020 06:11:45 -0600 Subject: [PATCH 18/42] Fix data type of haslake variable This is a double precision variable on the file, so I'm making it double precision in the code. I think ncdio_pio used to handle data type conversions like this, but now it seems to be corrupting memory: the original code led to errors in decompInitMod, due to what appears to be corruption in some unrelated variable. --- src/main/surfrdMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index f73265cd92..af092daa6b 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -995,7 +995,7 @@ subroutine surfrd_lakemask(begg, endg, ncid) ! ! !LOCAL VARIABLES: logical :: readvar - integer ,pointer :: haslake_id(:) + real(r8),pointer :: haslake_id(:) ! ! character(len=*), parameter :: subname = 'surfrd_lakemask' @@ -1009,7 +1009,7 @@ subroutine surfrd_lakemask(begg, endg, ncid) 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.) + where (haslake_id == 1._r8) haslake = .true. elsewhere haslake = .false. From f9f87b61b35ac7e74cef51d2508ed080fe9148f0 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 12 Aug 2020 14:25:22 -0600 Subject: [PATCH 19/42] Add new runtime consistency checks for do_transient_lakes I'm treating do_transient_lakes similar to do_transient_pfts and do_transient_crops in this respect. 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. --- src/dyn_subgrid/dynSubgridControlMod.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index 4cd78c7c1f..a2ee0ef17d 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -248,9 +248,11 @@ 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 + 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 @@ -258,7 +260,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, & From e18234149300c7051b027d6e110352a4ecf52971 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 12 Aug 2020 14:45:54 -0600 Subject: [PATCH 20/42] do_transient_lakes in namelist: set default value and do error checks The new function setup_logic_do_transient_lakes is based on setup_logic_do_transient_crops. As with my previous commit: 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. --- bld/CLMBuildNamelist.pm | 69 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 4ee6a65eda..de106915ec 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2372,6 +2372,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'); @@ -2537,6 +2538,74 @@ 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 + # + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + + my $var = 'do_transient_lakes'; + + # Start by assuming a default value of '.true.'. Then check a number of + # conditions under which do_transient_lakes cannot be true. Under these + # conditions: (1) set default value to '.false.'; (2) make sure that the + # value is indeed false (e.g., that the user didn't try to set it to true). + + my $default_val = ".true."; + + # 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) { + $default_val = ".false."; + } + + 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, val=>$default_val); + } + + # 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 From 5a64ed84bd9a9a777d4bbe848292e3eafd3ed36b Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 12 Aug 2020 20:57:56 -0600 Subject: [PATCH 21/42] Get do_transient_lakes default from namelist_defaults I didn't see any reason why the default value should be set in the code rather than coming from namelist_defaults_ctsm.xml. Getting it from the latter will make it easier for the default value to differ with different physics versions. --- bld/CLMBuildNamelist.pm | 13 +------------ bld/namelist_files/namelist_defaults_ctsm.xml | 1 + 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index de106915ec..5d79b50ddd 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2550,13 +2550,6 @@ sub setup_logic_do_transient_lakes { my $var = 'do_transient_lakes'; - # Start by assuming a default value of '.true.'. Then check a number of - # conditions under which do_transient_lakes cannot be true. Under these - # conditions: (1) set default value to '.false.'; (2) make sure that the - # value is indeed false (e.g., that the user didn't try to set it to true). - - my $default_val = ".true."; - # 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 @@ -2576,15 +2569,11 @@ sub setup_logic_do_transient_lakes { $cannot_be_true = "$var can only be set to true when running a transient case (flanduse_timeseries non-blank)"; } - if ($cannot_be_true) { - $default_val = ".false."; - } - 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, val=>$default_val); + 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 diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index e9dae4e762..5918e30cee 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2965,6 +2965,7 @@ lnd/clm2/surfdata_map/landuse.timeseries_ne30np4_hist_16pfts_Irrig_CMIP6_simyr18 +.false. .false. From 273b75ee959427fde57d4a9af9abbdd3f4796457 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sun, 23 Aug 2020 06:35:33 -0600 Subject: [PATCH 22/42] Add some comments --- bld/CLMBuildNamelist.pm | 6 ++++++ src/dyn_subgrid/dynSubgridControlMod.F90 | 7 +++++++ 2 files changed, 13 insertions(+) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 5d79b50ddd..7a3d6b4193 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2546,6 +2546,12 @@ sub setup_logic_do_transient_lakes { # 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'; diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index a2ee0ef17d..6d9a703b31 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -248,6 +248,13 @@ subroutine check_namelist_consistency end if end if + ! 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 From 9f8d2df69ec58051cc0fb464309cbcfd9938988c Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sun, 23 Aug 2020 06:51:23 -0600 Subject: [PATCH 23/42] Remove unnecessary use statements and unnecessary comments --- src/dyn_subgrid/dynSubgridDriverMod.F90 | 7 +------ src/dyn_subgrid/dynlakeFileMod.F90 | 19 ++++--------------- 2 files changed, 5 insertions(+), 21 deletions(-) diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index 7edcec8f3b..676e14da81 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -20,7 +20,6 @@ module dynSubgridDriverMod use dynpftFileMod , only : dynpft_init, dynpft_interp use dyncropFileMod , only : dyncrop_init, dyncrop_interp use dynHarvestMod , only : dynHarvest_init, dynHarvest_interp - ! use new lake file module use dynlakeFileMod , only : dynlake_init, dynlake_interp use dynLandunitAreaMod , only : update_landunit_weights use subgridWeightsMod , only : compute_higher_order_weights, set_subgrid_diagnostic_fields @@ -44,10 +43,7 @@ module dynSubgridDriverMod use CropType , only : crop_type use glc2lndMod , only : glc2lnd_type use filterMod , only : filter, filter_inactive_and_active - - - - + ! ! !PUBLIC MEMBER FUNCTIONS: implicit none private @@ -128,7 +124,6 @@ subroutine dynSubgrid_init(bounds_proc, glc_behavior, crop_inst) ! Initialize stuff for prescribed transient lakes - ! error if keyword argument like for dynamical land units above: dynlake_filename=get_flanduse_timeseries() if (get_do_transient_lakes()) then call dynlake_init(bounds_proc, dynlake_filename=get_flanduse_timeseries()) end if diff --git a/src/dyn_subgrid/dynlakeFileMod.F90 b/src/dyn_subgrid/dynlakeFileMod.F90 index 304c191a02..17b21484fd 100644 --- a/src/dyn_subgrid/dynlakeFileMod.F90 +++ b/src/dyn_subgrid/dynlakeFileMod.F90 @@ -6,8 +6,6 @@ module dynlakeFileMod ! ! !USES: - ! check this still on unneccesary parts. - #include "shr_assert.h" use shr_log_mod , only : errMsg => shr_log_errMsg use shr_kind_mod , only : r8 => shr_kind_r8 @@ -15,12 +13,9 @@ module dynlakeFileMod use dynFileMod , only : dyn_file_type use dynVarTimeUninterpMod , only : dyn_var_time_uninterp_type use clm_varctl , only : iulog - use clm_varcon , only : grlnd, namec + use clm_varcon , only : grlnd use abortutils , only : endrun - use spmdMod , only : masterproc, mpicom - use LandunitType , only : lun - use ColumnType , only : col - use PatchType , only : patch + use spmdMod , only : masterproc ! !PUBLIC MEMBER FUNCTIONS: implicit none @@ -35,8 +30,6 @@ module dynlakeFileMod ! Names of variables on file character(len=*), parameter :: lake_varname = 'PCT_LAKE' - - character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -75,8 +68,8 @@ subroutine dynlake_init(bounds, dynlake_filename) ! 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 crop updates determined - ! prognostically, if crop areas are ever determined prognostically rather than + ! 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) @@ -92,7 +85,6 @@ subroutine dynlake_init(bounds, dynlake_filename) 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 @@ -136,9 +128,6 @@ subroutine dynlake_interp(bounds) end do deallocate(wtlake_cur) -! also account for landunits, patches and columns? - - end subroutine dynlake_interp end module dynlakeFileMod From 9976b45b44cb716cb0b575888f0af9688cf76c2e Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sun, 23 Aug 2020 07:31:14 -0600 Subject: [PATCH 24/42] Whitespace cleanup --- src/biogeophys/TotalWaterAndHeatMod.F90 | 14 ++++++-------- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 6 ++---- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index e817cb6511..e5f5b17fa7 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -881,9 +881,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & ! ! Note: Changes to this routine - for anything involving liquid or ice - should ! generally be accompanied by similar changes to ComputeLiqIceMassLake - - - + ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_lakec @@ -974,13 +972,14 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & end do end do - do fc = 1, num_lakec + 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) 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. + ! 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) @@ -1079,8 +1078,8 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & ! calculate heat content of lake itself do j = 1, nlevlak do fc = 1, num_c - c = filter_c(fc) - ! liquid heat + 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), & @@ -1265,7 +1264,6 @@ subroutine AccumulateLiquidWaterHeat(temp, h2o, & latent_heat_liquid = latent_heat_liquid + h2o*hfus end subroutine AccumulateLiquidWaterHeat - !----------------------------------------------------------------------- pure function TempToHeat(temp, cv) result(heat) ! diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 9947203a3f..d8bb809bba 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -203,17 +203,15 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & ! set baselines for lake columns - ! Calculate the total water volume of the lake column + ! Calculate the total water volume of the lake column call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, & 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) - + 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 associate From e44c19d5959d7d1cdfe82b7047fcca263a6e0280 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Mon, 24 Aug 2020 06:16:59 -0600 Subject: [PATCH 25/42] Remove temperature_inst%lake_heat as not used, remove comments on saving lake heat in LakeTemperatureMod --- src/biogeophys/LakeTemperatureMod.F90 | 6 +----- src/biogeophys/TemperatureType.F90 | 3 --- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/src/biogeophys/LakeTemperatureMod.F90 b/src/biogeophys/LakeTemperatureMod.F90 index 6a20856f67..280020b6ef 100644 --- a/src/biogeophys/LakeTemperatureMod.F90 +++ b/src/biogeophys/LakeTemperatureMod.F90 @@ -248,7 +248,6 @@ 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) - lake_heat => temperature_inst%lake_heat , & ! Output: [real(r8) (:) ] col lake heat (J/m²) 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 @@ -1004,10 +1003,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la end do end do - ! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts - - ! IV: currently commented out: caused crash. To do: look at this part of the code!!! - ! lake_heat(c) = ncvts(c) + call waterstatebulk_inst%CalculateTotalH2osno(bounds, num_lakec, filter_lakec, & caller = 'LakeTemperature-2', & diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 5dacf3fa55..b165ed0796 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -116,9 +116,6 @@ 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 - - ! lake heat - real(r8), pointer :: lake_heat (:) ! total heat of lake water (J/m²) contains From 2bfb8d7d39cd07469919dd04c402df5cb3e43a82 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Mon, 24 Aug 2020 07:36:17 -0600 Subject: [PATCH 26/42] Add description on dynamical lakes for AdjustDeltaHeatForDeltaLiq module and clean up --- src/biogeophys/TotalWaterAndHeatMod.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index e5f5b17fa7..6e4ae8624d 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -1093,7 +1093,7 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & end do end do - ! add ice heat and liquid heat together (look at this part) + ! add ice heat and liquid heat together do fc = 1, num_c c = filter_c(fc) heat(c) = heat(c) + lake_heat_ice(c) + & @@ -1130,11 +1130,15 @@ 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. ! - ! ADD HERE NOTES ABOUT LAKES!!! lake water not accounted. because baselines ! ! Sign convention: delta_liq and delta_heat are positive if the post-landcover change ! value is greater than the pre-landcover change value. From e37d269fe0e98e1aa0fd223196ea02e87fc59328 Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Mon, 24 Aug 2020 08:49:02 -0600 Subject: [PATCH 27/42] remove lake_heat allocation and history field --- src/biogeophys/TemperatureType.F90 | 7 ------- src/biogeophys/TotalWaterAndHeatMod.F90 | 1 - 2 files changed, 8 deletions(-) diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index b165ed0796..d9ef9f5cb6 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -279,7 +279,6 @@ subroutine InitAllocate(this, bounds) allocate(this%fact_col (begc:endc, -nlevsno+1:nlevgrnd)) ; this%fact_col (:,:) = nan allocate(this%c_h2osfc_col (begc:endc)) ; this%c_h2osfc_col (:) = nan - allocate(this%lake_heat (begc:endc)) ; this%lake_heat (:) = nan end subroutine InitAllocate !------------------------------------------------------------------------ @@ -621,12 +620,6 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) ptr_patch=this%t_veg10_night_patch, default='inactive') endif - ! add lake heat history field here - this%lake_heat(begc:endc) = spval - call hist_addfld1d (fname='LAKE_HEAT', units='J/m^2', & - avgflag='A', long_name='Heat content of gridcell lake water', & - ptr_col=this%lake_heat, default='active') - end subroutine InitHistory !----------------------------------------------------------------------- diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 6e4ae8624d..b774fe9761 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -919,7 +919,6 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & csol => soilstate_inst%csol_col, & ! heat capacity, soil solids (J/m**3/Kelvin) 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) - lake_heat => temperature_inst%lake_heat, & ! total heat of lake water (J/m²) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col, & ! liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! frozen water (kg/m2) ) From 30394f88b806ddc6f349572efc45e5fe000a83df Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Mon, 24 Aug 2020 10:09:06 -0600 Subject: [PATCH 28/42] Move code to read lakemask to surfd_lakemask routine and clean up --- src/main/surfrdMod.F90 | 75 +++++++++++++++++++++--------------------- 1 file changed, 38 insertions(+), 37 deletions(-) diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index af092daa6b..b9d2e1bb3d 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -292,9 +292,10 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) 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 landunit_varcon , only : max_lunit, istsoil, isturb_MIN, isturb_MAX use dynSubgridControlMod, only : get_flanduse_timeseries - use clm_varctl , only : fname_len + use dynSubgridControlMod, only : get_do_transient_lakes + ! ! !ARGUMENTS: integer, intent(in) :: begg, endg, actual_numcft @@ -312,10 +313,8 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) logical :: readvar ! true => variable is on dataset real(r8) :: rmaxlon,rmaxlat ! local min/max vars type(file_desc_t) :: ncid ! netcdf id - type(file_desc_t) :: ncid_dynuse ! netcdf id for landuse timeseries file logical :: istype_domain ! true => input file is of type domain logical :: isgrid2d ! true => intut grid is 2d - character(len=fname_len) :: fdynuse ! landuse.timeseries filename character(len=32) :: subname = 'surfrd_get_data' ! subroutine name !----------------------------------------------------------------------- @@ -450,33 +449,11 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) 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 - ! IV: add here call to subroutine to read in lake mask (necessary for initialisation of dynamical lakes) - ! First open landuse.timeseries file for this. - - if (masterproc) then - write(iulog,*) 'Attempting to read landuse.timeseries data .....' - if (lfsurdat == ' ') then - write(iulog,*)'fdynuse must be specified' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - endif - - - ! open landuse_timeseries file - fdynuse = get_flanduse_timeseries() - - call getfil(fdynuse, locfn, 0 ) - - call ncd_pio_openfile (ncid_dynuse, trim(locfn), 0) - - ! read the lakemask - call surfrd_lakemask(begg, endg, ncid_dynuse) - - ! close landuse_timeseries file again - call ncd_pio_closefile(ncid_dynuse) - - end subroutine surfrd_get_data !----------------------------------------------------------------------- @@ -979,7 +956,7 @@ subroutine surfrd_veg_dgvm(begg, endg) end subroutine surfrd_veg_dgvm !----------------------------------------------------------------------- - subroutine surfrd_lakemask(begg, endg, ncid) + subroutine surfrd_lakemask(begg, endg) ! ! !DESCRIPTION: ! Reads the lake mask, indicating where lakes are and will grow @@ -987,25 +964,47 @@ subroutine surfrd_lakemask(begg, endg, ncid) ! Necessary for the initialisation of the lake land units ! ! !USES: - use clm_instur , only : haslake + 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 - type(file_desc_t), intent(inout) :: ncid ! ! ! !LOCAL VARIABLES: - logical :: readvar - real(r8),pointer :: haslake_id(:) - ! + 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)) - call ncd_io(ncid=ncid, varname='HASLAKE' , flag='read', data=haslake_id, & + ! 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__)) @@ -1015,6 +1014,8 @@ subroutine surfrd_lakemask(begg, endg, ncid) haslake = .false. end where + ! close landuse_timeseries file again + call ncd_pio_closefile(ncid_dynuse) deallocate(haslake_id) From a6726928f50f10268d9659130fd4f349d9d1a34c Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Thu, 27 Aug 2020 13:50:41 -0600 Subject: [PATCH 29/42] Adjust ice to liquid density in ice mass calculation, as lake layer is not adjusted --- src/biogeophys/TotalWaterAndHeatMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index b774fe9761..87d3256bbf 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -518,7 +518,7 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & c = filter_c(fc) ! calculate lake liq and ice content per lake layer first h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) - h2olak_ice = dz_lake(c,j) * denice * lake_icefrac(c,j) + h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) ! use water density of liquid water as layer depth is not adjusted liquid_mass(c) = liquid_mass(c) + h2olak_liq ice_mass(c) = ice_mass(c) + h2olak_ice From ffc879f88203767e8fa4d70d69e0b192d73deb8d Mon Sep 17 00:00:00 2001 From: ivanderkelen Date: Fri, 28 Aug 2020 00:50:23 -0600 Subject: [PATCH 30/42] Adjust ice density to liquid density in AccumulateHeatLake routine --- src/biogeophys/TotalWaterAndHeatMod.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 87d3256bbf..79ee0a77d6 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -1086,7 +1086,8 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & heat_liquid = lake_heat_liquid(c), & latent_heat_liquid = lake_latent_heat_liquid(c)) ! ice heat - h2olak_ice = dz_lake(c,j) * denice * lake_icefrac(c,j) + ! use water density as lake layer does not adjust + 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 From 82f39d74448d03d1fccf5ea8e383c7624a03a253 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2020 08:51:15 -0600 Subject: [PATCH 31/42] Expand on some comments a bit --- src/biogeophys/TotalWaterAndHeatMod.F90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 79ee0a77d6..fcf81941fc 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -518,7 +518,10 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & c = filter_c(fc) ! calculate lake liq and ice content per lake layer first h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) - h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) ! use water density of liquid water as layer depth is not adjusted + + ! 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) liquid_mass(c) = liquid_mass(c) + h2olak_liq ice_mass(c) = ice_mass(c) + h2olak_ice @@ -1085,10 +1088,11 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & h2o = h2olak_liq, & heat_liquid = lake_heat_liquid(c), & latent_heat_liquid = lake_latent_heat_liquid(c)) - ! ice heat - ! use water density as lake layer does not adjust - h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) - lake_heat_ice(c) = lake_heat_ice(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 From 5658ddec583cd6183d4b282906c8062e774b2171 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2020 11:00:03 -0600 Subject: [PATCH 32/42] Tweak some comments --- src/biogeophys/TotalWaterAndHeatMod.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index fcf81941fc..3bec595bf0 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -107,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) @@ -159,8 +159,8 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & 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. + ! 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) @@ -215,8 +215,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) @@ -408,8 +408,8 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & 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. + ! 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) From a9fa8754cf01e39a129ba62369358e2e10219e48 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2020 11:10:14 -0600 Subject: [PATCH 33/42] Do not add lake water to begwb and endwb, and thus also TWS See https://github.com/ESCOMP/CTSM/issues/659#issuecomment-682867045 for details. --- src/biogeophys/BalanceCheckMod.F90 | 2 +- src/biogeophys/LakeHydrologyMod.F90 | 2 +- src/biogeophys/TotalWaterAndHeatMod.F90 | 30 +++++++++++++++--------- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 2 +- 4 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index 65f102a548..69f7cdb167 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -216,7 +216,7 @@ subroutine BeginWaterBalanceSingle(bounds, & call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & waterstate_inst, lakestate_inst, & - subtract_dynbal_baselines = .false., & + 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/LakeHydrologyMod.F90 b/src/biogeophys/LakeHydrologyMod.F90 index c078e958e1..f6f83d8956 100644 --- a/src/biogeophys/LakeHydrologyMod.F90 +++ b/src/biogeophys/LakeHydrologyMod.F90 @@ -650,7 +650,7 @@ subroutine LakeHydrology(bounds, & call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & b_waterstate_inst, lakestate_inst, & - subtract_dynbal_baselines = .false., & + add_lake_water_and_subtract_dynbal_baselines = .false., & water_mass = endwb(bounds%begc:bounds%endc)) do j = 1, nlevgrnd diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 3bec595bf0..9cb1d948d6 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -143,7 +143,7 @@ end subroutine ComputeWaterMassNonLake !----------------------------------------------------------------------- subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & waterstate_inst, lakestate_inst, & - subtract_dynbal_baselines, & + add_lake_water_and_subtract_dynbal_baselines, & water_mass) ! ! !DESCRIPTION: @@ -159,9 +159,12 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & class(waterstate_type) , intent(in) :: waterstate_inst type(lakestate_type) , intent(in) :: lakestate_inst + ! 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) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass + logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines real(r8) , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2) ! @@ -181,7 +184,7 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & filter_lakec = filter_lakec, & waterstate_inst = waterstate_inst, & lakestate_inst = lakestate_inst, & - subtract_dynbal_baselines = subtract_dynbal_baselines, & + 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)) @@ -388,7 +391,7 @@ end subroutine AccumulateSoilLiqIceMassNonLake !----------------------------------------------------------------------- subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & waterstate_inst, lakestate_inst,& - subtract_dynbal_baselines, & + add_lake_water_and_subtract_dynbal_baselines, & liquid_mass, ice_mass) ! ! !DESCRIPTION: @@ -408,9 +411,12 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & type(lakestate_type) , intent(in) :: lakestate_inst + ! 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) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass + 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) @@ -449,11 +455,13 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) end do end do - - ! Lake water content - call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, & - lakestate_inst, liquid_mass, ice_mass) - + + if (add_lake_water_and_subtract_dynbal_baselines) then + ! Lake water content + call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, & + lakestate_inst, liquid_mass, ice_mass) + end if + ! Soil water content of the soil under the lake do j = 1, nlevgrnd do fc = 1, num_lakec @@ -463,7 +471,7 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & end do end do - if (subtract_dynbal_baselines) then + if (add_lake_water_and_subtract_dynbal_baselines) then ! Subtract baselines set in initialization do fc = 1, num_lakec c = filter_lakec(fc) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index d8bb809bba..9dc9b46002 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -645,7 +645,7 @@ subroutine dyn_water_content(bounds, & call ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & waterstate_inst, & lakestate_inst, & - subtract_dynbal_baselines = .true., & + 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)) From 52105c4cd824764e1e2f862eacd7ace8696a76f3 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 1 Sep 2020 16:28:48 -0600 Subject: [PATCH 34/42] Use tracer ratio when counting lake water This is needed for water tracer masses to be counted correctly --- src/biogeophys/TotalWaterAndHeatMod.F90 | 19 ++++++++++++------- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 5 +++-- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 9cb1d948d6..ba794f67a4 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -390,7 +390,7 @@ end subroutine AccumulateSoilLiqIceMassNonLake !----------------------------------------------------------------------- subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & - waterstate_inst, lakestate_inst,& + waterstate_inst, lakestate_inst, & add_lake_water_and_subtract_dynbal_baselines, & liquid_mass, ice_mass) ! @@ -410,7 +410,6 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & class(waterstate_type), intent(in) :: waterstate_inst type(lakestate_type) , intent(in) :: lakestate_inst - ! 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 ! @@ -459,7 +458,10 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & if (add_lake_water_and_subtract_dynbal_baselines) then ! Lake water content call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, & - lakestate_inst, liquid_mass, ice_mass) + 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 ! Soil water content of the soil under the lake @@ -486,7 +488,7 @@ end subroutine ComputeLiqIceMassLake !----------------------------------------------------------------------- subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & - lakestate_inst, liquid_mass, ice_mass) + lakestate_inst, tracer_ratio, liquid_mass, ice_mass) ! ! !DESCRIPTION: ! Accumulate lake water mass of lake columns, separated into liquid and ice. @@ -501,6 +503,7 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & 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) ! @@ -524,12 +527,14 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, & do j = 1, nlevlak do fc = 1, num_c c = filter_c(fc) - ! calculate lake liq and ice content per lake layer first - h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) + ! 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) + 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 diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 9dc9b46002..6e42b405eb 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -205,8 +205,9 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & ! Calculate the total water volume of the lake column call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, & - liquid_mass = lake_liquid_mass_col(bounds%begc:bounds%endc), & - ice_mass = lake_ice_mass_col(bounds%begc:bounds%endc)) + 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) From de3e12c83ebfc2b825908d1e9af6da3cbae4d421 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2020 11:11:37 -0600 Subject: [PATCH 35/42] Reset lake dynbal baselines if using an old initial conditions file We have changed the definition of total column water for lake columns, so the baseline values for lakes are incorrect on old initial conditions files. This commit adds some code to check if we're using an old initial conditions file, and if so, resets the dynbal baseline values for lakes to use the new definition. --- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 251 ++++++++++++++--------- src/main/clm_initializeMod.F90 | 66 +++--- src/main/restFileMod.F90 | 78 ++++++- 3 files changed, 269 insertions(+), 126 deletions(-) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 6e42b405eb..e3b25b1322 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -12,6 +12,8 @@ module dynConsBiogeophysMod use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type + use spmdMod , only : masterproc + use clm_varctl , only : iulog use UrbanParamsType , only : urbanparams_type use EnergyFluxType , only : energyflux_type use SoilHydrologyType , only : soilhydrology_type @@ -66,7 +68,8 @@ module dynConsBiogeophysMod !----------------------------------------------------------------------- subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & num_lakec, filter_lakec, & - urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst) + urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, & + reset_all_baselines, reset_lake_baselines, print_resets) ! ! !DESCRIPTION: ! Set start-of-run baseline values for heat and water content in some columns. @@ -79,13 +82,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 @@ -98,12 +106,28 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & 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 + + ! Whether to print information about the reset flags + logical, intent(in) :: print_resets + ! ! !LOCAL VARIABLES: integer :: i @@ -112,6 +136,18 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & character(len=*), parameter :: subname = 'dyn_hwcontent_set_baselines' !----------------------------------------------------------------------- + if (masterproc .and. print_resets) then + if (reset_all_baselines) then + write(iulog,*) ' ' + write(iulog,*) 'Resetting dynbal baselines for all columns' + write(iulog,*) ' ' + else if (reset_lake_baselines) then + write(iulog,*) ' ' + write(iulog,*) 'Resetting dynbal baselines for lake columns' + write(iulog,*) ' ' + end if + end if + ! Note that we include inactive points in the following. This could be important if ! an inactive point later becomes active, so that we have an appropriate baseline ! value for that point. @@ -125,22 +161,26 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & call dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & num_icemecc, filter_icemecc, num_lakec, filter_lakec, & - bulk_or_tracer%waterstate_inst, lakestate_inst) + 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, num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, lakestate_inst, water_inst%waterstatebulk_inst, & - temperature_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, num_lakec, filter_lakec, & - waterstate_inst, lakestate_inst) + 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 @@ -155,9 +195,13 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & 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 + 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 @@ -174,47 +218,50 @@ 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 - lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8 - lake_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)) - - - ! set baselines for lake columns - - ! 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 + 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 + + 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 end subroutine dyn_water_content_set_baselines @@ -223,7 +270,8 @@ end subroutine dyn_water_content_set_baselines subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & num_icemecc, filter_icemecc, num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, lakestate_inst, waterstatebulk_inst, & - temperature_inst) + temperature_inst, & + reset_all_baselines, reset_lake_baselines) ! ! !DESCRIPTION: ! Set start-of-run baseline values for heat content. @@ -243,6 +291,10 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & 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: integer :: c, fc ! indices @@ -258,51 +310,52 @@ 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 - - lake_heat_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)) - - - ! Set baselines for lake columns - call AccumulateHeatLake(bounds, num_lakec, filter_lakec, & - temperature_inst, lakestate_inst, & - heat = lake_heat_col) - - do fc = 1, num_lakec - c = filter_lakec(fc) + 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 - dynbal_baseline_heat(c) = lake_heat_col(c) - - end do - end associate end subroutine dyn_heat_content_set_baselines diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index c07a47c7ed..b8f5d0aa95 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -339,6 +339,7 @@ subroutine initialize2( ) logical :: lexist integer :: closelatidx,closelonidx real(r8) :: closelat,closelon + logical :: reset_dynbal_baselines_lake_columns integer :: begp, endp integer :: begc, endc integer :: begl, endl @@ -490,7 +491,13 @@ subroutine initialize2( ) 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) + 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., & + ! no need to print information about these resets for this initial resetting, + ! which is always done + print_resets = .false.) end do ! ------------------------------------------------------------------------ @@ -543,6 +550,7 @@ subroutine initialize2( ) is_cold_start = .false. is_interpolated_start = .false. + reset_dynbal_baselines_lake_columns = .false. if (nsrest == nsrStartup) then @@ -563,7 +571,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 @@ -571,7 +580,13 @@ 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) + ! But for a continue or branch run, it seems safest to NOT reset lake dynbal + ! baselines. In nearly all cases, this will be irrelevant, because the restart + ! file will have been a recent one, where reset_dynbal_baselines_lake_columns + ! should be false already. + reset_dynbal_baselines_lake_columns = .false. end if @@ -598,7 +613,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) @@ -613,33 +629,31 @@ 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 + if (nsrest == nsrStartup) then - !$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) - end do - else if (nsrest == nsrBranch) then + !$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 = get_reset_dynbal_baselines(), & + reset_lake_baselines = reset_dynbal_baselines_lake_columns, & + print_resets = .true.) + end do + else if (nsrest == nsrBranch) then + if (get_reset_dynbal_baselines()) 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 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 ! ------------------------------------------------------------------------ ! Initialize nitrogen deposition diff --git a/src/main/restFileMod.F90 b/src/main/restFileMod.F90 index 83be13835b..4302cbc1c0 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) ! From acf09849ccbf6e1f3f4eade1ad959a7fb055ee0b Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2020 12:05:49 -0600 Subject: [PATCH 36/42] Reset lake baselines in branch or continue run, too I went back and forth about whether we should do this, but I actually feel that it's best if we do reset the lake baselines in a branch or continue run, if using an older restart file. If we didn't do this, we'd want to add some logic for writing out the issue-fixed metadata for any further restart files written from these runs, to note that this issue isn't actually fixed yet on these restart files. --- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 19 +------ src/main/clm_initializeMod.F90 | 71 +++++++++++++----------- 2 files changed, 41 insertions(+), 49 deletions(-) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index e3b25b1322..089715233d 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -12,8 +12,6 @@ module dynConsBiogeophysMod use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type - use spmdMod , only : masterproc - use clm_varctl , only : iulog use UrbanParamsType , only : urbanparams_type use EnergyFluxType , only : energyflux_type use SoilHydrologyType , only : soilhydrology_type @@ -69,7 +67,7 @@ module dynConsBiogeophysMod subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, & - reset_all_baselines, reset_lake_baselines, print_resets) + reset_all_baselines, reset_lake_baselines) ! ! !DESCRIPTION: ! Set start-of-run baseline values for heat and water content in some columns. @@ -125,9 +123,6 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & ! for lakes on older initial conditions files. logical, intent(in) :: reset_lake_baselines - ! Whether to print information about the reset flags - logical, intent(in) :: print_resets - ! ! !LOCAL VARIABLES: integer :: i @@ -136,18 +131,6 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, & character(len=*), parameter :: subname = 'dyn_hwcontent_set_baselines' !----------------------------------------------------------------------- - if (masterproc .and. print_resets) then - if (reset_all_baselines) then - write(iulog,*) ' ' - write(iulog,*) 'Resetting dynbal baselines for all columns' - write(iulog,*) ' ' - else if (reset_lake_baselines) then - write(iulog,*) ' ' - write(iulog,*) 'Resetting dynbal baselines for lake columns' - write(iulog,*) ' ' - end if - end if - ! Note that we include inactive points in the following. This could be important if ! an inactive point later becomes active, so that we have an appropriate baseline ! value for that point. diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index b8f5d0aa95..771b88e5da 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -339,6 +339,7 @@ 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 @@ -494,10 +495,7 @@ subroutine initialize2( ) 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., & - ! no need to print information about these resets for this initial resetting, - ! which is always done - print_resets = .false.) + reset_lake_baselines = .false.) end do ! ------------------------------------------------------------------------ @@ -582,12 +580,6 @@ subroutine initialize2( ) end if call restFile_read(bounds_proc, fnamer, glc_behavior, & reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns) - ! But for a continue or branch run, it seems safest to NOT reset lake dynbal - ! baselines. In nearly all cases, this will be irrelevant, because the restart - ! file will have been a recent one, where reset_dynbal_baselines_lake_columns - ! should be false already. - reset_dynbal_baselines_lake_columns = .false. - end if ! ------------------------------------------------------------------------ @@ -629,31 +621,48 @@ subroutine initialize2( ) ! interpolated restart file, if applicable). ! ------------------------------------------------------------------------ - if (nsrest == nsrStartup) then - - !$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 = get_reset_dynbal_baselines(), & - reset_lake_baselines = reset_dynbal_baselines_lake_columns, & - print_resets = .true.) - end do - else if (nsrest == nsrBranch) then - if (get_reset_dynbal_baselines()) 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 + 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 - ! nsrContinue not explicitly handled: it's okay for reset_dynbal_baselines to - ! remain set in a continue run, but it has no effect + + !$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 ! ------------------------------------------------------------------------ ! Initialize nitrogen deposition From 8088c3c5ae7c0f391fa67ddf1367727eaec5ddbd Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2020 13:44:54 -0600 Subject: [PATCH 37/42] Add missing 'OMP END PARALLEL DO' statements These seem to have been missing for a while (forever?). --- src/main/clm_initializeMod.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 771b88e5da..912bbb293b 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -497,6 +497,7 @@ subroutine initialize2( ) ! 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) @@ -663,6 +664,7 @@ subroutine initialize2( ) 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 From a31875d11287bc40c2645198ab39ca39234ec914 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2020 14:19:00 -0600 Subject: [PATCH 38/42] Reorder accumulations to avoid loss of precision --- src/biogeophys/TotalWaterAndHeatMod.F90 | 117 ++++++++++++++---------- 1 file changed, 69 insertions(+), 48 deletions(-) diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index ba794f67a4..cd7316af81 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -444,18 +444,22 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & ice_mass(c) = 0._r8 end do - ! Snow water content - do fc = 1, num_lakec - c = filter_lakec(fc) + ! ------------------------------------------------------------------------ + ! 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.) + ! ------------------------------------------------------------------------ - ice_mass(c) = ice_mass(c) + h2osno_no_layers(c) - do j = snl(c)+1,0 - liquid_mass(c) = liquid_mass(c) + h2osoi_liq(c,j) - ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) + 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 - end do - if (add_lake_water_and_subtract_dynbal_baselines) then ! Lake water content call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, & lakestate_inst, & @@ -464,23 +468,29 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & ice_mass = ice_mass(bounds%begc:bounds%endc)) end if - ! Soil water content of the soil under the lake - do j = 1, nlevgrnd - do fc = 1, num_lakec - c = filter_lakec(fc) + ! ------------------------------------------------------------------------ + ! Now add some other terms + ! ------------------------------------------------------------------------ + + ! Snow water content + do fc = 1, num_lakec + c = filter_lakec(fc) + + ice_mass(c) = ice_mass(c) + h2osno_no_layers(c) + do j = snl(c)+1,0 liquid_mass(c) = liquid_mass(c) + h2osoi_liq(c,j) ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) end do end do - if (add_lake_water_and_subtract_dynbal_baselines) then - ! Subtract baselines set in initialization + ! Soil water content of the soil under the lake + do j = 1, nlevgrnd 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) + liquid_mass(c) = liquid_mass(c) + h2osoi_liq(c,j) + ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j) end do - end if + end do end associate @@ -949,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) @@ -989,41 +1035,16 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, & 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 - ! 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) - - ! 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) = heat(c) - dynbal_baseline_heat(c) - end do - end associate end subroutine ComputeHeatLake !----------------------------------------------------------------------- subroutine AccumulateHeatLake(bounds, num_c, filter_c, & - temperature_inst, lakestate_inst, & + temperature_inst, lakestate_inst, & heat) ! ! !DESCRIPTION: @@ -1110,11 +1131,11 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, & end do end do - ! add ice heat and liquid heat together + ! 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) + heat(c) = heat(c) + (lake_heat_ice(c) + & + lake_heat_liquid(c) + lake_latent_heat_liquid(c)) end do end associate From 3ab74de5a347518f13129645b16ad35681218e60 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 25 Sep 2020 11:04:33 -0600 Subject: [PATCH 39/42] Add a single-point test that exercises dynamic lakes --- cime_config/testdefs/testlist_clm.xml | 10 ++++++++++ .../clm/smallville_dynlakes_monthly/include_user_mods | 1 + .../clm/smallville_dynlakes_monthly/user_nl_clm | 9 +++++++++ 3 files changed, 20 insertions(+) create mode 100644 cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/include_user_mods create mode 100644 cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 1d6552493f..520fb03fa1 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1392,6 +1392,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..49f6dc9e8b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm @@ -0,0 +1,9 @@ +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)' landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_c160127.nc landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_dynLakes_c200923.nc +! Key points are that lake area starts as 0, increases after the first +! year, then decreases after the second year. +! 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_c200923.nc' From 8ff6ca324634026f963506c6353568bf42f0c3a3 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 25 Sep 2020 14:49:11 -0600 Subject: [PATCH 40/42] Turn off methane in dynlakes test I was getting a gridcell-level methane balance check error when running ERS_Lm25.1x1_smallvilleIA.IHistClm50BgcCropQianRsGs.cheyenne_gnu.clm-smallville_dynlakes_monthly. I'm thinking it was due to ESCOMP/CTSM#43. So for now, I'm disabling methane in this testmod. (I thought about changing it to an Sp case, but (1) currently you're not allowed to run smallville in an Sp case, and (2) it would be nice to exercise as much of the BGC code as possible in this test, to exercise the dynamic landunit conservation code for BGC variables.) --- .../testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm | 4 ++++ src/biogeochem/ch4Mod.F90 | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) 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 index 49f6dc9e8b..7102f7052d 100644 --- 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 @@ -7,3 +7,7 @@ do_transient_lakes = .true. ! 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_c200923.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 From f2e2011652e9a99955c0f4b4353f61039d666f7c Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 28 Sep 2020 15:21:57 -0600 Subject: [PATCH 41/42] Point to new version of landuse timeseries file for dyn lakes test The important part of this is that it has PCT_LAKE + PCT_CROP <= 100% for all times. --- .../clm/smallville_dynlakes_monthly/user_nl_clm | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) 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 index 7102f7052d..b85d0ba6c1 100644 --- 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 @@ -1,12 +1,11 @@ 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)' landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_c160127.nc landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_dynLakes_c200923.nc -! Key points are that lake area starts as 0, increases after the first -! year, then decreases after the second year. -! 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_c200923.nc' +! 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. From d9b49729b7349f36d704ece617a13df4b99c0e63 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 28 Sep 2020 15:59:29 -0600 Subject: [PATCH 42/42] Fix unit tests --- src/biogeophys/CMakeLists.txt | 1 + .../test_dyn_cons_biogeophys.pf | 17 ++++++++++++----- 2 files changed, 13 insertions(+), 5 deletions(-) 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/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