diff --git a/cime_config/testdefs/testmods_dirs/clm/Hillslope_DominantPftLowland/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/Hillslope_DominantPftLowland/include_user_mods
new file mode 100644
index 0000000000..fa2e50a80d
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/Hillslope_DominantPftLowland/include_user_mods
@@ -0,0 +1 @@
+../Hillslope
diff --git a/cime_config/testdefs/testmods_dirs/clm/Hillslope_DominantPftLowland/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/Hillslope_DominantPftLowland/user_nl_clm
new file mode 100644
index 0000000000..4ab5c6cd56
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/Hillslope_DominantPftLowland/user_nl_clm
@@ -0,0 +1 @@
+hillslope_pft_distribution_method = 'DominantPftLowland'
diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90
index 42534d418c..03a0082e97 100644
--- a/src/main/surfrdMod.F90
+++ b/src/main/surfrdMod.F90
@@ -896,11 +896,12 @@ subroutine surfrd_hillslope(begg, endg, ncid, ns)
     ! !USES:
     use clm_instur, only : ncolumns_hillslope, wt_nat_patch
     use clm_varctl, only : nhillslope,max_columns_hillslope
-    use clm_varpar, only : natpft_size, natpft_lb
+    use clm_varpar, only : natpft_size, natpft_lb, natpft_ub
     use ncdio_pio,  only : ncd_inqdid, ncd_inqdlen
     use pftconMod , only : noveg
     use HillslopeHydrologyMod, only : pft_distribution_method, pft_standard, pft_from_file, pft_uniform_dominant_pft, pft_lowland_dominant_pft, pft_lowland_upland
     use array_utils, only: find_k_max_indices
+    use surfrdUtilsMod, only: collapse_to_dominant
 
     !
     ! !ARGUMENTS:
@@ -916,6 +917,8 @@ subroutine surfrd_hillslope(begg, endg, ncid, ns)
     logical  :: readvar                        ! is variable on dataset
     integer,pointer :: arrayl(:)               ! local array (needed because ncd_io expects a pointer)
     character(len=32) :: subname = 'surfrd_hillslope'  ! subroutine name
+    logical, allocatable :: do_not_collapse(:)
+    integer :: n_dominant
     !-----------------------------------------------------------------------
 
     ! number of hillslopes per landunit
@@ -960,52 +963,33 @@ subroutine surfrd_hillslope(begg, endg, ncid, ns)
           endif
        enddo
 
-    ! pft_uniform_dominant_pft uses the patch with the
-    ! largest weight for all hillslope columns in the gridcell
-    else if (pft_distribution_method == pft_uniform_dominant_pft) then
-       allocate(max_indices(1))
-       do g = begg, endg
-          ! If hillslopes will be used in a gridcell, modify wt_nat_patch,
-          ! otherwise use original patch distribution
-          if(ncolumns_hillslope(g) > 0) then
-
-             call find_k_max_indices(wt_nat_patch(g,:),natpft_lb,1,max_indices)
-             wt_nat_patch(g,:) = 0._r8
-             wt_nat_patch(g,max_indices(1)) = 100._r8
+    else if (pft_distribution_method == pft_uniform_dominant_pft &
+        .or. pft_distribution_method == pft_lowland_dominant_pft) then
 
-          endif
-       enddo
-       deallocate(max_indices)
-
-    ! pft_lowland_dominant_pft uses the two patches with the
-    ! largest weights for the hillslope columns in the gridcell
-    else if (pft_distribution_method == pft_lowland_dominant_pft) then
-       allocate(max_indices(2))
+       ! If hillslopes will be used in a gridcell, modify wt_nat_patch,
+       ! otherwise use original patch distribution
+       allocate(do_not_collapse(begg:endg))
+       do_not_collapse(begg:endg) = .false.
        do g = begg, endg
-          ! If hillslopes will be used in a gridcell, modify wt_nat_patch, otherwise use original patch distribution
-          if(ncolumns_hillslope(g) > 0) then
+          if (ncolumns_hillslope(g) == 0) then
+             do_not_collapse(g) = .true.
+          end if
+       end do
 
-             ! Preserve the relative weights of the largest and
-             ! next largest weights using arbitrarily chosen values
-             ! (i.e. 1 should be larger than 2) that sum to 100.
-             ! This will minimize memory usage while still allowing
-             ! HillslopeDominantLowlandPft to pick out the two largest patch types.
-
-             call find_k_max_indices(wt_nat_patch(g,:),natpft_lb,2,max_indices)
-             ! check that 2nd index weight is non-zero
-             if (wt_nat_patch(g,max_indices(2)) > 0._r8) then
-                wt_nat_patch(g,:) = 0._r8
-                wt_nat_patch(g,max_indices(1)) = 75._r8
-                wt_nat_patch(g,max_indices(2)) = 25._r8
-             else
-                ! if only one pft exists, set its weight to 100 per cent
-                wt_nat_patch(g,:) = 0._r8
-                wt_nat_patch(g,max_indices(1)) = 100._r8
-             endif
+       if (pft_distribution_method == pft_uniform_dominant_pft) then
+         ! pft_uniform_dominant_pft uses the patch with the
+         ! largest weight for all hillslope columns in the gridcell
+         n_dominant = 1
+       else if (pft_distribution_method == pft_lowland_dominant_pft) then
+         ! pft_lowland_dominant_pft uses the two patches with the
+         ! largest weights for the hillslope columns in the gridcell
+         n_dominant = 2
+       else
+          call endrun( msg=' ERROR: unrecognized hillslope_pft_distribution_method'//errMsg(sourcefile, __LINE__))
+       end if
 
-          endif
-       enddo
-       deallocate(max_indices)
+       call collapse_to_dominant(wt_nat_patch(begg:endg,:), natpft_lb, natpft_ub, begg, endg, n_dominant, do_not_collapse)
+       deallocate(do_not_collapse)
 
     else if (pft_distribution_method /= pft_standard) then
       call endrun( msg=' ERROR: unrecognized hillslope_pft_distribution_method'//errMsg(sourcefile, __LINE__))