import numpy as np from util_subs import (kappa, visc_air) # --------------------------------------------------------------------- def cdn_calc(u10n, usr, Ta, grav, meth): """ Calculate neutral drag coefficient. Parameters ---------- u10n : float neutral 10m wind speed [m/s] usr : float friction velocity [m/s] Ta : float air temperature [K] grav : float gravity [m/s^2] meth : str Returns ------- cdn : float zo : float """ cdn = np.zeros(Ta.shape)*np.nan if meth == "S80": # eq. 14 Smith 1980 cdn = np.maximum((0.61+0.063*u10n)*0.001, (0.61+0.063*6)*0.001) elif meth == "LP82": # Large & Pond 1981 u10n <11m/s & eq. 21 Large & Pond 1982 cdn = np.where(u10n < 11, 1.2*0.001, (0.49+0.065*u10n)*0.001) elif meth in ["S88", "UA", "ecmwf", "C30", "C35", "Beljaars"]: cdn = cdn_from_roughness(u10n, usr, Ta, grav, meth) elif meth == "YT96": # convert usr in eq. 21 to cdn to expand for low wind speeds cdn = np.power((0.10038+u10n*2.17e-3+np.power(u10n, 2)*2.78e-3 - np.power(u10n, 3)*4.4e-5)/u10n, 2) elif meth == "NCAR": # eq. 11 Large and Yeager 2009 cdn = np.where(u10n > 0.5, (0.142+2.7/u10n+u10n/13.09 - 3.14807e-10*np.power(u10n, 6))*1e-3, (0.142+2.7/0.5+0.5/13.09 - 3.14807e-10*np.power(0.5, 6))*1e-3) cdn = np.where(u10n > 33, 2.34e-3, np.copy(cdn)) cdn = np.maximum(np.copy(cdn), 0.1e-3) else: raise ValueError("Unknown method cdn: "+meth) zo = 10/np.exp(kappa/np.sqrt(cdn)) return cdn, zo # --------------------------------------------------------------------- def cdn_from_roughness(u10n, usr, Ta, grav, meth): """ Calculate neutral drag coefficient from roughness length. Parameters ---------- u10n : float neutral 10m wind speed [m/s] usr : float friction velocity [m/s] Ta : float air temperature [K] grav : float [m/s] gravity meth : str Returns ------- cdn : float """ # cdn = (0.61+0.063*u10n)*0.001 zo, zc, zs = np.zeros(Ta.shape), np.zeros(Ta.shape), np.zeros(Ta.shape) for it in range(5): if meth == "S88": # Charnock roughness length (eq. 4 in Smith 88) zc = 0.011*np.power(usr, 2)/grav # smooth surface roughness length (eq. 6 in Smith 88) zs = 0.11*visc_air(Ta)/usr zo = zc + zs # eq. 7 & 8 in Smith 88 elif meth == "UA": # valid for 0 10, 0.011+(u10n-10)*(0.018-0.011)/(18-10), np.where(u10n > 18, 0.018, a)) zo = a*np.power(usr, 2)/grav+0.11*visc_air(Ta)/usr elif meth == "C35": # eq.6-11 Edson et al. (2013) zo = (0.11*visc_air(Ta)/usr + np.minimum(0.0017*19-0.0050, 0.0017*u10n-0.0050) * np.power(usr, 2)/grav) elif meth in ["ecmwf", "Beljaars"]: # eq. (3.26) p.38 over sea IFS Documentation cy46r1 zo = 0.018*np.power(usr, 2)/grav+0.11*visc_air(Ta)/usr else: raise ValueError("Unknown method for cdn_from_roughness "+meth) cdn = np.power(kappa/np.log(10/zo), 2) return cdn # --------------------------------------------------------------------- def cd_calc(cdn, hin, hout, psim): """ Calculate drag coefficient at reference height. Parameters ---------- cdn : float neutral drag coefficient hin : float wind speed height [m] hout : float reference height [m] psim : float momentum stability function Returns ------- cd : float """ cd = (cdn/np.power(1+(np.sqrt(cdn)*(np.log(hin/hout)-psim))/kappa, 2)) return cd # --------------------------------------------------------------------- def ctqn_calc(corq, zol, cdn, usr, zo, Ta, meth): """ Calculate neutral heat and moisture exchange coefficients. Parameters ---------- corq : flag to select "ct" or "cq" zol : float height over MO length cdn : float neutral drag coefficient usr : float friction velocity [m/s] zo : float surface roughness [m] Ta : float air temperature [K] meth : str Returns ------- ctqn : float neutral heat exchange coefficient zotq : float roughness length for t or q """ if meth in ["S80", "S88", "YT96"]: cqn = np.ones(Ta.shape)*1.20*0.001 # from S88 ctn = np.ones(Ta.shape)*1.00*0.001 zot = 10/(np.exp(np.power(kappa, 2) / (ctn*np.log(10/zo)))) zoq = 10/(np.exp(np.power(kappa, 2) / (cqn*np.log(10/zo)))) elif meth == "LP82": cqn = np.where((zol <= 0), 1.15*0.001, 1*0.001) ctn = np.where((zol <= 0), 1.13*0.001, 0.66*0.001) zot = 10/(np.exp(np.power(kappa, 2) / (ctn*np.log(10/zo)))) zoq = 10/(np.exp(np.power(kappa, 2) / (cqn*np.log(10/zo)))) elif meth == "NCAR": cqn = np.maximum(34.6*0.001*np.sqrt(cdn), 0.1e-3) ctn = np.maximum(np.where(zol <= 0, 32.7*0.001*np.sqrt(cdn), 18*0.001*np.sqrt(cdn)), 0.1e-3) zot = 10/(np.exp(np.power(kappa, 2) / (ctn*np.log(10/zo)))) zoq = 10/(np.exp(np.power(kappa, 2) / (cqn*np.log(10/zo)))) elif meth == "UA": # Zeng et al. 1998 (25) rr = usr*zo/visc_air(Ta) zoq = zo/np.exp(2.67*np.power(rr, 1/4)-2.57) zot = np.copy(zoq) cqn = np.power(kappa, 2)/(np.log(10/zo)*np.log(10/zoq)) ctn = np.power(kappa, 2)/(np.log(10/zo)*np.log(10/zoq)) elif meth == "C30": rr = zo*usr/visc_air(Ta) zoq = np.minimum(5e-5/np.power(rr, 0.6), 1.15e-4) # moisture roughness zot = np.copy(zoq) # temperature roughness cqn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zoq) ctn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zot) elif meth == "C35": rr = zo*usr/visc_air(Ta) zoq = np.minimum(5.8e-5/np.power(rr, 0.72), 1.6e-4) # moisture rough. zot = np.copy(zoq) # temperature roughness cqn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zoq) ctn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zot) elif meth in ["ecmwf", "Beljaars"]: # eq. (3.26) p.38 over sea IFS Documentation cy46r1 zot = 0.40*visc_air(Ta)/usr zoq = 0.62*visc_air(Ta)/usr cqn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zoq) ctn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zot) else: raise ValueError("Unknown method ctqn: "+meth) if corq == "ct": ctqn = ctn zotq = zot elif corq == "cq": ctqn = cqn zotq = zoq else: raise ValueError("Unknown flag - should be ct or cq: "+corq) return ctqn, zotq # --------------------------------------------------------------------- def ctq_calc(cdn, cd, ctqn, hin, hout, psitq): """ Calculate heat and moisture exchange coefficients at reference height. Parameters ---------- cdn : float neutral drag coefficient cd : float drag coefficient at reference height ctqn : float neutral heat or moisture exchange coefficient hin : float original temperature or humidity sensor height [m] hout : float reference height [m] psitq : float heat or moisture stability function Returns ------- ctq : float heat or moisture exchange coefficient """ ctq = (ctqn*np.sqrt(cd/cdn) / (1+ctqn*((np.log(hin/hout)-psitq)/(kappa*np.sqrt(cdn))))) return ctq # --------------------------------------------------------------------- def get_stabco(meth): r""" Give the coefficients $\alpha$, $\beta$, $\gamma$ for stability functions. Parameters ---------- meth : str Returns ------- coeffs : float """ alpha, beta, gamma = 0, 0, 0 if meth in ["S80", "S88", "NCAR", "UA", "ecmwf", "C30", "C35", "Beljaars"]: alpha, beta, gamma = 16, 0.25, 5 # Smith 1980, from Dyer (1974) elif meth == "LP82": alpha, beta, gamma = 16, 0.25, 7 elif meth == "YT96": alpha, beta, gamma = 20, 0.25, 5 else: raise ValueError("Unknown method stabco: "+meth) coeffs = np.zeros(3) coeffs[0] = alpha coeffs[1] = beta coeffs[2] = gamma return coeffs # --------------------------------------------------------------------- def psim_calc(zol, meth): """ Calculate momentum stability function. Parameters ---------- zol : float height over MO length meth : str Returns ------- psim : float """ if meth == "ecmwf": psim = psim_ecmwf(zol) elif meth in ["C30", "C35"]: psim = psiu_26(zol, meth) elif meth == "Beljaars": # Beljaars (1997) eq. 16, 17 psim = np.where(zol < 0, psim_conv(zol, meth), psi_Bel(zol)) else: psim = np.where(zol < 0, psim_conv(zol, meth), psim_stab(zol, meth)) return psim # --------------------------------------------------------------------- def psit_calc(zol, meth): """ Calculate heat stability function. Parameters ---------- zol : float height over MO length meth : str parameterisation method Returns ------- psit : float """ if meth == "ecmwf": psit = np.where(zol < 0, psi_conv(zol, meth), psi_ecmwf(zol)) elif meth in ["C30", "C35"]: psit = psit_26(zol) elif meth == "Beljaars": # Beljaars (1997) eq. 16, 17 psit = np.where(zol < 0, psi_conv(zol, meth), psi_Bel(zol)) else: psit = np.where(zol < 0, psi_conv(zol, meth), psi_stab(zol, meth)) return psit # --------------------------------------------------------------------- def psi_Bel(zol): """ Calculate momentum/heat stability function. Parameters ---------- zol : float height over MO length meth : str parameterisation method Returns ------- psit : float """ a, b, c, d = 0.7, 0.75, 5, 0.35 psi = -(a*zol+b*(zol-c/d)*np.exp(-d*zol)+b*c/d) return psi # --------------------------------------------------------------------- def psi_ecmwf(zol): """ Calculate heat stability function for stable conditions. For method ecmwf Parameters ---------- zol : float height over MO length Returns ------- psit : float """ # eq (3.22) p. 37 IFS Documentation cy46r1 a, b, c, d = 1, 2/3, 5, 0.35 psit = -b*(zol-c/d)*np.exp(-d*zol)-np.power(1+(2/3)*a*zol, 1.5)-(b*c)/d+1 return psit # --------------------------------------------------------------------- def psit_26(zol): """ Compute temperature structure function as in C35. Parameters ---------- zol : float height over MO length Returns ------- psi : float """ b, d = 2/3, 0.35 dzol = np.minimum(d*zol, 50) psi = -1*((1+b*zol)**1.5+b*(zol-14.28)*np.exp(-dzol)+8.525) k = np.where(zol < 0) x = np.sqrt(1-15*zol[k]) psik = 2*np.log((1+x)/2) x = np.power(1-34.15*zol[k], 1/3) psic = (1.5*np.log((1+x+np.power(x, 2))/3)-np.sqrt(3) * np.arctan((1+2*x)/np.sqrt(3))+4*np.arctan(1)/np.sqrt(3)) f = np.power(zol[k], 2)/(1+np.power(zol[k], 2)) psi[k] = (1-f)*psik+f*psic return psi # --------------------------------------------------------------------- def psi_conv(zol, meth): """ Calculate heat stability function for unstable conditions. Parameters ---------- zol : float height over MO length meth : str parameterisation method Returns ------- psit : float """ coeffs = get_stabco(meth) alpha, beta = coeffs[0], coeffs[1] xtmp = np.power(1-alpha*zol, beta) psit = 2*np.log((1+np.power(xtmp, 2))*0.5) return psit # --------------------------------------------------------------------- def psi_stab(zol, meth): """ Calculate heat stability function for stable conditions. Parameters ---------- zol : float height over MO length meth : str parameterisation method Returns ------- psit : float """ coeffs = get_stabco(meth) gamma = coeffs[2] psit = -gamma*zol return psit # --------------------------------------------------------------------- def psim_ecmwf(zol): """ Calculate momentum stability function for method ecmwf. Parameters ---------- zol : float height over MO length Returns ------- psim : float """ # eq (3.20, 3.22) p. 37 IFS Documentation cy46r1 coeffs = get_stabco("ecmwf") alpha, beta = coeffs[0], coeffs[1] xtmp = np.power(1-alpha*zol, beta) a, b, c, d = 1, 2/3, 5, 0.35 psim = np.where(zol < 0, np.pi/2-2*np.arctan(xtmp) + np.log((np.power(1+xtmp, 2)*(1+np.power(xtmp, 2)))/8), -b*(zol-c/d)*np.exp(-d*zol)-a*zol-(b*c)/d) return psim # --------------------------------------------------------------------- def psiu_26(zol, meth): """ Compute velocity structure function C35. Parameters ---------- zol : float height over MO length Returns ------- psi : float """ if meth == "C30": dzol = np.minimum(0.35*zol, 50) # stable psi = -1*((1+zol)+0.6667*(zol-14.28)*np.exp(-dzol)+8.525) k = np.where(zol < 0) # unstable x = (1-15*zol[k])**0.25 psik = (2*np.log((1+x)/2)+np.log((1+x*x)/2)-2*np.arctan(x) + 2*np.arctan(1)) x = (1-10.15*zol[k])**(1/3) psic = (1.5*np.log((1+x+x*x)/3) - np.sqrt(3)*np.arctan((1+2*x)/np.sqrt(3)) + 4*np.arctan(1)/np.sqrt(3)) f = zol[k]**2/(1+zol[k]**2) psi[k] = (1-f)*psik+f*psic elif meth == "C35": dzol = np.minimum(50, 0.35*zol) # stable a, b, c, d = 0.7, 3/4, 5, 0.35 psi = -1*(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d) k = np.where(zol < 0) # unstable x = np.power(1-15*zol[k], 1/4) psik = 2*np.log((1+x)/2)+np.log((1+x*x)/2) - \ 2*np.arctan(x)+2*np.arctan(1) x = np.power(1-10.15*zol[k], 1/3) psic = (1.5*np.log((1+x+np.power(x, 2))/3)-np.sqrt(3) * np.arctan((1+2*x)/np.sqrt(3))+4*np.arctan(1)/np.sqrt(3)) f = np.power(zol[k], 2)/(1+np.power(zol[k], 2)) psi[k] = (1-f)*psik+f*psic return psi # ---------------------------------------------------------------------------- def psim_conv(zol, meth): """ Calculate momentum stability function for unstable conditions. Parameters ---------- zol : float height over MO length meth : str parameterisation method Returns ------- psim : float """ coeffs = get_stabco(meth) alpha, beta = coeffs[0], coeffs[1] xtmp = np.power(1-alpha*zol, beta) psim = (2*np.log((1+xtmp)*0.5)+np.log((1+np.power(xtmp, 2))*0.5) - 2*np.arctan(xtmp)+np.pi/2) return psim # --------------------------------------------------------------------- def psim_stab(zol, meth): """ Calculate momentum stability function for stable conditions. Parameters ---------- zol : float height over MO length meth : str parameterisation method Returns ------- psim : float """ coeffs = get_stabco(meth) gamma = coeffs[2] psim = -gamma*zol return psim # --------------------------------------------------------------------- def get_gust(beta, Ta, usr, tsrv, zi, grav): """ Compute gustiness. Parameters ---------- beta : float constant Ta : float air temperature [K] usr : float friction velocity [m/s] tsrv : float star virtual temperature of air [K] zi : int scale height of the boundary layer depth [m] grav : float gravity Returns ------- ug : float [m/s] """ if np.nanmax(Ta) < 200: # convert to K if in Celsius Ta = Ta+273.16 # minus sign to allow cube root Bf = (-grav/Ta)*usr*tsrv ug = np.ones(np.shape(Ta))*0.2 ug = np.where(Bf > 0, beta*np.power(Bf*zi, 1/3), 0.2) return ug # --------------------------------------------------------------------- def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, cd, ct, cq, meth): """ Calculate star wind speed, temperature and specific humidity. Parameters ---------- hin : float sensor heights [m] monob : float M-O length [m] wind : float wind speed [m/s] zo : float momentum roughness length [m] zot : float temperature roughness length [m] zoq : float moisture roughness length [m] dt : float temperature difference [K] dq : float specific humidity difference [g/kg] ct : float temperature exchange coefficient cq : float moisture exchange coefficient meth : str bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA", "NCAR", "C30", "C35", "ecmwf", "Beljaars" Returns ------- usr : float friction wind speed [m/s] tsr : float star temperature [K] qsr : float star specific humidity [g/kg] """ if meth == "UA": # Zeng et al. 1998 # away from extremes UA follows e.g. S80 usr = wind*np.sqrt(cd) tsr = ct*wind*dt/usr qsr = cq*wind*dq/usr # momentum hol0 = hin[0]/np.copy(monob) # very unstable (Zeng et al. 1998 eq 7) usr = np.where( hol0 <= -1.574, wind*kappa/(np.log(-1.574*monob/zo) - psim_calc(-1.574, meth) + psim_calc(zo/monob, meth) + 1.14*(np.power(-hin[0]/monob, 1/3) - np.power(1.574, 1/3))), usr) # very stable (Zeng et al. 1998 eq 10) usr = np.where( hol0 > 1, wind*kappa/(np.log(monob/zo)+5-5*zo/monob + 5*np.log(hin[0]/monob)+hin[0]/monob-1), usr) # temperature hol1 = hin[1]/np.copy(monob) # very unstable (Zeng et al. 1998 eq 11) tsr = np.where( hol1 < -0.465, kappa*dt/(np.log((-0.465*monob)/zot) - psit_calc(-0.465, meth) + 0.8*(np.power(0.465, -1/3) - np.power(-hin[1]/monob, -1/3))), tsr) # very stable (Zeng et al. 1998 eq 14) tsr = np.where( hol1 > 1, kappa*(dt)/(np.log(monob/zot)+5-5*zot/monob + 5*np.log(hin[1]/monob)+hin[1]/monob-1), tsr) # humidity hol2 = hin[2]/monob # very unstable (Zeng et al. 1998 eq 11) qsr = np.where( hol2 < -0.465, kappa*dq/(np.log((-0.465*monob)/zoq) - psit_calc(-0.465, meth) + psit_calc(zoq/monob, meth) + 0.8*(np.power(0.465, -1/3) - np.power(-hin[2]/monob, -1/3))), qsr) # very stable (Zeng et al. 1998 eq 14) qsr = np.where(hol2 > 1, kappa*dq/(np.log(monob/zoq)+5-5*zoq/monob + 5*np.log(hin[2]/monob) + hin[2]/monob-1), qsr) else: usr = wind*np.sqrt(cd) tsr = ct*wind*dt/usr qsr = cq*wind*dq/usr return usr, tsr, qsr # --------------------------------------------------------------------- def get_tsrv(tsr, qsr, Ta, qair): """ Calculate virtual star temperature. Parameters ---------- tsr : float star temperature (K) qsr : float star specific humidity (g/kg) Ta : float air temperature (K) qair : float air specific humidity (g/kg) Returns ------- tsrv : float virtual star temperature (K) """ # as in aerobulk One_on_L in mod_phymbl.f90 tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr return tsrv # --------------------------------------------------------------------- def get_Rb(grav, usr, hin_u, hin_t, tv, dtv, wind, monob, meth): """ Calculate bulk Richardson number. Parameters ---------- grav : float acceleration due to gravity (m/s2) usr : float friction wind speed (m/s) hin_u : float u sensor height (m) hin_t : float t sensor height (m) tv : float virtual temperature (K) dtv : float virtual temperature difference, air and sea (K) wind : float wind speed (m/s) monob : float Monin-Obukhov length from previous iteration step (m) psim : float momentum stability function meth : str bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA", "NCAR", "C30", "C35", "ecmwf", "Beljaars" Returns ------- Rb : float Richardson number """ # now input dtv # tvs = sst*(1+0.6077*qsea) # virtual SST # dtv = tv - tvs # virtual air - sea temp. diff # adjust wind to t measurement height uz = (wind-usr/kappa*(np.log(hin_u/hin_t)-psim_calc(hin_u/monob, meth) + psim_calc(hin_t/monob, meth))) Rb = grav*dtv*hin_t/(tv*uz*uz) return Rb # --------------------------------------------------------------------- def get_LRb(Rb, hin_t, monob, zo, zot, meth): """ Calculate Monin-Obukhov length following ecmwf (IFS Documentation cy46r1). default for methods ecmwf and Beljaars Parameters ---------- Rb : float Richardson number hin_t : float t sensor height (m) monob : float Monin-Obukhov length from previous iteration step (m) zo : float surface roughness (m) zot : float temperature roughness length (m) meth : str bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA", "NCAR", "C30", "C35", "ecmwf", "Beljaars" Returns ------- monob : float M-O length (m) """ zol = Rb*(np.power( np.log((hin_t+zo)/zo)-psim_calc((hin_t+zo)/monob, meth) + psim_calc(zo/monob, meth), 2)/(np.log((hin_t+zo)/zot) - psit_calc((hin_t+zo)/monob, meth) + psit_calc(zot/monob, meth))) monob = hin_t/zol return monob # --------------------------------------------------------------------- def get_Ltsrv(tsrv, grav, tv, usr): """ Calculate Monin-Obukhov length from tsrv. Parameters ---------- tsrv : float virtual star temperature (K) grav : float acceleration due to gravity (m/s2) tv : float virtual temperature (K) usr : float friction wind speed (m/s) Returns ------- monob : float M-O length (m) """ # tmp = (g*kappa*tsrv / # np.maximum(np.power(usr, 2)*Ta*(1+0.6077*qair), 1e-9)) # tmp = np.minimum(np.abs(tmp), 200)*np.sign(tmp) # monob = 1/np.copy(tmp) tsrv = np.maximum(np.abs(tsrv), 1e-9)*np.sign(tsrv) monob = (np.power(usr, 2)*tv)/(grav*kappa*tsrv) return monob # ---------------------------------------------------------------------