Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add subgrid hillslope hydrology (NGEE Arctic IM2) #6718

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions cime_config/tests.py
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@
"SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-per_crop",
"SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-fan",
"SMS.r05_r05.IELM.elm-topounit",
"SMS.r05_r05.IELM.elm-topounit_im2",
"ERS.ELM_USRDAT.I1850ELM.elm-usrdat",
"ERS.r05_r05.IELM.elm-lnd_rof_2way",
"ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_default_I1850CNPRDCTCBC",
Expand Down
9 changes: 9 additions & 0 deletions components/elm/bld/namelist_files/namelist_definition.xml
Original file line number Diff line number Diff line change
Expand Up @@ -1953,6 +1953,15 @@ If TRUE setup up lateral grid connectivity in CLM.
Specifies the method for decomposing CLM grids across processors.
</entry>

<!-- ======================================================================================== -->
<!-- Namelist options related to within-gridcell hillslope hydrologic connectivity -->
<!-- ======================================================================================== -->
<entry id="use_IM2_hillslope_hydrology" type="logical" category="default_settings"
group="elm_inparm" valid_values=".true.,.false" value=".false." >
If TRUE use within-gridcell hillslope hydrologic connectivity based on drainage terms at the
topounit level. Based on the IM2 implementation from NGEE Arctic project.
</entry>

<!-- ======================================================================================== -->
<!-- Namelist options related to the bgc & pflotran interface -->
<!-- ======================================================================================== -->
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
#!/bin/bash
./xmlchange --id ELM_BLDNML_OPTS --val "-bgc sp -topounit"
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata/half_degree_merge_surfdata_0.5x0.5_simyr2000_c190418.with_aveDTB.20201222.nc'
use_IM2_hillslope_hydrology = .true.
29 changes: 23 additions & 6 deletions components/elm/src/biogeophys/BalanceCheckMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ module BalanceCheckMod
use GridcellType , only : grc_pp
use TopounitType , only : top_pp
use GridcellDataType , only : grc_ef, grc_wf, grc_ws
use TopounitDataType , only : top_af ! atmospheric flux variables
use TopounitDataType , only : top_af, top_ws
use LandunitType , only : lun_pp
use ColumnType , only : col_pp
use ColumnDataType , only : col_ef, col_ws, col_wf
Expand Down Expand Up @@ -168,7 +168,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall
use column_varcon , only : icol_road_perv, icol_road_imperv
use landunit_varcon , only : istice_mec, istdlak, istsoil,istcrop,istwet
use elm_varctl , only : create_glacier_mec_landunit
use elm_varctl , only : create_glacier_mec_landunit, use_IM2_hillslope_hydrology
use elm_initializeMod , only : surfalb_vars
use CanopyStateType , only : canopystate_type
use subgridAveMod
Expand Down Expand Up @@ -232,6 +232,8 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
qflx_h2osfc_to_ice => col_wf%qflx_h2osfc_to_ice , & ! Input: [real(r8) (:) ] conversion of h2osfc to ice
qflx_drain_perched => col_wf%qflx_drain_perched , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s)
qflx_floodc => col_wf%qflx_floodc , & ! Input: [real(r8) (:) ] total runoff due to flooding
qflx_from_uphill => col_wf%qflx_from_uphill , & ! Input: [real(r8) (:) ] received from uphill topounit
qflx_to_downhill => col_wf%qflx_to_downhill , & ! Input: [real(r8) (:) ] sent to downhill topounit
qflx_h2osfc_surf => col_wf%qflx_h2osfc_surf , & ! Input: [real(r8) (:) ] surface water runoff (mm/s)
qflx_snow_melt => col_wf%qflx_snow_melt , & ! Input: [real(r8) (:) ] snow melt (net)
qflx_surf => col_wf%qflx_surf , & ! Input: [real(r8) (:) ] surface runoff (mm H2O /s)
Expand Down Expand Up @@ -308,17 +310,30 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
forc_rain_col(c) = forc_rain(t)
forc_snow_col(c) = forc_snow(t)
end if

! update the topounit-level from_uphill water state
! when using fraction_from_uphill < 1.0, this state can get very small during recession
! which can lead to negative state from roundoff error. Trap this with a max(), and force to zero
! when the state gets too small.
if (use_IM2_hillslope_hydrology) then
top_ws%from_uphill(t) = max(0._r8, top_ws%from_uphill(t) - (col_wf%qflx_from_uphill(c) * dtime))
if (top_ws%from_uphill(t) < 1.e-20_r8) then
top_ws%from_uphill(t) = 0._r8
endif
endif
end do

! Water balance check

do c = bounds%begc, bounds%endc

! add qflx_drain_perched and qflx_flood
! add qflx_from_uphill and qflx_to_downhill
if (col_pp%active(c)) then
errh2o(c) = endwb(c) - begwb(c) &
- (forc_rain_col(c) + forc_snow_col(c) + qflx_floodc(c) + qflx_surf_irrig_col(c) + qflx_over_supply_col(c) &
- qflx_evap_tot(c) - qflx_surf(c) - qflx_h2osfc_surf(c) &
- (forc_rain_col(c) + forc_snow_col(c) + qflx_floodc(c) + qflx_from_uphill(c) &
+ qflx_surf_irrig_col(c) + qflx_over_supply_col(c) &
- qflx_evap_tot(c) - qflx_surf(c) - qflx_h2osfc_surf(c) - qflx_to_downhill(c) &
- qflx_qrgwl(c) - qflx_drain(c) - qflx_drain_perched(c) - qflx_snwcp_ice(c) &
- qflx_lateral(c) + qflx_h2orof_drain(c)) * dtime
dwb(c) = (endwb(c)-begwb(c))/dtime
Expand Down Expand Up @@ -960,6 +975,8 @@ subroutine GridBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
qflx_h2osfc_to_ice => col_wf%qflx_h2osfc_to_ice , & ! Input: [real(r8) (:) ] conversion of h2osfc to ice
qflx_drain_perched => col_wf%qflx_drain_perched , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s)
qflx_floodc => col_wf%qflx_floodc , & ! Input: [real(r8) (:) ] total runoff due to flooding
qflx_from_uphill => col_wf%qflx_from_uphill , & ! Input: [real(r8) (:) ] received from uphill topounit (mm/s)
qflx_to_downhill => col_wf%qflx_to_downhill , & ! Input: [real(r8) (:) ] sent to downhill topounit (mm/s)
qflx_h2osfc_surf => col_wf%qflx_h2osfc_surf , & ! Input: [real(r8) (:) ] surface water runoff (mm/s)
qflx_snow_melt => col_wf%qflx_snow_melt , & ! Input: [real(r8) (:) ] snow melt (net)
qflx_surf => col_wf%qflx_surf , & ! Input: [real(r8) (:) ] surface runoff (mm H2O /s)
Expand Down Expand Up @@ -1053,8 +1070,8 @@ subroutine GridBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
if (col_pp%active(c)) then

qflx_net_col(c) = &
- forc_rain_col(c) - forc_snow_col(c) - qflx_floodc(c) - qflx_irrig(c) &
+ qflx_evap_tot(c) + qflx_surf(c) + qflx_h2osfc_surf(c) &
- forc_rain_col(c) - forc_snow_col(c) - qflx_floodc(c) - qflx_from_uphill(c) - qflx_irrig(c) &
+ qflx_evap_tot(c) + qflx_surf(c) + qflx_h2osfc_surf(c) + qflx_to_downhill(c) &
+ qflx_qrgwl(c) + qflx_drain(c) + qflx_drain_perched(c) + qflx_snwcp_ice(c) &
+ qflx_lateral(c)

Expand Down
4 changes: 2 additions & 2 deletions components/elm/src/biogeophys/CanopyHydrologyMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -282,8 +282,8 @@ subroutine CanopyHydrology(bounds, &
ltype(l)==istcrop) then

qflx_candrip(p) = 0._r8 ! rate of canopy runoff
qflx_through_snow(p) = 0._r8 ! rain precipitation direct through canopy
qflx_through_rain(p) = 0._r8 ! snow precipitation direct through canopy
qflx_through_snow(p) = 0._r8 ! snow precipitation direct through canopy
qflx_through_rain(p) = 0._r8 ! rain precipitation direct through canopy
qflx_prec_intr(p) = 0._r8 ! total intercepted precipitation
fracsnow(p) = 0._r8 ! fraction of input precip that is snow
fracrain(p) = 0._r8 ! fraction of input precip that is rain
Expand Down
43 changes: 39 additions & 4 deletions components/elm/src/biogeophys/HydrologyDrainageMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -49,15 +49,16 @@ subroutine HydrologyDrainage(bounds, &
!$acc routine seq
use landunit_varcon , only : istice, istwet, istsoil, istice_mec, istcrop
use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv, icol_sunwall, icol_shadewall
use elm_varcon , only : denh2o, denice, secspday
use elm_varcon , only : denh2o, denice, secspday, frac_to_downhill
use elm_varctl , only : glc_snow_persistence_max_days, use_vichydro, use_betr
!use domainMod , only : ldomain
use elm_varsur , only : f_surf
use TopounitType , only : top_pp
use TopounitDataType , only : top_ws
use atm2lndType , only : atm2lnd_type
use elm_varpar , only : nlevgrnd, nlevurb, nlevsoi
use SoilHydrologyMod , only : ELMVICMap, Drainage
use elm_varctl , only : use_vsfm
use elm_varctl , only : use_vsfm, use_IM2_hillslope_hydrology
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
Expand All @@ -77,7 +78,8 @@ subroutine HydrologyDrainage(bounds, &
!
! !LOCAL VARIABLES:
real(r8) :: dtime
integer :: g,t,l,c,j,fc,tpu_ind ! indices
real(r8) :: temp_to_downhill, temp_mass
integer :: g,t,l,c,j,fc,tpu_ind, downhill_t ! indices
!-----------------------------------------------------------------------

associate( &
Expand Down Expand Up @@ -120,7 +122,8 @@ subroutine HydrologyDrainage(bounds, &
qflx_runoff_r => col_wf%qflx_runoff_r , & ! Output: [real(r8) (:) ] Rural total runoff (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s)
qflx_snwcp_ice => col_wf%qflx_snwcp_ice , & ! Output: [real(r8) (:) ] excess snowfall due to snow capping (mm H2O /s) [+]`
qflx_glcice => col_wf%qflx_glcice , & ! Output: [real(r8) (:) ] flux of new glacier ice (mm H2O /s)
qflx_glcice_frz => col_wf%qflx_glcice_frz & ! Output: [real(r8) (:) ] ice growth (positive definite) (mm H2O/s)
qflx_glcice_frz => col_wf%qflx_glcice_frz , & ! Output: [real(r8) (:) ] ice growth (positive definite) (mm H2O/s)
qflx_to_downhill => col_wf%qflx_to_downhill & ! Output: [real(r8) (:) ] flux transferred to downhill topounit (mm H2O/s)
)

! Determine time step and step size
Expand Down Expand Up @@ -281,6 +284,38 @@ subroutine HydrologyDrainage(bounds, &

end if

! if using topounit hillslope hydrology, fractions of qflx_surf, qflx_drain_perched, and qflx_h2osfc
! are passed to the from_uphill water state on the downhill topounit, via qflx_to_downhill
! only shift positive fluxes, and only set fluxes if there is a downhill topounit
if (use_IM2_hillslope_hydrology) then
downhill_t = top_pp%downhill_ti(t)
if (downhill_t /= -1) then
! shift a fixed fraction of qflx_surf
temp_to_downhill = max(0._r8, frac_to_downhill * qflx_surf(c))
qflx_to_downhill(c) = temp_to_downhill
qflx_surf(c) = qflx_surf(c) - temp_to_downhill

! shift a fixed fraction of qflx_drain_perched
temp_to_downhill = max(0._r8, frac_to_downhill * qflx_drain_perched(c))
qflx_to_downhill(c) = qflx_to_downhill(c) + temp_to_downhill
qflx_drain_perched(c) = qflx_drain_perched(c) - temp_to_downhill

! shift a fixed fraction of qflx_h2osfc_surf
temp_to_downhill = max(0._r8, frac_to_downhill * qflx_h2osfc_surf(c))
qflx_to_downhill(c) = qflx_to_downhill(c) + temp_to_downhill
qflx_h2osfc_surf(c) = qflx_h2osfc_surf(c) - temp_to_downhill

! update the downhill topounit water state
! relative weight of this column to downhill topounit on the gridcell is used to
! scale the mass of water from this column to a potentially larger or smaller
! downhill topounit
temp_mass = (qflx_to_downhill(c) * dtime) * (col_pp%wtgcell(c) / top_pp%wtgcell(downhill_t))
top_ws%from_uphill(downhill_t) = top_ws%from_uphill(downhill_t) + temp_mass
else
qflx_to_downhill(c) = 0._r8
endif
endif

qflx_runoff(c) = qflx_drain(c) + qflx_surf(c) + qflx_h2osfc_surf(c) + qflx_qrgwl(c) + qflx_drain_perched(c)

if ((lun_pp%itype(l)==istsoil .or. lun_pp%itype(l)==istcrop) .and. col_pp%active(c)) then
Expand Down
32 changes: 29 additions & 3 deletions components/elm/src/biogeophys/SoilHydrologyMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ module SoilHydrologyMod
use SoilHydrologyType , only : soilhydrology_type
use SoilStateType , only : soilstate_type
use WaterfluxType , only : waterflux_type
use TopounitType , only : top_pp
use TopounitDataType , only : top_ws
use LandunitType , only : lun_pp
use ColumnType , only : col_pp
use ColumnDataType , only : col_es, col_ws, col_wf
Expand Down Expand Up @@ -47,12 +49,12 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
!
! !USES:
!$acc routine seq
use elm_varcon , only : denice, denh2o, wimp, pondmx_urban
use elm_varcon , only : denice, denh2o, wimp, pondmx_urban, frac_from_uphill
use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall
use column_varcon , only : icol_road_imperv, icol_road_perv
use elm_varpar , only : nlevsoi, nlevgrnd, maxpatch_pft
use elm_varpar , only : nlayer, nlayert
use elm_varctl , only : use_var_soil_thick
use elm_varctl , only : use_var_soil_thick, use_IM2_hillslope_hydrology
use SoilWaterMovementMod, only : zengdecker_2009_with_var_soil_thick
!
! !ARGUMENTS:
Expand All @@ -66,7 +68,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
real(r8), intent(in) :: dtime
!
! !LOCAL VARIABLES:
integer :: c,j,fc,g,l,i !indices
integer :: c,j,fc,g,l,t,i !indices
integer :: nlevbed !# levels to bedrock
real(r8) :: xs(bounds%begc:bounds%endc) !excess soil water above urban ponding limit
real(r8) :: vol_ice(bounds%begc:bounds%endc,1:nlevgrnd) !partial volume of ice lens in layer
Expand Down Expand Up @@ -101,6 +103,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
qflx_floodc => col_wf%qflx_floodc , & ! Input: [real(r8) (:) ] column flux of flood water from RTM
qflx_evap_grnd => col_wf%qflx_evap_grnd , & ! Input: [real(r8) (:) ] ground surface evaporation rate (mm H2O/s) [+]
qflx_top_soil => col_wf%qflx_top_soil , & ! Output: [real(r8) (:) ] net water input into soil from top (mm/s)
qflx_from_uphill => col_wf%qflx_from_uphill , & ! Output: [real(r8) (:) ] water received from uphill topounit (mm/s)
qflx_surf => col_wf%qflx_surf , & ! Output: [real(r8) (:) ] surface runoff (mm H2O /s)
qflx_irrig => col_wf%qflx_irrig , & ! Input: [real(r8) (:) ] irrigation flux (mm H2O /s)
irrig_rate => veg_wf%irrig_rate , & ! Input: [real(r8) (:) ] current irrigation rate (applied if !n_irrig_steps_left > 0) [mm/s]
Expand Down Expand Up @@ -241,11 +244,34 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &

end do

! calculate the sum of column weights on each topounit for columns in the hydrologyc filter
! This will be the istsoil, istcrop, and icol_road_perv subset of urban columns
! First zero the topounit sum of weights. This zeros some multiple times, but no harm done.
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
t = col_pp%topounit(c)
top_pp%uphill_wt(t) = 0._r8
end do
! Next sum the weights
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
t = col_pp%topounit(c)
top_pp%uphill_wt(t) = top_pp%uphill_wt(t) + col_pp%wttopounit(c)
end do
Comment on lines +250 to +260
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Couldn't these weights be computed during initialization?


! remove stormflow and snow on h2osfc from qflx_top_soil
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
t = col_pp%topounit(c)
! add flood water flux to qflx_top_soil
qflx_top_soil(c) = qflx_top_soil(c) + qflx_snow_h2osfc(c) + qflx_floodc(c)

! flow from uphill topounit goes to top of soil for soil, crop, and pervious road columns
if (use_IM2_hillslope_hydrology) then
qflx_from_uphill(c) = (col_pp%wttopounit(c)/top_pp%uphill_wt(t)) * (frac_from_uphill * top_ws%from_uphill(t)) / dtime
qflx_top_soil(c) = qflx_top_soil(c) + qflx_from_uphill(c)
endif

end do
Comment on lines 263 to 275
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about avoiding the if (use_IM2_hillslope_hydrology) within the do-loop by leaving the original code as is and then adding a separate if-loop like

do fc = 1, num_hydrologyc
   ... original code
enddo

if (use_IM2_hillslope_hydrology) then
   .... new code
   do fc = 1, ...
   enddo
endif


end associate
Expand Down
14 changes: 14 additions & 0 deletions components/elm/src/data_types/ColumnDataType.F90
Original file line number Diff line number Diff line change
Expand Up @@ -515,6 +515,8 @@ module ColumnDataType
real(r8), pointer :: qflx_irr_demand (:) => null() ! col surface irrigation demand (mm H2O /s)
real(r8), pointer :: qflx_over_supply (:) => null() ! col over supplied irrigation
real(r8), pointer :: qflx_h2orof_drain (:) => null() ! drainage from floodplain inundation volume (mm H2O/s))
real(r8), pointer :: qflx_from_uphill (:) => null() ! input to top soil layer from uphill topounit(s) (mm H2O/s))
real(r8), pointer :: qflx_to_downhill (:) => null() ! output from column to the downhill topounit (mm H2O/s))

real(r8), pointer :: mflx_infl_1d (:) => null() ! infiltration source in top soil control volume (kg H2O /s)
real(r8), pointer :: mflx_dew_1d (:) => null() ! liquid+snow dew source in top soil control volume (kg H2O /s)
Expand Down Expand Up @@ -5737,6 +5739,8 @@ subroutine col_wf_init(this, begc, endc)
allocate(this%qflx_over_supply (begc:endc)) ; this%qflx_over_supply (:) = spval
allocate(this%qflx_irr_demand (begc:endc)) ; this%qflx_irr_demand (:) = spval
allocate(this%qflx_h2orof_drain (begc:endc)) ; this%qflx_h2orof_drain (:) = spval
allocate(this%qflx_from_uphill (begc:endc)) ; this%qflx_from_uphill (:) = spval
allocate(this%qflx_to_downhill (begc:endc)) ; this%qflx_to_downhill (:) = spval

!VSFM variables
ncells = endc - begc + 1
Expand Down Expand Up @@ -5843,6 +5847,14 @@ subroutine col_wf_init(this, begc, endc)
avgflag='A', long_name='column-integrated snow freezing rate', &
ptr_col=this%qflx_snofrz, set_lake=spval, c2l_scale_type='urbanf', default='inactive')

call hist_addfld1d (fname='QFROM_UPHILL', units='mm/s', &
avgflag='A', long_name='input to top layer soil from uphill topounit(s)', &
ptr_col=this%qflx_from_uphill, c2l_scale_type='urbanf')

call hist_addfld1d (fname='QTO_DOWNHILL', units='mm/s', &
avgflag='A', long_name='output from column to downhill topounit', &
ptr_col=this%qflx_to_downhill, c2l_scale_type='urbanf')

Comment on lines +5850 to +5857
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If these fields are added by default, let's only add them if the new physics is turned on.

if (create_glacier_mec_landunit) then
this%qflx_glcice(begc:endc) = spval
call hist_addfld1d (fname='QICE', units='mm/s', &
Expand Down Expand Up @@ -5890,6 +5902,8 @@ subroutine col_wf_init(this, begc, endc)
this%qflx_grnd_irrig(begc:endc) = 0._r8
this%qflx_over_supply(begc:endc) = 0._r8
this%qflx_h2orof_drain(begc:endc)= 0._r8
this%qflx_from_uphill(begc:endc) = 0._r8
this%qflx_to_downhill(begc:endc) = 0._r8
! needed for CNNLeaching
do c = begc, endc
l = col_pp%landunit(c)
Expand Down
Loading