diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py index 068dbfc52dcf0c54acd25c1407dd2933223afe54..1c6287f2dc037f30f27feb4bd56881b2a7976a8d 100644 --- a/AirSeaFluxCode.py +++ b/AirSeaFluxCode.py @@ -422,6 +422,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, tref = tref-(CtoK+tlapse*h_out[1]) qref = (qair-qsr/kappa*(np.log(h_in[2]/h_out[2]) - psit+psit_calc(h_out[2]/monob, meth))) + if (wl == 0): + dtwl = np.zeros(T.shape) # reset to zero if not used flag = np.where((q10n < 0) & (flag == "n"), "q", np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"], flag)) @@ -444,8 +446,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, np.where(((u10n < 3) | (u10n > 26)) & (flag != "n"), flag+[","]+["o"], flag)) elif (meth == "UA"): - flag = np.where(((u10n < 0.5) | (u10n > 18)) & (flag == "n"), "o", - np.where(((u10n < 0.5) | (u10n > 18)) & (flag != "n"), + flag = np.where((u10n > 18) & (flag == "n"), "o", + np.where((u10n > 18) & (flag != "n"), flag+[","]+["o"], flag)) elif (meth == "LY04"): flag = np.where((u10n < 0.5) & (flag == "n"), "o", diff --git a/flux_subs.py b/flux_subs.py index f8e079bc7f81f46d8207abc59796acd1efea9915..3175abd5f51564f5cbc7c010c8ff6fc9752323f8 100755 --- a/flux_subs.py +++ b/flux_subs.py @@ -33,7 +33,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"): np.where((u10n < 11) & (u10n >= 4), 1.2*0.001, (0.49+0.065*u10n)*0.001)) elif (meth == "S88" or meth == "UA" or meth == "ecmwf" or meth == "C30" or - meth == "C35" or meth == "Beljaars"): + meth == "C35" or meth == "Beljaars"): cdn = cdn_from_roughness(u10n, Ta, None, lat, meth) elif (meth == "YT96"): # for u<3 YT96 convert usr in eq. 21 to cdn @@ -96,9 +96,9 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"): np.where(u10n > 18, 0.018, a)) zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr elif (meth == "C35"): - a = 0.011*np.ones(Ta.shape) - a = np.where(u10n > 19, 0.0017*19-0.0050, 0.0017*u10n-0.0050) - zo = 0.11*visc_air(Ta)/usr+a*np.power(usr, 2)/g + zo = (0.11*visc_air(Ta)/usr + + np.minimum(0.0017*19-0.0050, 0.0017*u10n-0.0050) * + np.power(usr, 2)/g) elif ((meth == "ecmwf" or meth == "Beljaars")): # eq. (3.26) p.38 over sea IFS Documentation cy46r1 zo = 0.018*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr @@ -174,9 +174,9 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): re=usr*zo/visc_air(Ta) zoq = zo/np.exp(2.67*np.power(re, 1/4)-2.57) zot = zoq - cqn = np.where((u10n > 0.5) & (u10n < 18), np.power(kappa, 2) / + cqn = np.where(u10n < 18, np.power(kappa, 2) / (np.log(10/zo)*np.log(10/zoq)), np.nan) - ctn = np.where((u10n > 0.5) & (u10n < 18), np.power(kappa, 2) / + ctn = np.where(u10n < 18, np.power(kappa, 2) / (np.log(10/zo)*np.log(10/zoq)), np.nan) elif (meth == "C30"): usr = np.sqrt(cdn*np.power(u10n, 2))