Merge branch 'thorntonpe/lnd/IM2' (PR #6718)

This new capability builds on the existing topographic unit (topounit) subgrid scheme by
implementing transport of water within a gridcell from higher to lower elevation topounits.
Assuming that the topounit information on an existing surface dataset includes mean
elevation for each topounit on a gridcell, the code will calculate topounit connectivity (topology)
 at run time during model initialization. On each time step a user-specified fraction of the
surface runoff, ponded water runoff, and perched water table drainage terms are shifted
to the next-lower topounit, where that water is distributed to the columns on the topounit
with a user-defined decay term.

[BFB]
[FCC]
[NML]

Conflicts:
	components/elm/src/biogeophys/HydrologyDrainageMod.F90
	components/elm/src/data_types/ColumnDataType.F90

cime_config/tests.py

@@ -95,6 +95,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.r05_r05.IELM.elm-V2_ELM_MOSART_features",

components/elm/bld/namelist_files/namelist_definition.xml

@@ -1959,6 +1959,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                                  -->
 <!-- ========================================================================================  -->

components/elm/cime_config/testdefs/testmods_dirs/elm/topounit_im2/shell_commands

@@ -0,0 +1,2 @@
+#!/bin/bash
+./xmlchange --id ELM_BLDNML_OPTS --val "-bgc sp -topounit"

components/elm/cime_config/testdefs/testmods_dirs/elm/topounit_im2/user_nl_elm

@@ -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.

components/elm/src/biogeophys/BalanceCheckMod.F90

@@ -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
@@ -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
@@ -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)
@@ -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
@@ -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)
@@ -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)
 

components/elm/src/biogeophys/CanopyHydrologyMod.F90

@@ -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

components/elm/src/biogeophys/HydrologyDrainageMod.F90

@@ -49,15 +49,16 @@ subroutine HydrologyDrainage(bounds,    &
       !$acc routine seq
     use landunit_varcon  , only : istice, istwet, istsoil, istice_mec, istcrop, istice
     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
@@ -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(                                                                  &
@@ -122,7 +124,8 @@ subroutine HydrologyDrainage(bounds,    &
          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_diag       => col_wf%qflx_glcice_diag        , & ! Output: [real(r8) (:)   ]  flux of new glacier ice (mm H2O/s) - diagnostic, no MECs or GLC
-         qflx_glcice_frz_diag   => col_wf%qflx_glcice_frz_diag      & ! Output: [real(r8) (:)   ]  ice growth (positive definite) (mm H2O/s)) - diagnostic, no MECs or GLC
+         qflx_glcice_frz_diag   => col_wf%qflx_glcice_frz_diag    , & ! Output: [real(r8) (:)   ]  ice growth (positive definite) (mm H2O/s)) - diagnostic, no MECs or GLC
+         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
@@ -289,6 +292,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

components/elm/src/biogeophys/SoilHydrologyMod.F90

@@ -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
@@ -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:
@@ -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
@@ -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]
@@ -244,10 +247,39 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
       ! 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)
       end do
 
+      ! when using the subgrid hillslope lateral flow mechanism (IM2 from NGEE Arctic):
+      ! 1. calculate fraction of topounit flow that goes to each column
+      ! 2. Calculate the weighted flux from uphill (qflx_from_uphill) and add to qflx_top_soil 
+      if (use_IM2_hillslope_hydrology) then
+         ! 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 relevant column 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
+
+         ! flow from uphill topounit goes to top of soil for soil, crop, and pervious road columns
+         do fc = 1, num_hydrologyc
+            c = filter_hydrologyc(fc)
+            t = col_pp%topounit(c)
+            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)
+         end do
+      endif
+
     end associate
 
    end subroutine SurfaceRunoff