Skip to content

Commit

Permalink
Add lake water and heat content to gridcell total water and heat comp…
Browse files Browse the repository at this point in the history 
…utation
  • Loading branch information
@Ivanderkelen @billsacks
Ivanderkelen authored and billsacks committed Jun 19, 2020
1 parent dc6c885 commit 2ea9eee
Show file tree
Hide file tree
Showing 6 changed files with 87 additions and 89 deletions.
6 changes: 4 additions & 2 deletions src/biogeophys/BalanceCheckMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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
!
Expand Down Expand Up @@ -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))

Expand Down
4 changes: 2 additions & 2 deletions src/biogeophys/LakeHydrologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))

Expand Down
130 changes: 57 additions & 73 deletions src/biogeophys/TotalWaterAndHeatMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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 = &
Expand Down Expand Up @@ -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)
!
Expand All @@ -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.
Expand All @@ -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))
Expand Down Expand Up @@ -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)
!
Expand All @@ -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.
Expand All @@ -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'
!-----------------------------------------------------------------------

Expand All @@ -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)
)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -816,22 +836,19 @@ 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
integer , intent(in) :: 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:
Expand All @@ -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'
!-----------------------------------------------------------------------

Expand All @@ -853,15 +871,16 @@ 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)
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)
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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 2ea9eee

Please sign in to comment.