Skip to content

GitLab

Commit a10a41da authored by sbiri's avatar sbiri
Browse files

1. humidity option changed to strings ('rh', 'q', 'Td') instead of integer (0, 1, 2) 

2. added break if method not in options
3. added break if qmeth not in options
4. added break if gust not a 3x1 array
5. added break if tol option is not in options
6. added restriction if abs(L) < 1
parent 30625bc5
Hide whitespace changes
Inline Side-by-side
Showing with 49 additions and 31 deletions
AirSeaFluxCode.py
View file @ a10a41da
...@@ -7,9 +7,9 @@ from flux_subs import (kappa, CtoK, get_heights, cdn_calc, cd_calc, get_skin, ...@@ -7,9 +7,9 @@ from flux_subs import (kappa, CtoK, get_heights, cdn_calc, cd_calc, get_skin,
def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
hin=18, hout=10, Rl=None, Rs=None, jcool=None, hin=18, hout=10, Rl=None, Rs=None, jcool=None,
gust=None, meth="S80", qmeth="Buck2", tol=None, n=10, gust=None, meth="S80", qmeth="Buck2", tol=None, n=10,
out=0): out=0):
""" Calculates momentum and heat fluxes using different parameterizations """ Calculates momentum and heat fluxes using different parameterizations
Parameters Parameters
...@@ -25,9 +25,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -25,9 +25,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
latitude (deg), default 45deg latitude (deg), default 45deg
hum : float hum : float
humidity input switch 2x1 [x, values] default is relative humidity humidity input switch 2x1 [x, values] default is relative humidity
x=0 : relative humidity in % x='rh' : relative humidity in %
x=1 : specific humidity (g/kg) x='q' : specific humidity (g/kg)
x=2 : dew point temperature (K) x='Td' : dew point temperature (K)
P : float P : float
air pressure (hPa), default 1013hPa air pressure (hPa), default 1013hPa
hin : float hin : float
...@@ -110,31 +110,38 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -110,31 +110,38 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
h_in = get_heights(hin, len(spd)) # heights of input measurements/fields h_in = get_heights(hin, len(spd)) # heights of input measurements/fields
h_out = get_heights(hout, 1) # desired height of output variables h_out = get_heights(hout, 1) # desired height of output variables
# if input values are nan break # if input values are nan break
if meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
"C40","ERA5"]:
sys.exit("unknown method")
if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler", "CIMO",
"GoffGratch", "MagnusTetens", "Buck", "Buck2", "WMO",
"WMO2000", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
sys.exit("unknown q-method")
if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))): if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))):
sys.exit("input wind, T or SST is empty") sys.exit("input wind, T or SST is empty")
logging.debug('all spd or T or SST input is nan') logging.debug('all spd or T or SST input is nan')
if (np.all(lat) == None): # set latitude to 45deg if empty if (np.all(lat == None)): # set latitude to 45deg if empty
lat = 45*np.ones(spd.shape) lat = 45*np.ones(spd.shape)
elif ((np.all(lat) != None) and (np.size(lat) == 1)): elif ((np.all(lat != None)) and (np.size(lat) == 1)):
lat = np.ones(spd.shape)*np.copy(lat) lat = np.ones(spd.shape)*np.copy(lat)
g = gc(lat, None) # acceleration due to gravity g = gc(lat, None) # acceleration due to gravity
if ((np.all(P) == None) and (meth == "C30" or meth == "C40")): if ((np.all(P == None)) and (meth == "C30" or meth == "C40")):
P = np.ones(spd.shape)*1015 # if empty set to default for COARE3.0 P = np.ones(spd.shape)*1015 # if empty set to default for COARE3.0
elif ((np.all(P) == None) or np.all(np.isnan(P))): elif ((np.all(P == None)) or np.all(np.isnan(P))):
P = np.ones(spd.shape)*1013 P = np.ones(spd.shape)*1013
logging.debug('input P is empty and set to 1013hPa') logging.debug('input P is empty and set to 1013hPa')
elif (((np.all(P) != None) or np.all(~np.isnan(P))) and np.size(P) == 1): elif (((np.all(P != None)) or np.all(~np.isnan(P))) and np.size(P) == 1):
P = np.ones(spd.shape)*np.copy(P) P = np.ones(spd.shape)*np.copy(P)
if ((np.all(Rl) == None or np.all(np.isnan(Rl))) and meth == "C30"): if ((np.all(Rl == None) or np.all(np.isnan(Rl))) and meth == "C30"):
Rl = np.ones(spd.shape)*150 # set to default for COARE3.0 Rl = np.ones(spd.shape)*150 # set to default for COARE3.0
elif (((np.all(Rl) == None or np.all(np.isnan(Rl))) and meth == "C35") or elif (((np.all(Rl == None) or np.all(np.isnan(Rl))) and meth == "C35") or
((np.all(Rl) == None or np.all(np.isnan(Rl))) and meth == "C40")): ((np.all(Rl == None) or np.all(np.isnan(Rl))) and meth == "C40")):
Rl = np.ones(spd.shape)*370 # set to default for COARE3.5 Rl = np.ones(spd.shape)*370 # set to default for COARE3.5
elif (np.all(Rl) == None or np.all(np.isnan(Rl))): elif (np.all(Rl == None) or np.all(np.isnan(Rl))):
Rl = np.ones(spd.shape)*370 # set to default for COARE3.5 Rl = np.ones(spd.shape)*370 # set to default for COARE3.5
if ((np.all(Rs) == None or np.all(np.isnan(Rs))) and meth == "C30"): if ((np.all(Rs == None) or np.all(np.isnan(Rs))) and meth == "C30"):
Rs = np.ones(spd.shape)*370 # set to default for COARE3.0 Rs = np.ones(spd.shape)*370 # set to default for COARE3.0
elif (np.all(Rs) == None or np.all(np.isnan(Rs))): elif (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 ((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]
...@@ -142,17 +149,21 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -142,17 +149,21 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
gust = [1, 1, 1000] gust = [1, 1, 1000]
elif (gust == None): elif (gust == None):
gust = [1, 1.2, 800] gust = [1, 1.2, 800]
elif (np.size(gust) < 3):
sys.exit("gust input must be a 3x1 array")
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']):
sys.exit("unknown tolerance input")
if ((jcool == None) and (meth == "S80" or meth == "S88" or meth == "LP82" if ((jcool == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96")): or meth == "YT96")):
jcool = 0 jcool = 0
elif ((jcool == None) and (meth == "UA" or meth == "LY04" or meth == "C30" elif ((jcool == None) and (meth == "UA" or meth == "LY04" or meth == "C30"
or meth == "C35" or meth == "C40" or meth == "C35" or meth == "C40"
or meth == "ERA5")): or meth == "ERA5")):
jcool = 1 jcool = 1
logging.info('method %s, inputs: h_in: %s | lat: %s | P: %s | Rl: %s |' logging.info('method %s, inputs: lat: %s | P: %s | Rl: %s |'
' Rs: %s | gust: %s | jcool: %s', meth, h_in[:, 1], ' Rs: %s | gust: %s | jcool: %s', meth,
np.nanmedian(lat), np.nanmedian(P), np.nanmedian(Rl), np.nanmedian(lat), np.nanmedian(P), np.nanmedian(Rl),
np.nanmedian(Rs), gust, jcool) np.nanmedian(Rs), gust, jcool)
#### ####
...@@ -166,16 +177,18 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -166,16 +177,18 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
RH = np.ones(spd.shape)*80 RH = np.ones(spd.shape)*80
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg) qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg) qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg)
elif (hum[0] == 0): elif (hum[0] not in ['rh', 'q', 'Td']):
sys.exit("unknown humidity input")
elif (hum[0] == 'rh'):
RH = hum[1] RH = hum[1]
if (np.all(RH) < 1): if (np.all(RH < 1)):
sys.exit("input relative humidity units should be \%") sys.exit("input relative humidity units should be \%")
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg) qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg) qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg)
elif (hum[0] == 1): elif (hum[0] == 'q'):
qair = hum[1] qair = hum[1]
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg) qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
elif (hum[0] == 2): elif (hum[0] == 'Td'):
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))
...@@ -255,7 +268,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -255,7 +268,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
psit_calc(zot/monob, meth)))) psit_calc(zot/monob, meth))))
monob = h_in[0]/zol monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, ' logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, Rb:%s, monob:%s', meth, 'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo), np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob)) np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
elif (meth == "C30" or meth == "C35" or meth == "C40"): elif (meth == "C30" or meth == "C35" or meth == "C40"):
...@@ -411,8 +424,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -411,8 +424,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
else: else:
usr[ind] = (wind[ind]*kappa/(np.log(h_in[0, ind]/zo[ind]) - usr[ind] = (wind[ind]*kappa/(np.log(h_in[0, ind]/zo[ind]) -
psim_calc(h_in[0, ind]/monob[ind], meth))) psim_calc(h_in[0, ind]/monob[ind], meth)))
tsr[ind] = ct[ind]*wind[ind]*(dt[ind]+dter[ind]*jcool)/usr[ind]
qsr[ind] = cq[ind]*wind[ind]*(dq[ind]+dqer[ind]*jcool)/usr[ind] qsr[ind] = cq[ind]*wind[ind]*(dq[ind]+dqer[ind]*jcool)/usr[ind]
tsr[ind] = ct[ind]*wind[ind]*(dt[ind]+dter[ind]*jcool)/usr[ind]
if (jcool == 1): if (jcool == 1):
dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind], dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
rho[ind], Rl[ind], rho[ind], Rl[ind],
...@@ -458,7 +471,11 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -458,7 +471,11 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
monob[ind] = h_in[0, ind]/zol[ind] monob[ind] = h_in[0, ind]/zol[ind]
else: else:
tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind] tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
monob[ind] = (tv10n[ind]*usr[ind]**2)/(g[ind]*kappa*tsrv[ind]) monob[ind] = ((tv10n[ind]*np.power(usr[ind], 2)) /
(g[ind]*kappa*tsrv[ind]))
monob[ind] = np.where(np.fabs(monob[ind]) < 1,
np.where(monob[ind] < 0, -1, 1),
monob[ind])
psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth) psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth) psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth) psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
...@@ -518,9 +535,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -518,9 +535,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
# calculate output parameters # calculate output parameters
rho = (0.34838*P)/(tv10n) rho = (0.34838*P)/(tv10n)
t10n = t10n-(273.16+tlapse*ref_ht) t10n = t10n-(273.16+tlapse*ref_ht)
# zo = ref_ht/np.exp(kappa/cdn**0.5) zo = ref_ht/np.exp(kappa/cdn**0.5)
# zot = ref_ht/(np.exp(kappa**2/(ctn*np.log(ref_ht/zo)))) zot = ref_ht/(np.exp(kappa**2/(ctn*np.log(ref_ht/zo))))
# zoq = ref_ht/(np.exp(kappa**2/(cqn*np.log(ref_ht/zo)))) zoq = ref_ht/(np.exp(kappa**2/(cqn*np.log(ref_ht/zo))))
uref = (spd-usr/kappa*(np.log(h_in[0]/h_out[0])-psim + uref = (spd-usr/kappa*(np.log(h_in[0]/h_out[0])-psim +
psim_calc(h_out[0]/monob, meth))) psim_calc(h_out[0]/monob, meth)))
tref = (Ta-tsr/kappa*(np.log(h_in[1]/h_out[1])-psit + tref = (Ta-tsr/kappa*(np.log(h_in[1]/h_out[1])-psit +
...@@ -561,5 +578,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, ...@@ -561,5 +578,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
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.where(spd < 0, np.nan, res) res = np.where(spd < 0, np.nan, res)
return res, ind return res, ind
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