Functionize & unit-test HillslopeSoilThicknessProfile_linear().

Squash-"merge" of hillslope_hydrology-functionize_linear into hillslope_hydrology-ssr3

src/biogeophys/CMakeLists.txt

@@ -8,6 +8,7 @@ list(APPEND clm_sources
   CanopyStateType.F90
   EnergyFluxType.F90
   GlacierSurfaceMassBalanceMod.F90
+  HillslopeHydrologyUtilsMod.F90
   HumanIndexMod.F90
   InfiltrationExcessRunoffMod.F90
   IrrigationMod.F90

src/biogeophys/HillslopeHydrologyMod.F90

@@ -14,6 +14,7 @@ module HillslopeHydrologyMod
   use clm_varctl     , only : use_hillslope_routing
   use decompMod      , only : bounds_type
   use clm_varcon     , only : rpi
+  use HillslopeHydrologyUtilsMod, only : HillslopeSoilThicknessProfile_linear
 
   ! !PUBLIC TYPES:
   implicit none
@@ -703,7 +704,6 @@ subroutine HillslopeSoilThicknessProfile(bounds,&
     real(r8)              :: soil_depth_upland
     real(r8), parameter   :: soil_depth_lowland_default = 8.0
     real(r8), parameter   :: soil_depth_upland_default  = 8.0
-    real(r8), parameter   :: toosmall_distance  = 1e-6
 
     character(len=*), parameter :: subname = 'HillslopeSoilThicknessProfile'
 
@@ -746,30 +746,7 @@ subroutine HillslopeSoilThicknessProfile(bounds,&
        end do
     ! Linear soil thickness profile
     else if (soil_profile_method == soil_profile_linear) then
-       do l = bounds%begl,bounds%endl
-          min_hill_dist = minval(col%hill_distance(lun%coli(l):lun%colf(l)))
-          max_hill_dist = maxval(col%hill_distance(lun%coli(l):lun%colf(l)))
-
-          if (abs(max_hill_dist - min_hill_dist) > toosmall_distance) then
-             m = (soil_depth_lowland - soil_depth_upland)/ &
-                  (max_hill_dist - min_hill_dist)
-          else
-             m = 0._r8
-          end if
-          b = soil_depth_upland
-
-          do c =  lun%coli(l), lun%colf(l)
-             if (col%is_hillslope_column(c) .and. col%active(c)) then
-                soil_depth_col = m*(max_hill_dist - col%hill_distance(c)) + b
-
-                do j = 1,nlevsoi
-                   if ((zisoi(j-1) <  soil_depth_col) .and. (zisoi(j) >= soil_depth_col)) then
-                      col%nbedrock(c) = j
-                   end if
-                enddo
-             end if
-          enddo
-       enddo
+       call HillslopeSoilThicknessProfile_linear(bounds, soil_depth_lowland, soil_depth_upland)
     else if (masterproc) then
        call endrun( 'ERROR:: invalid soil_profile_method.'//errmsg(sourcefile, __LINE__) )
     end if

src/biogeophys/HillslopeHydrologyUtilsMod.F90

@@ -0,0 +1,76 @@
+module HillslopeHydrologyUtilsMod
+
+  !-----------------------------------------------------------------------
+  ! !DESCRIPTION:
+  ! Utilities used in HillslopeHydrologyMod
+  !
+  ! !USES:
+#include "shr_assert.h"
+  use decompMod      , only : bounds_type
+  use shr_kind_mod   , only : r8 => shr_kind_r8
+  use shr_log_mod    , only : errMsg => shr_log_errMsg
+  use spmdMod        , only : masterproc, iam
+  use abortutils     , only : endrun
+  use clm_varctl     , only : iulog
+
+  ! !PUBLIC TYPES:
+  implicit none
+
+  private
+  save
+
+  ! !PUBLIC MEMBER FUNCTIONS:
+  public HillslopeSoilThicknessProfile_linear
+
+contains
+
+  !------------------------------------------------------------------------
+  subroutine HillslopeSoilThicknessProfile_linear(bounds, soil_depth_lowland, soil_depth_upland)
+    !
+    ! !DESCRIPTION:
+    ! Modify soil thickness across hillslope by changing
+    ! nbedrock according to the "Linear" method
+    !
+    ! !USES:
+    use LandunitType    , only : lun
+    use ColumnType      , only : col
+    use clm_varpar      , only : nlevsoi
+    use clm_varcon      , only : zisoi
+    !
+    ! !ARGUMENTS:
+    type(bounds_type), intent(in) :: bounds
+    real(r8), intent(in) :: soil_depth_lowland, soil_depth_upland
+    !
+    ! !LOCAL VARIABLES
+    real(r8) :: min_hill_dist, max_hill_dist
+    real(r8) :: soil_depth_col
+    real(r8) :: m, b
+    integer :: c, j, l
+    real(r8), parameter :: toosmall_distance  = 1e-6
+
+    do l = bounds%begl,bounds%endl
+       min_hill_dist = minval(col%hill_distance(lun%coli(l):lun%colf(l)))
+       max_hill_dist = maxval(col%hill_distance(lun%coli(l):lun%colf(l)))
+
+       if (abs(max_hill_dist - min_hill_dist) > toosmall_distance) then
+          m = (soil_depth_lowland - soil_depth_upland)/ &
+               (max_hill_dist - min_hill_dist)
+       else
+          m = 0._r8
+       end if
+       b = soil_depth_upland
+
+       do c =  lun%coli(l), lun%colf(l)
+          write(iulog, *) 'c = ',c
+          if (col%is_hillslope_column(c) .and. col%active(c)) then
+             soil_depth_col = m*(max_hill_dist - col%hill_distance(c)) + b
+             do j = 1,nlevsoi
+               if ((zisoi(j-1) <  soil_depth_col) .and. (zisoi(j) >= soil_depth_col)) then
+                  col%nbedrock(c) = j
+               end if
+             enddo
+          end if
+       enddo
+    enddo
+   end subroutine HillslopeSoilThicknessProfile_linear
+end module HillslopeHydrologyUtilsMod

src/biogeophys/test/CMakeLists.txt

@@ -1,6 +1,7 @@
 add_subdirectory(Daylength_test)
 add_subdirectory(Irrigation_test)
 add_subdirectory(HumanStress_test)
+add_subdirectory(HillslopeHydrology_test)
 add_subdirectory(SnowHydrology_test)
 add_subdirectory(Photosynthesis_test)
 add_subdirectory(Balance_test)

src/biogeophys/test/HillslopeHydrology_test/CMakeLists.txt

@@ -0,0 +1,6 @@
+set (pfunit_sources
+  test_hillslopehydrologyUtils.pf)
+
+add_pfunit_ctest(HillslopeHydrologyUtils
+  TEST_SOURCES "${pfunit_sources}"
+  LINK_LIBRARIES clm csm_share esmf_wrf_timemgr)

src/biogeophys/test/HillslopeHydrology_test/test_hillslopehydrologyUtils.pf

@@ -0,0 +1,168 @@
+module test_hillslopehydrologyUtils
+
+  ! Tests of the HillslopeHydrologyUtils module
+
+  use funit
+  use unittestSubgridMod
+  use ColumnType      , only : col
+  use LandunitType    , only : lun
+  use landunit_varcon , only : istwet
+  use decompMod       , only : bounds_type
+  use clm_varpar      , only : nlevgrnd
+  use shr_kind_mod    , only : r8 => shr_kind_r8
+  use HillslopeHydrologyUtilsMod, only : HillslopeSoilThicknessProfile_linear
+
+  implicit none
+
+  ! From clm_instInit
+  real(r8), parameter :: soil_depth_lowland = 8.5_r8
+  real(r8), parameter :: soil_depth_upland = 2._r8
+
+  integer, parameter :: nbedrock_dummy_value = 9999
+
+  @TestCase
+  type, extends(TestCase) :: TestInit
+   contains
+     procedure :: setUp
+     procedure :: tearDown
+  end type TestInit
+
+contains
+
+  subroutine setUp(this)
+    ! Set up variables needed for tests: various subgrid type variables, along with
+    ! bounds.
+    !
+    use clm_varcon, only: clm_varcon_init, zisoi
+    use clm_varpar, only: nlevsoi, nlevgrnd
+    class(TestInit), intent(inout) :: this
+    integer :: g, l, c
+    real(r8), allocatable :: my_zisoi(:)
+
+    ! Set up subgrid structure
+    ! The weights (of both landunits and columns) and column types in the following are
+    ! arbitrary, since they are not important for these tests
+
+    call unittest_subgrid_setup_start()
+
+    ! Set up gridcell with one landunit and two columns
+    call unittest_add_gridcell()
+    call unittest_add_landunit(my_gi=gi, ltype=istwet, wtgcell=0.25_r8)
+    call unittest_add_column(my_li=li, ctype=1, wtlunit=0.5_r8)
+    call unittest_add_column(my_li=li, ctype=1, wtlunit=0.5_r8)
+
+    call unittest_subgrid_setup_end()
+
+    ! These will be enabled by specific tests
+    col%active(begc:endc) = .false.
+    col%is_hillslope_column(begc:endc) = .false.
+
+    ! Set up ground/soil structure
+    nlevsoi = 5
+    allocate(my_zisoi(1:nlevsoi))
+    my_zisoi = [0.01_r8, 0.02_r8, 2._r8, 4._r8, 6._r8]
+    nlevgrnd = size(my_zisoi)
+    call clm_varcon_init( is_simple_buildtemp = .true.)
+    zisoi(0) = 0._r8
+    zisoi(1:nlevgrnd) = my_zisoi(:)
+    col%nbedrock(bounds%begc:bounds%endc) = nbedrock_dummy_value
+
+    ! Set up hill_distance
+    l = bounds%begl
+    do c = lun%coli(l), lun%colf(l)
+       col%hill_distance(c) = real(c, kind=r8)
+    end do
+
+    deallocate(my_zisoi)
+
+  end subroutine setUp
+
+  subroutine tearDown(this)
+    ! clean up stuff set up in setup()
+    use clm_varcon, only: clm_varcon_clean
+    class(TestInit), intent(inout) :: this
+
+    call unittest_subgrid_teardown()
+    call clm_varcon_clean()
+
+  end subroutine tearDown
+
+  @Test
+  subroutine test_HillslopeSoilThicknessProfile_linear(this)
+    class(TestInit), intent(inout) :: this
+    integer, allocatable :: nbedrock_expected(:)
+    integer :: l, c
+
+    l = bounds%begl
+
+    col%active(bounds%begc:bounds%endc) = .true.
+    col%is_hillslope_column(bounds%begc:bounds%endc) = .true.
+
+    ! Get expected values
+    ! Layer 1 soil_depth_col = 8.5
+    ! Layer 2 soil_depth_col = 2.0
+    allocate(nbedrock_expected(bounds%begc:bounds%endc))
+    nbedrock_expected(lun%coli(l)) = 9999
+    nbedrock_expected(lun%coli(l) + 1) = 3
+
+    call HillslopeSoilThicknessProfile_linear(bounds, soil_depth_lowland, soil_depth_upland)
+
+    @assertEqual(nbedrock_expected(lun%coli(l):lun%colf(l)), col%nbedrock(lun%coli(l):lun%colf(l)))
+
+    deallocate(nbedrock_expected)
+
+  end subroutine test_HillslopeSoilThicknessProfile_linear
+
+  @Test
+  subroutine test_HillslopeSoilThicknessProfile_linear_inactive(this)
+    class(TestInit), intent(inout) :: this
+    integer, allocatable :: nbedrock_expected(:)
+    integer :: l, c
+
+    l = bounds%begl
+
+    col%active(bounds%begc:bounds%endc) = .false.
+    col%is_hillslope_column(bounds%begc:bounds%endc) = .true.
+
+    ! Get expected values
+    ! Layer 1 soil_depth_col = 8.5
+    ! Layer 2 soil_depth_col = 2.0, but not active
+    allocate(nbedrock_expected(bounds%begc:bounds%endc))
+    nbedrock_expected(lun%coli(l)) = 9999
+    nbedrock_expected(lun%coli(l) + 1) = 9999
+
+    call HillslopeSoilThicknessProfile_linear(bounds, soil_depth_lowland, soil_depth_upland)
+
+    @assertEqual(nbedrock_expected(lun%coli(l):lun%colf(l)), col%nbedrock(lun%coli(l):lun%colf(l)))
+
+    deallocate(nbedrock_expected)
+
+  end subroutine test_HillslopeSoilThicknessProfile_linear_inactive
+
+  @Test
+  subroutine test_HillslopeSoilThicknessProfile_linear_nohillslope(this)
+    class(TestInit), intent(inout) :: this
+    integer, allocatable :: nbedrock_expected(:)
+    integer :: l, c
+
+    l = bounds%begl
+
+    col%active(bounds%begc:bounds%endc) = .true.
+    col%is_hillslope_column(bounds%begc:bounds%endc) = .false.
+
+    ! Get expected values
+    ! Layer 1 soil_depth_col = 8.5
+    ! Layer 2 soil_depth_col = 2.0, but not is_hillslope_column
+    allocate(nbedrock_expected(bounds%begc:bounds%endc))
+    nbedrock_expected(lun%coli(l)) = 9999
+    nbedrock_expected(lun%coli(l) + 1) = 9999
+
+    call HillslopeSoilThicknessProfile_linear(bounds, soil_depth_lowland, soil_depth_upland)
+
+    @assertEqual(nbedrock_expected(lun%coli(l):lun%colf(l)), col%nbedrock(lun%coli(l):lun%colf(l)))
+
+    deallocate(nbedrock_expected)
+
+  end subroutine test_HillslopeSoilThicknessProfile_linear_nohillslope
+
+end module test_hillslopehydrologyUtils