Skip to content

GitLab

Commit b4c11efd authored by sbiri's avatar sbiri
Browse files

Update AirSeaFluxCode.py, flux_subs.py, get_init.py, hum_subs.py files

parent 3d302aa5
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 161 additions and 210 deletions
AirSeaFluxCode.py
View file @ b4c11efd
...@@ -71,10 +71,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -71,10 +71,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
set 1 to keep points set 1 to keep points
L : int L : int
Monin-Obukhov length definition options Monin-Obukhov length definition options
0 : default for S80, S88, LP82, YT96 and LY04 "S80" : default for S80, S88, LP82, YT96 and LY04
1 : following UA (Zeng et al., 1998), default for UA "ERA5" : following ERA5 (IFS Documentation cy46r1), default for ERA5
2 : following ERA5 (IFS Documentation cy46r1), default for ERA5
3 : COARE3.5 (Edson et al., 2013), default for C30, C35 and C40
Returns Returns
------- -------
res : array that contains res : array that contains
...@@ -113,6 +111,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -113,6 +111,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
""" """
logging.basicConfig(filename='flux_calc.log', logging.basicConfig(filename='flux_calc.log',
format='%(asctime)s %(message)s',level=logging.INFO) format='%(asctime)s %(message)s',level=logging.INFO)
logging.captureWarnings(True)
# check input values and set defaults where appropriate
lat, P, Rl, Rs, cskin, gust, tol, L = get_init(spd, T, SST, lat, P, Rl, Rs, lat, P, Rl, Rs, cskin, gust, tol, L = get_init(spd, T, SST, lat, P, Rl, Rs,
cskin, gust, L, tol, meth, cskin, gust, L, tol, meth,
qmeth) qmeth)
...@@ -123,7 +123,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -123,7 +123,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
' 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.nanmedian(P), np.nanmedian(Rl),
np.nanmedian(Rs), gust, cskin, L) np.nanmedian(Rs), gust, cskin, L)
#### # 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),
np.copy(T)*np.power(1000/P,287.1/1004.67)) # potential T np.copy(T)*np.power(1000/P,287.1/1004.67)) # potential T
...@@ -135,9 +135,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -135,9 +135,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
np.nanmedian(qsea), np.nanmedian(qair)) np.nanmedian(qsea), np.nanmedian(qair))
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")
# first guesses
dt = Ta - sst dt = Ta - sst
dq = qair - qsea dq = qair - qsea
# first guesses
t10n, q10n = np.copy(Ta), np.copy(qair) t10n, q10n = np.copy(Ta), np.copy(qair)
tv10n = t10n*(1 + 0.61*q10n) tv10n = t10n*(1 + 0.61*q10n)
# Zeng et al. 1998 # Zeng et al. 1998
...@@ -162,6 +162,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -162,6 +162,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
Rnl = 0.97*(5.67e-8*np.power(sst-0.3*cskin+CtoK, 4)-Rl) Rnl = 0.97*(5.67e-8*np.power(sst-0.3*cskin+CtoK, 4)-Rl)
dter = np.ones(T.shape)*0.3 dter = np.ones(T.shape)*0.3
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2)) dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
# gustiness adjustment
if (gust[0] == 1 and meth == "UA"): if (gust[0] == 1 and meth == "UA"):
wind = np.where(dtv >= 0, np.where(spd > 0.1, spd, 0.1), wind = np.where(dtv >= 0, np.where(spd > 0.1, spd, 0.1),
np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))) np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2)))
...@@ -169,6 +170,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -169,6 +170,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2)) wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
elif (gust[0] == 0): elif (gust[0] == 0):
wind = np.copy(spd) wind = np.copy(spd)
# stars and roughness lengths
usr = np.sqrt(cd*np.power(wind, 2)) usr = np.sqrt(cd*np.power(wind, 2))
zo = 0.0001*np.ones(spd.shape) zo = 0.0001*np.ones(spd.shape)
zot, zoq = 0.0001*np.ones(spd.shape), 0.0001*np.ones(spd.shape) zot, zoq = 0.0001*np.ones(spd.shape), 0.0001*np.ones(spd.shape)
...@@ -177,11 +179,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -177,11 +179,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
psit_calc(h_in[1]/monob, meth)) psit_calc(h_in[1]/monob, meth))
qsr = (dq+dqer*cskin)*kappa/(np.log(h_in[2]/zoq) - qsr = (dq+dqer*cskin)*kappa/(np.log(h_in[2]/zoq) -
psit_calc(h_in[2]/monob, meth)) psit_calc(h_in[2]/monob, meth))
# set-up to feed into iteration loop
it, ind = 0, np.where(spd > 0) it, ind = 0, np.where(spd > 0)
ii, itera = True, np.zeros(spd.shape)*np.nan ii, itera = True, np.zeros(spd.shape)*np.nan
tau = 0.05*np.ones(spd.shape) tau = 0.05*np.ones(spd.shape)
sensible = 5*np.ones(spd.shape) sensible = 5*np.ones(spd.shape)
latent = 65*np.ones(spd.shape) latent = 65*np.ones(spd.shape)
# iteration loop
while np.any(ii): while np.any(ii):
it += 1 it += 1
if it > n: if it > n:
...@@ -288,14 +292,14 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -288,14 +292,14 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
d = np.fabs(new-old) d = np.fabs(new-old)
if (tol[0] == 'flux'): if (tol[0] == 'flux'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) + ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3])) (d[2, :] > tol[3]))
elif (tol[0] == 'ref'): elif (tol[0] == 'ref'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) + ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3])) (d[2, :] > tol[3]))
elif (tol[0] == 'all'): elif (tol[0] == 'all'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) + ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3])+(d[3, :] > tol[4]) + (d[2, :] > tol[3])+(d[3, :] > tol[4]) +
(d[4, :] > tol[5])+(d[5, :] > tol[6])) (d[4, :] > tol[5])+(d[5, :] > tol[6]))
if (ind[0].size == 0): if (ind[0].size == 0):
ii = False ii = False
else: else:
......
flux_subs.py
View file @ b4c11efd
...@@ -651,10 +651,8 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -651,10 +651,8 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
---------- ----------
L : int L : int
Monin-Obukhov length definition options Monin-Obukhov length definition options
0 : default for S80, S88, LP82, YT96 and LY04 "S80" : default for S80, S88, LP82, YT96 and LY04
1 : following UA (Zeng et al., 1998), default for UA "ERA5" : following ERA5 (IFS Documentation cy46r1), default for ERA5
2 : following ERA5 (IFS Documentation cy46r1), default for ERA5
3 : COARE3.5 (Edson et al., 2013), default for C30, C35 and C40
lat : float lat : float
latitude latitude
usr : float usr : float
...@@ -702,19 +700,11 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -702,19 +700,11 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
""" """
g = gc(lat) g = gc(lat)
if (L == 0): if (L == "S80"):
tsrv = tsr+0.61*t10n*qsr tsrv = tsr+0.61*t10n*qsr
monob = ((tv10n*np.power(usr, 2))/(g*kappa*tsrv)) monob = ((tv10n*np.power(usr, 2))/(g*kappa*tsrv))
monob = np.where(np.fabs(monob) < 1, np.where(monob < 0, -1, 1), monob) monob = np.where(np.fabs(monob) < 1, np.where(monob < 0, -1, 1), monob)
elif (L == 1): elif (L == "ERA5"):
Rb = g*h_in[0]*dtv/(tv*np.power(wind, 2))
zol = np.where(Rb >= 0, Rb*np.log(h_in[0]/zo) /
(1-5*np.where(Rb < 0.19, Rb, 0.19)),
Rb*np.log(h_in[0]/zo))
monob = h_in[0]/zol
tsrv = tsr*(1.+0.61*qair)+0.61*th*qsr
# monob = ((tv*np.power(usr, 2))/(kappa*g*tsrv))
elif (L == 2):
tsrv = tsr+0.61*t10n*qsr tsrv = tsr+0.61*t10n*qsr
Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_in[0])+0.61*dq)) / Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_in[0])+0.61*dq)) /
np.power(wind, 2)) np.power(wind, 2))
...@@ -727,11 +717,6 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -727,11 +717,6 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
psit_calc((h_in[0]+zo)/monob, meth) + psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth)))) psit_calc(zot/monob, meth))))
monob = h_in[0]/zol monob = h_in[0]/zol
elif (L == 3):
tsrv = tsr+0.61*(T+CtoK)*qsr
zol = (kappa*g*h_in[0]/(T+CtoK)*(tsr+0.61*(T+CtoK)*qsr) /
np.power(usr, 2))
monob = h_in[0]/zol
return tsrv, monob return tsrv, monob
#------------------------------------------------------------------------------ #------------------------------------------------------------------------------
...@@ -784,48 +769,50 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq, ...@@ -784,48 +769,50 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq,
""" """
if (meth == "UA"): if (meth == "UA"):
usr = np.where(h_in[0]/monob < -1.574, kappa*wind / usr = np.where(h_in[0]/monob <= -1.574, kappa*wind /
(np.log(-1.574*monob/zo)-psim_calc(-1.574, meth) + (np.log(-1.574*monob/zo)-psim_calc(-1.574, meth) +
psim_calc(zo/monob, meth) + psim_calc(zo/monob, meth) +
1.14*(np.power(-h_in[0]/monob, 1/3) - 1.14*(np.power(-h_in[0]/monob, 1/3) -
np.power(1.574, 1/3))), np.power(1.574, 1/3))),
np.where((h_in[0]/monob > -1.574) & (h_in[0]/monob < 0), np.where(h_in[0]/monob < 0, kappa*wind /
kappa*wind/(np.log(h_in[0]/zo) - (np.log(h_in[0]/zo) -
psim_calc(h_in[0]/monob, meth) + psim_calc(h_in[0]/monob, meth) +
psim_calc(zo/monob, meth)), psim_calc(zo/monob, meth)),
np.where((h_in[0]/monob > 0) & np.where(h_in[0]/monob <= 1, kappa*wind /
(h_in[0]/monob < 1), (np.log(h_in[0]/zo) +
kappa*wind/(np.log(h_in[0]/zo) + 5*h_in[0]/monob-5*zo/monob),
5*h_in[0]/monob-5*zo/monob), kappa*wind/(np.log(monob/zo)+5 -
kappa*wind/(np.log(monob/zo)+5-5*zo/monob + 5*zo/monob +
5*np.log(h_in[0]/monob)+h_in[0]/monob-1)))) 5*np.log(h_in[0]/monob) +
h_in[0]/monob-1))))
# Zeng et al. 1998 (7-10) # Zeng et al. 1998 (7-10)
tsr = np.where(h_in[1]/monob < -0.465, kappa*(dt+dter*cskin) / tsr = np.where(h_in[1]/monob < -0.465, kappa*(dt+dter*cskin) /
(np.log((-0.465*monob)/zot) - (np.log((-0.465*monob)/zot) -
psit_calc(-0.465, meth)+0.8*(np.power(0.465, -1/3) - psit_calc(-0.465, meth)+0.8*(np.power(0.465, -1/3) -
np.power(-h_in[1]/monob, -1/3))), np.power(-h_in[1]/monob, -1/3))),
np.where((h_in[1]/monob > -0.465) & (h_in[1]/monob < 0), np.where(h_in[1]/monob < 0, kappa*(dt+dter*cskin) /
kappa*(dt+dter*cskin)/(np.log(h_in[1]/zot) - (np.log(h_in[1]/zot) -
psit_calc(h_in[1]/monob, meth) + psit_calc(h_in[1]/monob, meth) +
psit_calc(zot/monob, meth)), psit_calc(zot/monob, meth)),
np.where((h_in[1]/monob > 0) & (h_in[1]/monob < 1), np.where(h_in[1]/monob <= 1,
kappa*(dt+dter*cskin)/(np.log(h_in[1]/zot) + kappa*(dt+dter*cskin) /
5*h_in[1]/monob-5*zot/monob), (np.log(h_in[1]/zot) +
kappa*(dt+dter*cskin)/(np.log(monob/zot)+5 - 5*h_in[1]/monob-5*zot/monob),
5*zot/monob+5*np.log(h_in[1]/monob) + kappa*(dt+dter*cskin) /
h_in[1]/monob-1)))) (np.log(monob/zot)+5 -
5*zot/monob+5*np.log(h_in[1]/monob) +
h_in[1]/monob-1))))
# Zeng et al. 1998 (11-14) # Zeng et al. 1998 (11-14)
qsr = np.where(h_in[2]/monob < -0.465, kappa*(dq+dqer*cskin) / qsr = np.where(h_in[2]/monob < -0.465, kappa*(dq+dqer*cskin) /
(np.log((-0.465*monob)/zoq) - (np.log((-0.465*monob)/zoq) -
psit_calc(-0.465, meth)+psit_calc(zoq/monob, meth) + psit_calc(-0.465, meth)+psit_calc(zoq/monob, meth) +
0.8*(np.power(0.465, -1/3) - 0.8*(np.power(0.465, -1/3) -
np.power(-h_in[2]/monob, -1/3))), np.power(-h_in[2]/monob, -1/3))),
np.where((h_in[2]/monob > -0.465) & (h_in[2]/monob < 0), np.where(h_in[2]/monob < 0,
kappa*(dq+dqer*cskin)/(np.log(h_in[1]/zot) - kappa*(dq+dqer*cskin)/(np.log(h_in[1]/zot) -
psit_calc(h_in[2]/monob, meth) + psit_calc(h_in[2]/monob, meth) +
psit_calc(zoq/monob, meth)), psit_calc(zoq/monob, meth)),
np.where((h_in[2]/monob > 0) & np.where(h_in[2]/monob <= 1,
(h_in[2]/monob<1),
kappa*(dq+dqer*cskin) / kappa*(dq+dqer*cskin) /
(np.log(h_in[1]/zoq)+5*h_in[2]/monob - (np.log(h_in[1]/zoq)+5*h_in[2]/monob -
5*zoq/monob), 5*zoq/monob),
......
get_init.py
View file @ b4c11efd
...@@ -99,11 +99,11 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth): ...@@ -99,11 +99,11 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
if (np.all(Rs == None) or np.all(np.isnan(Rs))): if (np.all(Rs == None) or np.all(np.isnan(Rs))):
Rs = np.ones(spd.shape)*150 # set to default for COARE3.5 Rs = np.ones(spd.shape)*150 # set to default for COARE3.5
if ((cskin == None) and (meth == "S80" or meth == "S88" or meth == "LP82" if ((cskin == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96" or meth == "LY04")): or meth == "YT96" or meth == "UA" or
meth == "LY04")):
cskin = 0 cskin = 0
elif ((cskin == None) and (meth == "UA" or meth == "C30" elif ((cskin == None) and (meth == "C30" or meth == "C35" or meth == "C40"
or meth == "C35" or meth == "C40" or meth == "ERA5")):
or meth == "ERA5")):
cskin = 1 cskin = 1
if ((gust == None) and (meth == "C30" or meth == "C35" or meth == "C40")): if ((gust == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
gust = [1, 1.2, 600] gust = [1, 1.2, 600]
...@@ -113,19 +113,17 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth): ...@@ -113,19 +113,17 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
gust = [1, 1.2, 800] gust = [1, 1.2, 800]
elif (np.size(gust) < 3): elif (np.size(gust) < 3):
sys.exit("gust input must be a 3x1 array") sys.exit("gust input must be a 3x1 array")
if (L not in [None, 0, 1, 2, 3]): if (L not in [None, "S80", "ERA5"]):
sys.exit("L input must be either None, 0, 1, 2 or 3") sys.exit("L input must be either None, 0, 1, 2 or 3")
if ((L == None) and (meth == "S80" or meth == "S88" or meth == "LP82" if ((L == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96" or meth == "LY04")): or meth == "YT96" or meth == "LY04" or
L = 0 meth == "UA" or meth == "C30" or meth == "C35"
elif ((L == None) and (meth == "UA")): or meth == "C40")):
L = 1 L = "S80"
elif ((L == None) and (meth == "ERA5")): elif ((L == None) and (meth == "ERA5")):
L = 2 L = "ERA5"
elif ((L == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
L = 3
if (tol == None): if (tol == None):
tol = ['flux', 0.01, 1, 1] tol = ['flux', 0.01, 1, 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, P, Rl, Rs, cskin, gust, tol, L return lat, P, Rl, Rs, cskin, gust, tol, L
\ No newline at end of file
hum_subs.py
View file @ b4c11efd
This diff is collapsed.
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