Merge branch 'tillage' into tillage-and-residues3

src/soilbiogeochem/TillageMod.F90

@@ -247,31 +247,33 @@ subroutine get_tillage_multipliers(tillage_mults, idop)
   end subroutine get_tillage_multipliers
 
 
-  function get_fraction_tilled(zisoi, j, max_tillage_depth_gft) result(fraction_tilled)
+  function get_fraction_tilled(layer_bottom, layer_thickness, max_tillage_depth_gft) result(fraction_tilled)
     ! !ARGUMENTS
-    ! zisoi and max_tillage_depth are passed to this function, rather than using from clm_varcon or
-    ! parent scope (respectively), to enable unit testing.
-    real(r8) :: zisoi(:) ! Soil layer interface depths
-    integer  :: j        ! Soil layer
+    real(r8), intent(in) :: layer_bottom    ! Soil layer interface (between j and j+1) depth (zisoi)
+    real(r8), intent(in) :: layer_thickness ! Soil layer thickness (dzsoi_decomp)
     real(r8) :: max_tillage_depth_gft ! Maximum tillage depth
     ! !LOCAL VARIABLES
-    real(r8) :: layer_top ! Depth (m) of the top of this soil layer. zisoi is the depth of the bottom.
-    real(r8) :: layer_thickness ! Thickness of this soil layer (m)
+    real(r8) :: layer_top
     ! !RESULT
     real(r8) :: fraction_tilled ! Fraction of this layer that's within the tillage depth
 
-    if (j == 1) then
-        layer_top = 0._r8
-    else
-        layer_top = zisoi(j-1)
-    end if
-
+    ! If the top of the layer is below the max tillage depth, do not till.
+    layer_top = layer_bottom - layer_thickness
     if (layer_top > max_tillage_depth_gft) then
         fraction_tilled = 0._r8
         return
     end if
 
-    layer_thickness = zisoi(j) - layer_top
+    ! Handle zero-thickness layers. This may not be necessary.
+    if (layer_thickness == 0._r8) then
+        if (layer_bottom <= max_tillage_depth_gft) then
+            fraction_tilled = 1._r8
+        else
+            fraction_tilled = 0._r8
+        end if
+        return
+    end if
+
     fraction_tilled = max(0._r8, min(1._r8, (max_tillage_depth_gft - layer_top) / layer_thickness))
 
   end function get_fraction_tilled
@@ -285,7 +287,7 @@ subroutine get_apply_tillage_multipliers(idop, c, j, decomp_k)
     !
     ! !USES
     use pftconMod , only : nc3crop
-    use clm_varcon, only : zisoi
+    use clm_varcon, only : zisoi, dzsoi_decomp
     use landunit_varcon , only : istcrop
     use PatchType , only : patch
     !
@@ -303,7 +305,7 @@ subroutine get_apply_tillage_multipliers(idop, c, j, decomp_k)
     real(r8) :: fraction_tilled ! Fraction of this layer that's within the tillage depth
 
     ! Skip tillage if column is inactive or this layer doesn't get tilled
-    fraction_tilled = get_fraction_tilled(zisoi, j, max_tillage_depth)
+    fraction_tilled = get_fraction_tilled(zisoi(j), dzsoi_decomp(j), max_tillage_depth)
     if (.not. col%active(c) .or. fraction_tilled == 0._r8 .or. col%lun_itype(c) /= istcrop) then
         return
     end if

src/soilbiogeochem/test/tillage_test/test_tillage.pf

@@ -17,30 +17,38 @@ contains
   @Test
   subroutine test_get_fraction_tilled()
 
-    real(r8) :: zisoi(5)          ! Depth of soil interfaces (bottom of each layer)
+    integer, parameter :: nlayers = 5
+    real(r8) :: zisoi(nlayers)        ! Depth of soil interfaces (bottom of each layer)
+    real(r8) :: dzsoi_decomp(nlayers) ! Thickness of soil layers
     integer  :: j                 ! Soil layer
-    real(r8) :: max_tillage_depth
+    real(r8), parameter :: max_tillage_depth = 0.21_r8
 
     zisoi = [0.01_r8, 0.05_r8, 0.1_r8, 0.5_r8, 1._r8]
-
-    max_tillage_depth = 10._r8
-    do j = 1, 5
-      @assertEqual(1._r8, get_fraction_tilled(zisoi, j, max_tillage_depth))
-    end do
-
-    max_tillage_depth = 0._r8
-    do j = 1, 5
-      @assertEqual(0._r8, get_fraction_tilled(zisoi, j, max_tillage_depth))
-    end do
-
-    max_tillage_depth = 0.21_r8
-    @assertEqual(1._r8,    get_fraction_tilled(zisoi, 1, max_tillage_depth))
-    @assertEqual(1._r8,    get_fraction_tilled(zisoi, 2, max_tillage_depth))
-    @assertEqual(1._r8,    get_fraction_tilled(zisoi, 3, max_tillage_depth))
-    @assertEqual(0.275_r8, get_fraction_tilled(zisoi, 4, max_tillage_depth), tolerance=tol)
-    @assertEqual(0._r8,    get_fraction_tilled(zisoi, 5, max_tillage_depth))
+    dzsoi_decomp = [zisoi(1) - 0._r8, &
+                    zisoi(2) - zisoi(1), &
+                    zisoi(3) - zisoi(2), &
+                    zisoi(4) - zisoi(3), &
+                    zisoi(5) - zisoi(4)]
+
+    @assertEqual(1._r8,    get_fraction_tilled(zisoi(1), dzsoi_decomp(1), max_tillage_depth))
+    @assertEqual(1._r8,    get_fraction_tilled(zisoi(2), dzsoi_decomp(2), max_tillage_depth))
+    @assertEqual(1._r8,    get_fraction_tilled(zisoi(3), dzsoi_decomp(3), max_tillage_depth))
+    @assertEqual(0.275_r8, get_fraction_tilled(zisoi(4), dzsoi_decomp(4), max_tillage_depth), tolerance=tol)
+    @assertEqual(0._r8,    get_fraction_tilled(zisoi(5), dzsoi_decomp(5), max_tillage_depth))
 
   end subroutine test_get_fraction_tilled
 
 
+  @Test
+
+  subroutine test_get_fraction_tilled_0thickness()
+    real(r8), parameter :: max_tillage_depth = 0.5_r8
+
+    @assertEqual(1._r8, get_fraction_tilled(0.4_r8, 0._r8, max_tillage_depth))
+    @assertEqual(1._r8, get_fraction_tilled(0.5_r8, 0._r8, max_tillage_depth))
+    @assertEqual(0._r8, get_fraction_tilled(0.6_r8, 0._r8, max_tillage_depth))
+  end subroutine test_get_fraction_tilled_0thickness
+
+
+
 end module test_tillage