Add unit test exposing bug and fix it

Add a unit test that exposes the bug in #4 - in some cases,
points are given soil color 0 when they should be given a real
color. Also change code to fix this bug. (The new unit test fails
without the changes in mksoilUtilsMod.F90.)

(Addresses #4, but fully fixing that issue will require
recreating all surface datasets.)

tools/mksurfdata_map/src/mksoilUtilsMod.F90

@@ -55,7 +55,7 @@ subroutine dominant_soil_color(tgridmap, mask_i, soil_color_i, nsoicol, soil_col
     integer  :: k, n, ni, no, ns_i, ns_o
     real(r8) :: wt                            ! map overlap weight
     real(r8), allocatable :: wst(:,:)         ! overlap weights, by surface type
-    integer , allocatable :: color(:)         ! 0: none; 1: some
+    logical :: has_color                      ! whether this grid cell has non-zero color
     integer, parameter :: miss = 99999        ! missing data indicator
 
     character(len=*), parameter :: subname = 'dominant_soil_color'
@@ -74,8 +74,6 @@ subroutine dominant_soil_color(tgridmap, mask_i, soil_color_i, nsoicol, soil_col
     ns_o = size(soil_color_o)
     allocate(wst(0:nsoicol,ns_o))
     wst(0:nsoicol,:) = 0
-    allocate(color(ns_o))
-    color(:) = 0
 
     ! TODO: need to do a loop to determine
     ! the maximum number of over lap cells throughout the grid 
@@ -89,19 +87,24 @@ subroutine dominant_soil_color(tgridmap, mask_i, soil_color_i, nsoicol, soil_col
        wt = tgridmap%wovr(n)
        k  = soil_color_i(ni) * mask_i(ni)
        wst(k,no) = wst(k,no) + wt
-       if (k>0 .and. wst(k,no)>0.) then
-          color(no) = 1
-          wst(0,no) = 0.0
-       end if
     enddo
 
     soil_color_o(:) = 0
     do no = 1,ns_o
 
+       ! If the output cell has any non-zero-colored inputs, then set the weight of
+       ! zero-colored inputs to 0, to ensure that the zero-color is NOT dominant.
+       if (any(wst(1:nsoicol,no) > 0.)) then
+          has_color = .true.
+          wst(0,no) = 0.0
+       else
+          has_color = .false.
+       end if
+
        ! Rank non-zero weights by color type. wsti(1) is the most extensive
        ! color type. 
 
-       if (color(no) == 1) then
+       if (has_color) then
           call mkrank (nsoicol, wst(0:nsoicol,no), miss, wsti, num)
           soil_color_o(no) = wsti(1)
        end if
@@ -133,7 +136,6 @@ subroutine dominant_soil_color(tgridmap, mask_i, soil_color_i, nsoicol, soil_col
     end do
 
     deallocate (wst)
-    deallocate (color)
 
   end subroutine dominant_soil_color
 

tools/mksurfdata_map/src/test/mksoilUtils_test/test_dominant_soil_color.pf

@@ -119,4 +119,22 @@ contains
     @assertEqual(1, soil_color_o(1))
   end subroutine noColorInFirstPoints
 
+  @Test
+  subroutine noColorInLastPoints(this)
+    ! No color in the last points, but a color in the first point
+    class(tdsc), intent(inout) :: this
+    integer :: mask_i(4)
+    integer :: soil_color_i(4)
+    integer :: soil_color_o(1)
+
+    call this%createGridmap1dst([0.25_r8, 0.25_r8, 0.25_r8, 0.25_r8])
+    ! Some points are inside the mask with color = 0, other points are outside the mask
+    mask_i(:) = 1
+    soil_color_i(:) = [1, 0, 0, 0]
+
+    call dominant_soil_color(this%gridmap, mask_i, soil_color_i, 20, soil_color_o)
+
+    @assertEqual(1, soil_color_o(1))
+  end subroutine noColorInLastPoints
+
 end module test_dominant_soil_color