Another divide-by-zero (caught by CESM in debug mode)

[mnlevy1981]

In cvmix_kpp.F90, when running MatchBoth and computing the derivative of the shape function at the boundary layer depth, one of the terms is

second_term = real(5,cvmix_r8)*surf_buoy/(surf_fric**4)

This term is only included in derivative computation if the column is stable, but it is computed regardless (and when surf_fric=0 debug mode aborts because of the divide by zero). Our proposed fix was to change this to

if (lstable) then
  second_term = real(5,cvmix_r8)*surf_buoy/(surf_fric**4)
else
  second_term = cvmix_zero
end if

But we probably also want a check to make sure surf_fric.ne.cvmix_zero in the lstable=.true. case.

[mnlevy1981]

Note that with the fix above, we no longer need the if (lstable) checks when computing d[MTS]shapeAt1

if (lstable) &
  dMshapeAt1 = dMshapeAt1 + second_term*Mdiff_OBL

if (lstable) then
  dTshapeAt1 = dTshapeAt1 + second_term*Tdiff_OBL
  dSshapeAt1 = dSshapeAt1 + second_term*Sdiff_OBL
end if

because in unstable cases we would just be adding 0 to each term (I assume the logical check will be more of a performance hit than the multiply and add once we do a better job of vectorizing our data structures)