Minor bug in creation of soil color in mksurfdata_map: points can be given the default soil color, when they should have a real color

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2012-01-29 20:29:57 -0700
Bugzilla Id: 1457
Bugzilla CC: mvertens@ucar.edu, rfisher@ucar.edu,

I think there is a bug in mksoilcol in module mksoilMod in mksurfdata_map... I haven't done any testing to confirm this, but the logic looks wrong to me:

The relevant code is copied below.

It looks like the intention is that, in the first do loop, wst(0,no) should be 0 if there is ever a non-zero color in this output cell. However, it seems there is an unintentional order dependence here: Consider the following two cases, showing the colors in input grid cells corresponding to a single output grid cell; assume for simplicity that all weights are 0.25:

(1) 0, 0, 0, 1

(2) 1, 0, 0, 0

In case (1), wst(0,no) will be correctly set to 0.0 when the 1 value is read. This will result in soil_color_o = 1.

In case (2), wst(0,no) will be set to 0.0 at first, but upon reading the three 0 values, wst(0,no) will become 0.75. This will result in soil_color_o = 0, which will later be changed to either 4 or 15.

So, if I'm correct in my interpretation of this algorithm, soil_color_o sometimes depends on the order in which values appear in the source grid cells!

 do n = 1,tgridmap%ns
    ni = tgridmap%src_indx(n)
    no = tgridmap%dst_indx(n)
    wt = tgridmap%wovr(n)
    k  = soil_color_i(ni) * tdomain%mask(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

    ! Rank non-zero weights by color type. wsti(1) is the most extensive
    ! color type. 

    if (color(no) == 1) then
       call mkrank (nsoicol, wst(0:nsoicol,no), miss, wsti, num)
       soil_color_o(no) = wsti(1)
    end if

    ! If land but no color, set color to 15 (in older dataset generic 
    ! soil color 4)
    
    if (nsoicol == 8) then
       if (soil_color_o(no)==0) soil_color_o(no) = 4
    else if (nsoicol == 20) then
       if (soil_color_o(no)==0) soil_color_o(no) = 15
    end if

[billsacks]

I have verified that this is, indeed, a bug by isolating the relevant logic and wrapping it in a unit test (74b69e3). However, this has absolutely no impact on a 1-degree surface dataset (./mksurfdata.pl -no-crop -glc_nec 10 -y 1850 -res 0.9x1.25). Now that I understand the nature of the bug, I am not surprised by that finding: I think this will only change the outcome in uncommon circumstances (maybe more common for high-resolution datasets?).