import numpy as np CtoK = 273.16 # 273.15 """ Conversion factor for $^\circ\,$C to K """ kappa = 0.4 # NOTE: 0.41 """ von Karman's constant """ #------------------------------------------------------------------------------ def get_heights(h, dim_len): """ Reads input heights for velocity, temperature and humidity Parameters ---------- h : float input heights (m) dim_len : int length dimension Returns ------- hh : array """ hh = np.zeros((3, dim_len)) if (type(h) == float or type(h) == int): hh[0, :], hh[1, :], hh[2, :] = h, h, h elif (len(h) == 2 and np.ndim(h) == 1): hh[0, :], hh[1, :], hh[2, :] = h[0], h[1], h[1] elif (len(h) == 3 and np.ndim(h) == 1): hh[0, :], hh[1, :], hh[2, :] = h[0], h[1], h[2] elif (len(h) == 1 and np.ndim(h) == 2): hh = np.zeros((3, h.shape[1])) hh[0, :], hh[1, :], hh[2, :] = h[0, :], h[0, :], h[0, :] elif (len(h) == 2 and np.ndim(h) == 2): hh = np.zeros((3, h.shape[1])) hh[0, :], hh[1, :], hh[2, :] = h[0, :], h[1, :], h[1, :] elif (len(h) == 3 and np.ndim(h) == 2): hh = np.zeros((3, h.shape[1])) hh = np.copy(h) return hh # --------------------------------------------------------------------- def gc(lat, lon=None): """ Computes gravity relative to latitude Parameters ---------- lat : float latitude ($^\circ$) lon : float longitude ($^\circ$, optional) Returns ------- gc : float gravity constant (m/s^2) """ gamma = 9.7803267715 c1 = 0.0052790414 c2 = 0.0000232718 c3 = 0.0000001262 c4 = 0.0000000007 if lon is not None: lon_m, lat_m = np.meshgrid(lon, lat) else: lat_m = lat phi = lat_m*np.pi/180. xx = np.sin(phi) gc = (gamma*(1+c1*np.power(xx, 2)+c2*np.power(xx, 4)+c3*np.power(xx, 6) + c4*np.power(xx, 8))) return gc # --------------------------------------------------------------------- def visc_air(T): """ Computes the kinematic viscosity of dry air as a function of air temp. following Andreas (1989), CRREL Report 89-11. Parameters ---------- Ta : float air temperature ($^\circ$\,C) Returns ------- visa : float kinematic viscosity (m^2/s) """ T = np.asarray(T) if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius T = T-273.16 visa = 1.326e-5*(1+6.542e-3*T+8.301e-6*np.power(T, 2) - 4.84e-9*np.power(T, 3)) return visa # --------------------------------------------------------------------- def set_flag(miss, rh, u10n, q10n, t10n, Rb, hin, monob, itera, out=0): """ Set general flags. Parameters ---------- miss : int mask of missing input points rh : float relative humidity [%] u10n : float 10m neutral wind speed [ms^{-1}] q10n : float 10m neutral specific humidity [kg/kg] t10n : float 10m neutral air temperature [K] Rb : float bulk Richardson number hin : float measurement heights [m] monob : float Monin-Obukhov length [m] itera : int number of iteration out : int, optional output option for non converged points. The default is 0. Returns ------- flag : str """ # set maximum/minimum acceptable values u10max = 200 q10max = 40*0.001 t10min, t10max = 173, 373 Rbmin, Rbmax = -0.5, 0.2 flag = np.full(miss.shape, "n", dtype="object") flag = np.where(np.isnan(miss), "m", flag) # relative humidity flag flag = np.where(rh > 100, "r", flag) # u10n flag flag = np.where(((u10n < 0) | (u10n > u10max)) & (flag == "n"), "u", np.where(((u10n < 0) | (u10n > u10max)) & (np.char.find(flag.astype(str), 'u') == -1), flag+[","]+["u"], flag)) # q10n flag flag = np.where(((q10n < 0) | (q10n > q10max)) & (flag == "n"), "q", np.where(((q10n < 0) | (q10n > q10max)) & (flag != "n"), flag+[","]+["q"], flag)) # t10n flag flag = np.where(((t10n < t10min) | (t10n > t10max)) & (flag == "n"), "t", np.where( ((t10n < t10min) | (t10n > t10max)) & (flag != "n"), flag+[","]+["t"], flag)) # stability flag flag = np.where(((Rb < Rbmin) | (Rb > Rbmax) | ((hin[0]/monob) > 1000)) & (flag == "n"), "l", np.where(((Rb < Rbmin) | (Rb > Rbmax) | ((hin[0]/monob) > 1000)) & (flag != "n"), flag+[","]+["l"], flag)) if out == 1: flag = np.where((itera == -1) & (flag == "n"), "i", np.where( (itera == -1) & ((flag != "n") & ( np.char.find(flag.astype(str), 'm') == -1)), flag+[","]+["i"], flag)) else: flag = np.where((itera == -1) & (flag == "n"), "i", np.where( (itera == -1) & ((flag != "n") & ( np.char.find(flag.astype(str), 'm') == -1) & (np.char.find(flag.astype(str), 'u') == -1)), flag+[","]+["i"], flag)) return flag # --------------------------------------------------------------------- def get_outvars(out_var, cskin, gust): if out_var is None: # full output if cskin == 1 and gust[0] == 0: # skin ON and gust OFF res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n", "ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr", "usr", "psim", "psit", "psiq", "psim_ref", "psit_ref", "psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq", "uref", "tref", "qref", "dter", "dqer", "dtwl", "tkt", "qair", "qsea", "Rl", "Rs", "Rnl", "Rb", "rh", "rho", "cp", "lv", "theta", "itera") elif cskin == 0 and gust[0] != 0: # skin OFF and gust ON res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n", "ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr", "usr", "psim", "psit", "psiq", "psim_ref", "psit_ref", "psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq", "uref", "tref", "qref", "qair", "qsea", "ug", "usrGF", "GustFact", "Rb", "rh", "rho", "cp", "lv", "theta", "itera") else: res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n", "ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr", "usr", "psim", "psit", "psiq", "psim_ref", "psit_ref", "psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq", "uref", "tref", "qref", "dter", "dqer", "dtwl", "tkt", "qair", "qsea", "Rl", "Rs", "Rnl", "ug", "usrGF", "GustFact", "Rb", "rh", "rho", "cp", "lv", "theta", "itera") elif out_var == "limited": res_vars = ("tau", "sensible", "latent", "uref", "tref", "qref") else: res_vars = out_var return res_vars