Add rollover of surface fluxes (cdn) at high wind speed

[zarzycki]

Reopening this: ESMCI/cime#3149

From @vlarson: The current calculation of the neutral drag coefficient at 10 m (cdn) in subroutine shr_flux_atmOcn and subroutine shr_flux_atmOcn_diurnal does not include a moderating effect at wind speeds greater than 30 m/s. This rollover effect can be incorporated using equation (11a,b) of Large and Yeager (2009, Climate Dynamics).

From @zarzycki: Personally, I don't see any downside -- at lower surface wind speeds the LY09 equation collapses to what is currently in CIME. I have never done a gridpoint PDF of U10 at 1deg, but I can't fathom there being more than a handful of grid cells that ever see that wind speed so the solution should be effectively identical.

However, we have found benefits for tropical cyclones pressure-wind relationships at 0.25deg with the drag saturation (i.e., LY09).

FWIW, my patch is below (with some other Cd profiles I was playing around with as well...)

    qsat(Tk)   = 640380.0_R8 / exp(5107.4_R8/Tk)
-   cdn(Umps)  =   0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps
+!++CMZ
+! Large and Yeager 2009
+   cdn(Umps)  =  0.0027_R8 / min(33.0000_R8,Umps) + 0.000142_R8 + 0.0000764_R8 * min(33.0000_R8,Umps) - 3.14807e-13_r8 * min(33.0000_R8,Umps)**6
+! Capped Large and Pond by wind
+!   cdn(Umps)  =   0.0027_R8 / min(30.0_R8,Umps) + 0.000142_R8 + 0.0000764_R8 * min(30.0_R8,Umps) 
+! Capped Large and Pond by Cd
+!   cdn(Umps) = min(0.0025_R8, (0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps ))
+! Large and Pond
+!   cdn(Umps)  =   0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps
+!--CMZ
    psimhu(xd) = log((1.0_R8+xd*(2.0_R8+xd))*(1.0_R8+xd*xd)/8.0_R8) - 2.0_R8*atan(xd) + 1.571_R8
    psixhu(xd) = 2.0_R8 * log((1.0_R8 + xd*xd)/2.0_R8)

Relevant CMEPS lines: https://github.com/ESCOMP/CMEPS/blob/98dcf46c8886104b95cddfd5b02588b3dd9f6722/cesm/flux_atmocn/shr_flux_mod.F90#L262C1-L262C72