Commit 2ab01484 authored by sbiri's avatar sbiri
Browse files

- flag "u" set relative to u10n

- for first three iterations u10n<0 is set to 0.5 to allow convergence for small negative values
- added initial restriction for n not to be less than 5; if so it is set to 5
parent 91b4cd7d
...@@ -71,8 +71,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -71,8 +71,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
option : 'flux' to set tolerance limits for fluxes only lim1-3 option : 'flux' to set tolerance limits for fluxes only lim1-3
option : 'ref' to set tolerance limits for height adjustment lim-1-3 option : 'ref' to set tolerance limits for height adjustment lim-1-3
option : 'all' to set tolerance limits for both fluxes and height option : 'all' to set tolerance limits for both fluxes and height
adjustment lim1-6 ['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1] adjustment lim1-6
default is tol=['flux', 1e-3, 0.1, 0.1] default is tol=['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1]
n : int n : int
number of iterations (defautl = 10) number of iterations (defautl = 10)
out : int out : int
...@@ -132,16 +132,19 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -132,16 +132,19 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
2021 / Author S. Biri 2021 / Author S. Biri
""" """
logging.basicConfig(filename='flux_calc.log', logging.basicConfig(filename='flux_calc.log', filemode="w",
format='%(asctime)s %(message)s',level=logging.INFO) format='%(asctime)s %(message)s',level=logging.INFO)
logging.captureWarnings(True) logging.captureWarnings(True)
# check input values and set defaults where appropriate # check input values and set defaults where appropriate
lat, hum, P, Rl, Rs, cskin, skin, wl, gust, tol, L = get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, tol, L, n = get_init(spd, T,
lat, hum, P, SST, lat,
hum, P,
Rl, Rs, Rl, Rs,
cskin, skin, cskin,
wl, gust, L, skin,
tol, meth, wl, gust,
L, tol, n,
meth,
qmeth) qmeth)
flag = np.ones(spd.shape, dtype="object")*"n" flag = np.ones(spd.shape, dtype="object")*"n"
flag = np.where(np.isnan(spd+T+SST+lat+hum[1]+P+Rs), "m", flag) flag = np.where(np.isnan(spd+T+SST+lat+hum[1]+P+Rs), "m", flag)
...@@ -150,8 +153,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -150,8 +153,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
h_out = get_heights(hout, 1) # desired height of output variables h_out = get_heights(hout, 1) # desired height of output variables
logging.info('method %s, inputs: lat: %s | P: %s | Rl: %s |' logging.info('method %s, inputs: lat: %s | P: %s | Rl: %s |'
' Rs: %s | gust: %s | cskin: %s | L : %s', meth, ' Rs: %s | gust: %s | cskin: %s | L : %s', meth,
np.nanmedian(lat), np.nanmedian(P), np.nanmedian(Rl), np.nanmedian(lat), np.round(np.nanmedian(P), 2),
np.nanmedian(Rs), gust, cskin, L) np.round(np.nanmedian(Rl),2 ), np.round(np.nanmedian(Rs), 2),
gust, cskin, L)
# set up/calculate temperatures and specific humidities # set up/calculate temperatures and specific humidities
th = np.where(T < 200, (np.copy(T)+CtoK) * th = np.where(T < 200, (np.copy(T)+CtoK) *
np.power(1000/P,287.1/1004.67), np.power(1000/P,287.1/1004.67),
...@@ -164,7 +168,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -164,7 +168,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
Ta = np.where(T < 200, np.copy(T)+CtoK+tlapse*h_in[1], Ta = np.where(T < 200, np.copy(T)+CtoK+tlapse*h_in[1],
np.copy(T)+tlapse*h_in[1]) # convert to Kelvin if needed np.copy(T)+tlapse*h_in[1]) # convert to Kelvin if needed
logging.info('method %s and q method %s | qsea:%s, qair:%s', meth, qmeth, logging.info('method %s and q method %s | qsea:%s, qair:%s', meth, qmeth,
np.nanmedian(qsea), np.nanmedian(qair)) np.round(np.nanmedian(qsea), 7),
np.round(np.nanmedian(qair), 7))
if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))): if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
print("qsea and qair cannot be nan") print("qsea and qair cannot be nan")
...@@ -322,10 +327,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -322,10 +327,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
tkt[ind] = 0.001*np.ones(T[ind].shape) tkt[ind] = 0.001*np.ones(T[ind].shape)
logging.info('method %s | dter = %s | dqer = %s | tkt = %s | Rnl = %s ' logging.info('method %s | dter = %s | dqer = %s | tkt = %s | Rnl = %s '
'| usr = %s | tsr = %s | qsr = %s', meth, '| usr = %s | tsr = %s | qsr = %s', meth,
np.nanmedian(dter), np.nanmedian(dqer), np.round(np.nanmedian(dter), 2),
np.nanmedian(tkt), np.nanmedian(Rnl), np.round(np.nanmedian(dqer), 7),
np.nanmedian(usr), np.nanmedian(tsr), np.round(np.nanmedian(tkt), 2),
np.nanmedian(qsr)) np.round(np.nanmedian(Rnl), 2),
np.round(np.nanmedian(usr), 3),
np.round(np.nanmedian(tsr), 4),
np.round(np.nanmedian(qsr), 7))
Rnl[ind] = 0.97*(5.67e-8*np.power(sst[ind] - Rnl[ind] = 0.97*(5.67e-8*np.power(sst[ind] -
dter[ind]*cskin, 4)-Rl[ind]) dter[ind]*cskin, 4)-Rl[ind])
t10n[ind] = (Ta[ind] - t10n[ind] = (Ta[ind] -
...@@ -362,6 +370,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -362,6 +370,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
wind[ind] = np.copy(spd[ind]) wind[ind] = np.copy(spd[ind])
u10n[ind] = wind[ind]-usr[ind]/kappa*(np.log(h_in[0, ind]/10) - u10n[ind] = wind[ind]-usr[ind]/kappa*(np.log(h_in[0, ind]/10) -
psim[ind]) psim[ind])
if (it < 4): # make sure you allow small negative values convergence
u10n = np.where(u10n < 0, 0.5, u10n)
flag = np.where((u10n < 0) & (flag == "n"), "u", flag = np.where((u10n < 0) & (flag == "n"), "u",
np.where((u10n < 0) & (flag != "u"), np.where((u10n < 0) & (flag != "u"),
flag+[","]+["u"], flag)) flag+[","]+["u"], flag))
...@@ -398,7 +408,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -398,7 +408,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
itera[ind] = -1 itera[ind] = -1
# itera = np.where(itera > n, -1, itera) # itera = np.where(itera > n, -1, itera)
logging.info('method %s | # of iterations:%s', meth, it) logging.info('method %s | # of iterations:%s', meth, it)
logging.info('method %s | # of points that did not converge :%s', meth, logging.info('method %s | # of points that did not converge :%s \n', meth,
ind[0].size) ind[0].size)
# calculate output parameters # calculate output parameters
rho = (0.34838*P)/(tv10n) rho = (0.34838*P)/(tv10n)
...@@ -428,30 +438,43 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -428,30 +438,43 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"], np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"],
flag)) flag))
flag = np.where(((Rb < -0.5) | (Rb > 0.2)) & (flag == "n"), "l", flag = np.where(((Rb < -0.5) | (Rb > 0.2)) & (flag == "n"), "l",
np.where(((Rb < -0.5) | (Rb > 0.2)) & (flag != "n"), np.where(((Rb < -0.5) | (Rb > 0.2)) &
flag+[","]+["l"], flag)) ((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)), flag+[","]+["l"], flag))
flag = np.where((itera == -1) & (flag == "n"), "i", flag = np.where((itera == -1) & (flag == "n"), "i",
np.where((itera == -1) & (flag != "n"), flag+[","]+["i"], np.where((itera == -1) &
flag)) ((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)),
flag+[","]+["i"], flag))
if (meth == "S80"): if (meth == "S80"):
flag = np.where(((spd < 6) | (spd > 22)) & (flag == "n"), "o", flag = np.where(((u10n < 6) | (u10n > 22)) & (flag == "n"), "o",
np.where(((spd < 6) | (spd > 22)) & (flag != "n"), np.where(((u10n < 6) | (u10n > 22)) &
((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)),
flag+[","]+["o"], flag)) flag+[","]+["o"], flag))
elif (meth == "LP82"): elif (meth == "LP82"):
flag = np.where(((spd < 3) | (spd > 25)) & (flag == "n"), "o", flag = np.where(((u10n < 3) | (u10n > 25)) & (flag == "n"), "o",
np.where(((spd < 3) | (spd > 25)) & (flag != "n"), np.where(((u10n < 3) | (u10n > 25)) &
((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)),
flag+[","]+["o"], flag)) flag+[","]+["o"], flag))
elif (meth == "YT96"): elif (meth == "YT96"):
flag = np.where(((spd < 3) | (spd > 26)) & (flag == "n"), "o", flag = np.where(((u10n < 3) | (u10n > 26)) & (flag == "n"), "o",
np.where(((spd < 3) | (spd > 26)) & (flag != "n"), np.where(((u10n < 3) | (u10n > 26)) &
((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)),
flag+[","]+["o"], flag)) flag+[","]+["o"], flag))
elif (meth == "UA"): elif (meth == "UA"):
flag = np.where((spd > 18) & (flag == "n"), "o", flag = np.where((u10n > 18) & (flag == "n"), "o",
np.where((spd > 18) & (flag != "n"), np.where((u10n > 18) &
((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)),
flag+[","]+["o"], flag)) flag+[","]+["o"], flag))
elif (meth == "LY04"): elif (meth == "LY04"):
flag = np.where((spd < 0.5) & (flag == "n"), "o", flag = np.where((u10n < 0.5) & (flag == "n"), "o",
np.where((spd < 0.5) & (flag != "n"), np.where((u10n < 0.5) &
((flag != "n") & (("u" in flag) == False) &
(("q" in flag) == False)),
flag+[","]+["o"], flag)) flag+[","]+["o"], flag))
if (hum == None): if (hum == None):
rh = np.ones(sst.shape)*80 rh = np.ones(sst.shape)*80
...@@ -462,8 +485,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -462,8 +485,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
Td = hum[1] # dew point temperature (K) Td = hum[1] # dew point temperature (K)
Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td)) Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
T = np.where(T < 200, np.copy(T)+CtoK, np.copy(T)) T = np.where(T < 200, np.copy(T)+CtoK, np.copy(T))
esd = 611.21*np.exp(17.502*((Td-273.16)/(Td-32.19))) esd = 611.21*np.exp(17.502*((Td-CtoK)/(Td-32.19)))
es = 611.21*np.exp(17.502*((T-273.16)/(T-32.19))) es = 611.21*np.exp(17.502*((T-CtoK)/(T-32.19)))
rh = 100*esd/es rh = 100*esd/es
rh = np.where(rh > 100, np.nan, rh) rh = np.where(rh > 100, np.nan, rh)
...@@ -511,17 +534,17 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, ...@@ -511,17 +534,17 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
if (out == 0): if (out == 0):
res[:, ind] = np.nan res[:, ind] = np.nan
# set missing values where data have non acceptable values # set missing values where data have non acceptable values
res = np.asarray([np.where((spd < 0) | (q10n < 0), np.nan, res = np.asarray([np.where(q10n < 0, np.nan,
res[i][:]) for i in range(39)]) res[i][:]) for i in range(39)])
# output with pandas # output with pandas
resAll = pd.DataFrame(data=res.T, index=range(len(spd)), resAll = pd.DataFrame(data=res.T, index=range(len(spd)),
columns=["tau", "shf", "lhf", "L", "cd", "cdn", "ct", columns=["tau", "shf", "lhf", "L", "cd", "cdn", "ct",
"ctn", "cq", "cqn", "tsrv", "tsr", "qsr", "ctn", "cq", "cqn", "tsrv", "tsr", "qsr",
"usr", "psim", "psit","psiq", "u10n", "t10n", "usr", "psim", "psit","psiq", "u10n",
"tv10n", "q10n", "zo", "zot", "zoq", "uref", "t10n", "tv10n", "q10n", "zo", "zot", "zoq",
"tref", "qref", "iteration", "dter", "dqer", "uref", "tref", "qref", "iteration", "dter",
"dtwl", "qair", "qsea", "Rl", "Rs", "Rnl", "dqer", "dtwl", "qair", "qsea", "Rl", "Rs",
"ug", "Rib", "rh"]) "Rnl", "ug", "Rib", "rh"])
resAll["flag"] = flag resAll["flag"] = flag
return resAll return resAll
...@@ -2,7 +2,7 @@ import numpy as np ...@@ -2,7 +2,7 @@ import numpy as np
import sys import sys
def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
meth, qmeth): n, meth, qmeth):
""" """
Checks initial input values and sets defaults if needed Checks initial input values and sets defaults if needed
...@@ -17,7 +17,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, ...@@ -17,7 +17,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
sea surface temperature in K sea surface temperature in K
lat : float lat : float
latitude (deg), default 45deg latitude (deg), default 45deg
hum : array hum : float
relative humidity, if None is set to 80% relative humidity, if None is set to 80%
P : float P : float
air pressure (hPa), default 1013hPa air pressure (hPa), default 1013hPa
...@@ -40,12 +40,14 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, ...@@ -40,12 +40,14 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
default else [1, 1.2, 800] default else [1, 1.2, 800]
L : int L : int
Monin-Obukhov length definition options Monin-Obukhov length definition options
tol : array tol : float
4x1 or 7x1 [option, lim1-3 or lim1-6] 4x1 or 7x1 [option, lim1-3 or lim1-6]
option : 'flux' to set tolerance limits for fluxes only lim1-3 option : 'flux' to set tolerance limits for fluxes only lim1-3
option : 'ref' to set tolerance limits for height adjustment lim-1-3 option : 'ref' to set tolerance limits for height adjustment lim-1-3
option : 'all' to set tolerance limits for both fluxes and height option : 'all' to set tolerance limits for both fluxes and height
adjustment lim1-6 ['all', 0.01, 0.01, 5e-05, 1e-3, 0.1, 0.1] adjustment lim1-6 ['all', 0.01, 0.01, 5e-05, 1e-3, 0.1, 0.1]
n : int
number of iterations
meth : str meth : str
"S80","S88","LP82","YT96","UA","LY04","C30","C35","ecmwf", "S80","S88","LP82","YT96","UA","LY04","C30","C35","ecmwf",
"Beljaars" "Beljaars"
...@@ -78,7 +80,8 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, ...@@ -78,7 +80,8 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
tolerance limits tolerance limits
L : int L : int
MO length switch MO length switch
n : int
number of iterations
""" """
# check if input is correct (type, size, value) and set defaults # check if input is correct (type, size, value) and set defaults
if ((type(spd) != np.ndarray) or (type(T) != np.ndarray) or if ((type(spd) != np.ndarray) or (type(T) != np.ndarray) or
...@@ -143,15 +146,13 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, ...@@ -143,15 +146,13 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
sys.exit("gust input must be a 3x1 array") sys.exit("gust input must be a 3x1 array")
if (L not in [None, "S80", "ecmwf"]): if (L not in [None, "S80", "ecmwf"]):
sys.exit("L input must be either None, S80 or ecmwf") sys.exit("L input must be either None, S80 or ecmwf")
if ((L == None) and (meth == "S80" or meth == "S88" or meth == "LP82" if (L == None):
or meth == "YT96" or meth == "LY04" or
meth == "UA" or meth == "C30" or meth == "C35"
or meth == "Beljaars")):
L = "S80"
elif ((L == None) and (meth == "ecmwf")):
L = "ecmwf" L = "ecmwf"
if (tol == None): if (tol == None):
tol = ['flux', 1e-3, 0.1, 0.1] tol = ['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1]
elif (tol[0] not in ['flux', 'ref', 'all']): elif (tol[0] not in ['flux', 'ref', 'all']):
sys.exit("unknown tolerance input") sys.exit("unknown tolerance input")
return lat, hum, P, Rl, Rs, cskin, skin, wl, gust, tol, L if (n < 5):
n = 5
print("Number of iterations should not be less than 5; n is set to 5")
return lat, hum, P, Rl, Rs, cskin, skin, wl, gust, tol, L, n
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment