diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index b5cd378151da7291bb19078432c8ef95c8b0c3fa..62b29f713fcc662c5c37d5c9c68ba6c6d72633b6 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -6,20 +6,14 @@ from flux_subs import (kappa, CtoK, get_heights, cdn_calc, cd_calc, get_skin,
gc, qsat_sea, qsat_air, visc_air, psit_26, psiu_26)
-def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
- gust, meth, qmeth, n):
+def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
+ hin=18, hout=10, Rl=None, Rs=None, jcool=None,
+ gust=None, meth="S80", qmeth="Buck2", tol=None, n=10,
+ out=0):
""" Calculates momentum and heat fluxes using different parameterizations
Parameters
----------
- meth : str
- "S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
- qmeth : str
- is the saturation evaporation method to use amongst
- "HylandWexler","Hardy","Preining","Wexler","GoffGratch","CIMO",
- "MagnusTetens","Buck","Buck2","WMO","WMO2000","Sonntag","Bolton",
- "IAPWS","MurphyKoop"]
- default is Buck2
spd : float
relative wind speed in m/s (is assumed as magnitude difference
between wind and surface current vectors)
@@ -28,28 +22,52 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
SST : float
sea surface temperature in K (will convert if < 200)
lat : float
- latitude (deg)
- RH : float
- relative humidity in %
+ latitude (deg), default 45deg
+ hum : float
+ humidity input switch 2x1 [x, values] default is relative humidity
+ x=0 : relative humidity in %
+ x=1 : specific humidity (g/kg)
+ x=2 : dew point temperature (K)
P : float
- air pressure
+ air pressure (hPa), default 1013hPa
hin : float
- sensor heights in m (array 3x1 or 3xn) default 10m ([10, 10, 10])
+ sensor heights in m (array 3x1 or 3xn), default 18m
hout : float
output height, default is 10m
- zi : int
- PBL height (m)
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
jcool : int
- 0 if sst is true ocean skin temperature, else 1
+ 0 switch cool skin adjustment off, else 1
+ default is 1
gust : int
- 1 to include the effect of gustiness, else 0
+ 3x1 [x, beta, zi] x=1 to include the effect of gustiness, else 0
+ beta gustiness parameter, beta=1 for UA, beta=1.2 for COARE
+ zi PBL height (m) 600 for COARE, 1000 for UA and ERA5, 800 default
+ default for COARE [1, 1.2, 600]
+ default for UA, ERA5 [1, 1, 1000]
+ default else [1, 1.2, 800]
+ meth : str
+ "S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
+ qmeth : str
+ is the saturation evaporation method to use amongst
+ "HylandWexler","Hardy","Preining","Wexler","GoffGratch","CIMO",
+ "MagnusTetens","Buck","Buck2","WMO","WMO2000","Sonntag","Bolton",
+ "IAPWS","MurphyKoop"]
+ default is Buck2
+ tol : float
+ 4x1 or 7x1 [option, lim1-3 or lim1-6]
+ option : 'flux' to set tolerance limits for fluxes only lim1-3
+ option : 'ref' to set tolerance limits for height adjustment lim-1-3
+ option : 'all' to set tolerance limits for both fluxes and height
+ adjustment lim1-6 ['all', 0.01, 0.01, 5e-05, 0.01, 1, 1]
+ default is tol=['flux', 0.01, 1, 1]
n : int
- number of iterations
-
+ number of iterations (defautl = 10)
+ out : int
+ set 0 to set points that have not converged to missing (default)
+ set 1 to keep points
Returns
-------
res : array that contains
@@ -77,9 +95,9 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
22. surface roughness length (zo)
23. heat roughness length (zot)
24. moisture roughness length (zoq)
- 25. velocity at reference height (urefs)
- 26. temperature at reference height (trefs)
- 27. specific humidity at reference height (qrefs)
+ 25. velocity at reference height (uref)
+ 26. temperature at reference height (tref)
+ 27. specific humidity at reference height (qref)
28. number of iterations until convergence
ind : int
the indices in the matrix for the points that did not converge
@@ -91,43 +109,52 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
ref_ht, tlapse = 10, 0.0098 # reference height, lapse rate
h_in = get_heights(hin, len(spd)) # heights of input measurements/fields
h_out = get_heights(hout, 1) # desired height of output variables
- if np.all(np.isnan(lat)): # set latitude to 45deg if empty
- lat=45*np.ones(spd.shape)
- g = gc(lat, None) # acceleration due to gravity
# if input values are nan break
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")
- logging.debug('all input is nan')
- if (np.all(np.isnan(RH)) and meth == "C35"):
- RH = np.ones(spd.shape)*80 # if empty set to default for COARE3.5
- elif (np.all(np.isnan(RH))):
- sys.exit("input RH is empty")
- logging.debug('input RH is empty')
- if (np.all(np.isnan(P)) and (meth == "C30" or meth == "C40")):
+ logging.debug('all spd or T or SST input is nan')
+ if (np.all(lat) == None): # set latitude to 45deg if empty
+ lat = 45*np.ones(spd.shape)
+ elif ((np.all(lat) != None) and (np.size(lat) == 1)):
+ lat = np.ones(spd.shape)*np.copy(lat)
+ g = gc(lat, None) # acceleration due to gravity
+ 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
- elif ((np.all(np.isnan(P))) and meth == "C35"):
- P = np.ones(spd.shape)*1013 # if empty set to default for COARE3.5
- elif (np.all(np.isnan(P))):
- sys.exit("input P is empty")
- logging.debug('input P is empty')
- if (np.all(np.isnan(Rl)) and meth == "C30"):
+ elif ((np.all(P) == None) or np.all(np.isnan(P))):
+ P = np.ones(spd.shape)*1013
+ 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):
+ P = np.ones(spd.shape)*np.copy(P)
+ 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
- elif ((np.all(np.isnan(Rl)) and meth == "C35") or
- (np.all(np.isnan(Rl)) and meth == "C40")):
+ 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")):
Rl = np.ones(spd.shape)*370 # set to default for COARE3.5
- elif (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
- if (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
- elif ((np.all(np.isnan(Rs))) and (meth == "C35" or meth == "C40")):
+ elif (np.all(Rs) == None or np.all(np.isnan(Rs))):
Rs = np.ones(spd.shape)*150 # set to default for COARE3.5
- if ((np.all(np.isnan(zi))) and (meth == "C30" or meth == "C35" or
- meth == "C40")):
- zi = 600 # set to default for COARE3.5
- elif ((np.all(np.isnan(zi))) and (meth == "ERA5" or meth == "UA")):
- zi = 1000
- elif (np.all(np.isnan(zi))):
- zi = 800
+ if ((gust == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
+ gust = [1, 1.2, 600]
+ elif ((gust == None) and (meth == "UA" or meth == "ERA5")):
+ gust = [1, 1, 1000]
+ elif (gust == None):
+ gust = [1, 1.2, 800]
+ if (tol == None):
+ tol = ['flux', 0.01, 1, 1]
+ if ((jcool == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
+ or meth == "YT96")):
+ jcool = 0
+ elif ((jcool == None) and (meth == "UA" or meth == "LY04" or meth == "C30"
+ or meth == "C35" or meth == "C40"
+ or meth == "ERA5")):
+ jcool = 1
+ logging.info('method %s, inputs: h_in: %s | lat: %s | P: %s | Rl: %s |'
+ ' Rs: %s | gust: %s | jcool: %s', meth, h_in[:, 1],
+ np.nanmedian(lat), np.nanmedian(P), np.nanmedian(Rl),
+ np.nanmedian(Rs), gust, jcool)
####
th = np.where(T < 200, (np.copy(T)+CtoK) *
np.power(1000/P,287.1/1004.67),
@@ -135,22 +162,32 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
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
sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
- if qmeth:
+ if (hum == None):
+ RH = np.ones(spd.shape)*80
+ 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)
+ elif (hum[0] == 0):
+ RH = hum[1]
+ if (np.all(RH) < 1):
+ sys.exit("input relative humidity units should be \%")
+ 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)
+ elif (hum[0] == 1):
+ qair = hum[1]
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)
- logging.info('method %s and q method %s | qsea:%s, qair:%s', meth,
- qmeth, np.nanmedian(qsea), np.nanmedian(qair))
- else:
- qsea = qsat_sea(sst, P, "Buck2")/1000 # surface water q (g/kg)
- qair = qsat_air(T, P, RH, "Buck2")/1000 # q of air (g/kg)
- logging.info('method %s | qsea:%s, qair:%s', meth,
- np.nanmedian(qsea[~np.isnan(qsea)]),
- np.nanmedian(qair[~np.isnan(qair)]))
+ elif (hum[0] == 2):
+ Td = hum[1] # dew point temperature (K)
+ Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
+ 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)))
+ es = 611.21*np.exp(17.502*((T-273.16)/(T-32.19)))
+ RH = 100*esd/es
+ qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg)
+ qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
+ logging.info('method %s and q method %s | qsea:%s, qair:%s', meth, qmeth,
+ np.nanmedian(qsea), np.nanmedian(qair))
if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
print("qsea and qair cannot be nan")
- logging.info('method %s qsea and qair cannot be nan | sst:%s, Ta:%s,'
- 'P:%s, RH:%s', meth, np.nanmedian(sst), np.nanmedian(Ta),
- np.nanmedian(P), np.nanmedian(RH))
# first guesses
dt = Ta - sst
dq = qair - qsea
@@ -178,12 +215,12 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
Rnl = 0.97*(5.67e-8*np.power(sst-0.3*jcool+CtoK, 4)-Rl)
dter = np.ones(T.shape)*0.3
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
- if (gust == 1 and meth == "UA"):
+ if (gust[0] == 1 and meth == "UA"):
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)))
- elif (gust == 1):
+ elif (gust[0] == 1):
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
- elif (gust == 0):
+ elif (gust[0] == 0):
wind = np.copy(spd)
if (meth == "UA"):
@@ -193,55 +230,57 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
usr=kappa*wind/np.log(h_in[0]/zo)
Rb = g*h_in[0]*dtv/(tv*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))
+ (1-5*np.where(Rb < 0.19, Rb, 0.19)),
+ Rb*np.log(h_in[0]/zo))
monob = h_in[0]/zol
zo = 0.013*np.power(usr, 2)/g + 0.11*visc_air(Ta)/usr
zot = zo/np.exp(2.67*np.power(usr*zo/visc_air(Ta), 0.25)-2.57)
zoq = zot
logging.info('method %s | wind:%s, usr:%s, '
- 'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
- np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
- np.nanmedian(zot), np.nanmedian(zoq), np.nanmedian(Rb),
- np.nanmedian(monob))
+ 'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
+ np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
+ np.nanmedian(zot), np.nanmedian(zoq), np.nanmedian(Rb),
+ np.nanmedian(monob))
elif (meth == "ERA5"):
usr = np.sqrt(cd*np.power(wind, 2))
Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_in[0]) +
0.61*dq))/np.power(wind, 2))
zo = 0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g
zot = 0.40*visc_air(Ta)/usr
+ zoq = 0.62*visc_air(Ta)/usr
zol = (Rb*((np.log((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) /
- monob, meth)+psim_calc(zo/monob, meth)) /
- (np.log((h_in[0]+zo)/zot) -
- psit_calc((h_in[0]+zo)/monob, meth) +
- psit_calc(zot/monob, meth))))
+ monob, meth)+psim_calc(zo/monob, meth)) /
+ (np.log((h_in[0]+zo)/zot) -
+ psit_calc((h_in[0]+zo)/monob, meth) +
+ psit_calc(zot/monob, meth))))
monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
- 'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
- np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
- np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
+ 'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
+ np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
+ np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
elif (meth == "C30" or meth == "C35" or meth == "C40"):
usr = np.sqrt(cd*np.power(wind, 2))
a = 0.011*np.ones(T.shape)
a = np.where(wind > 10, 0.011+(wind-10)/(18-10)*(0.018-0.011),
- np.where(wind > 18, 0.018, a))
+ np.where(wind > 18, 0.018, a))
zo = a*np.power(usr, 2)/g+0.11*visc_air(T)/usr
rr = zo*usr/visc_air(T)
zoq = np.minimum(5e-5/np.power(rr, 0.6), 1.15e-4)
- zot = zoq
+ zot = np.copy(zoq)
Rb = g*h_in[0]*dtv/((T+CtoK)*np.power(wind, 2))
zol = (Rb*((np.log((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) /
- monob, meth)+psim_calc(zo/monob, meth)) /
- (np.log((h_in[0]+zo)/zot) -
- psit_calc((h_in[0]+zo)/monob, meth) +
- psit_calc(zot/monob, meth))))
+ monob, meth)+psim_calc(zo/monob, meth)) /
+ (np.log((h_in[0]+zo)/zot) -
+ psit_calc((h_in[0]+zo)/monob, meth) +
+ psit_calc(zot/monob, meth))))
monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
else:
- zo, zot = 0.0001*np.ones(spd.shape), 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)
usr = np.sqrt(cd*np.power(wind, 2))
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, monob:%s', meth,
@@ -249,18 +288,24 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
np.nanmedian(zot), np.nanmedian(monob))
tsr = (dt+dter*jcool)*kappa/(np.log(hin[1]/zot) -
psit_calc(h_in[1]/monob, meth))
- qsr = (dq+dqer*jcool)*kappa/(np.log(hin[2]/zot) -
+ qsr = (dq+dqer*jcool)*kappa/(np.log(hin[2]/zoq) -
psit_calc(h_in[2]/monob, meth))
- # tolerance for u,t,q,usr,tsr,qsr
- tol = np.array([0.01, 0.01, 5e-05, 0.005, 0.001, 5e-07])
it, ind = 0, np.where(spd > 0)
ii, itera = True, np.zeros(spd.shape)*np.nan
+ tau = 0.05*np.ones(spd.shape)
+ sensible = 5*np.ones(spd.shape)
+ latent = 65*np.ones(spd.shape)
while np.any(ii):
it += 1
if it > n:
break
- old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
- np.copy(usr), np.copy(tsr), np.copy(qsr)])
+ if (tol[0] == 'flux'):
+ old = np.array([np.copy(tau), np.copy(sensible), np.copy(latent)])
+ elif (tol[0] == 'ref'):
+ old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n)])
+ elif (tol[0] == 'all'):
+ old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
+ np.copy(tau), np.copy(sensible), np.copy(latent)])
cdn[ind] = cdn_calc(u10n[ind], Ta[ind], None, lat[ind], meth)
if (np.all(np.isnan(cdn))):
break
@@ -270,11 +315,15 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
cd[ind] = cd_calc(cdn[ind], h_in[0, ind], ref_ht, psim[ind])
ctn[ind], cqn[ind] = ctcqn_calc(h_in[1, ind]/monob[ind], cdn[ind],
u10n[ind], zo[ind], Ta[ind], meth)
+ zot[ind] = ref_ht/(np.exp(np.power(kappa, 2) /
+ (ctn[ind]*np.log(ref_ht/zo[ind]))))
+ zoq[ind] = ref_ht/(np.exp(np.power(kappa, 2) /
+ (cqn[ind]*np.log(ref_ht/zo[ind]))))
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
- h_in[1, ind], h_in[2, ind], ref_ht,
- psit[ind], psiq[ind])
+ h_in[1, ind], h_in[2, ind], ref_ht,
+ psit[ind], psiq[ind])
if (meth == "UA"):
usr[ind] = np.where(h_in[0, ind]/monob[ind] < -1.574,
kappa*wind[ind] /
@@ -351,10 +400,10 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
usr[ind] = (wind[ind]*kappa/(np.log(h_in[0, ind]/zo[ind]) -
psiu_26(h_in[0, ind]/monob[ind], meth)))
logging.info('method %s | dter = %s | Rnl = %s '
- '| usr = %s | tsr = %s | qsr = %s', meth,
- np.nanmedian(dter), np.nanmedian(Rnl),
- np.nanmedian(usr), np.nanmedian(tsr),
- np.nanmedian(qsr))
+ '| usr = %s | tsr = %s | qsr = %s', meth,
+ np.nanmedian(dter), np.nanmedian(Rnl),
+ np.nanmedian(usr), np.nanmedian(tsr),
+ np.nanmedian(qsr))
qsr[ind] = ((dq[ind]+dqer[ind]*jcool)*(kappa/(np.log(hin[2, ind] /
zoq[ind])-psit_26(hin[2, ind]/monob[ind]))))
tsr[ind] = ((dt[ind]+dter[ind]*jcool)*(kappa/(np.log(hin[1, ind] /
@@ -364,19 +413,24 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
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]
- dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
- rho[ind], Rl[ind],
- Rs[ind], Rnl[ind],
- cp[ind], lv[ind],
- np.copy(tkt[ind]),
- usr[ind], tsr[ind],
- qsr[ind], lat[ind])
+ if (jcool == 1):
+ dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
+ rho[ind], Rl[ind],
+ Rs[ind], Rnl[ind],
+ cp[ind], lv[ind],
+ np.copy(tkt[ind]),
+ usr[ind], tsr[ind],
+ qsr[ind], lat[ind])
+ else:
+ dter[ind] = np.zeros(sst[ind].shape)
+ dqer[ind] = np.zeros(sst[ind].shape)
+ tkt[ind] = np.zeros(sst[ind].shape)
Rnl[ind] = 0.97*(5.67e-8*np.power(SST[ind] -
- dter[ind]*jcool+CtoK, 4)-Rl[ind])
- t10n[ind] = (Ta[ind] +
- (tsr[ind]*(np.log(h_in[1, ind]/ref_ht)-psit[ind])/kappa))
- q10n[ind] = (qair[ind] +
- (qsr[ind]*(np.log(h_in[2, ind]/ref_ht)-psiq[ind])/kappa))
+ dter[ind]*jcool+CtoK, 4)-Rl[ind])
+ t10n[ind] = (Ta[ind] -
+ tsr[ind]/kappa*(np.log(h_in[1, ind]/ref_ht)-psit[ind]))
+ q10n[ind] = (qair[ind] -
+ qsr[ind]/kappa*(np.log(h_in[2, ind]/ref_ht)-psiq[ind]))
tv10n[ind] = t10n[ind]*(1+0.61*q10n[ind])
if (meth == "UA"):
tsrv[ind] = tsr[ind]*(1.+0.61*qair[ind])+0.61*th[ind]*qsr[ind]
@@ -389,11 +443,12 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
elif (meth == "ERA5"):
tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
Rb[ind] = ((g[ind]*h_in[0, ind]*((2*dt[ind])/(Ta[ind] +
- sst[ind]-g[ind]*h_in[0, ind])+0.61*dq[ind])) /
- np.power(wind[ind], 2))
+ sst[ind]-g[ind]*h_in[0, ind])+0.61*dq[ind])) /
+ np.power(wind[ind], 2))
zo[ind] = (0.11*visc_air(Ta[ind])/usr[ind]+0.018 *
np.power(usr[ind], 2)/g[ind])
zot[ind] = 0.40*visc_air(Ta[ind])/usr[ind]
+ zoq[ind] = 0.62*visc_air(Ta[ind])/usr[ind]
zol[ind] = (Rb[ind]*((np.log((h_in[0, ind]+zo[ind])/zo[ind]) -
psim_calc((h_in[0, ind]+zo[ind])/monob[ind], meth) +
psim_calc(zo[ind]/monob[ind], meth)) /
@@ -407,49 +462,52 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
psim[ind] = psim_calc(h_in[0, 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)
- if (gust == 1 and meth == "UA"):
+ if (gust[0] == 1 and meth == "UA"):
wind[ind] = np.where(dtv[ind] >= 0, np.where(spd[ind] > 0.1,
- spd[ind], 0.1),
- np.sqrt(np.power(np.copy(spd[ind]), 2) +
- np.power(get_gust(1, tv[ind], usr[ind],
- tsrv[ind], zi, lat[ind]), 2)))
- # Zeng et al. 1998 (20)
- elif (gust == 1 and (meth == "C30" or meth == "C35" or meth == "C40")):
+ spd[ind], 0.1),
+ np.sqrt(np.power(np.copy(spd[ind]), 2) +
+ np.power(get_gust(gust[1], tv[ind], usr[ind],
+ tsrv[ind], gust[2], lat[ind]), 2)))
+ # Zeng et al. 1998 (20)
+ elif (gust[0] == 1 and (meth == "C30" or meth == "C35" or
+ meth == "C40")):
wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
- np.power(get_gust(1.2, Ta[ind], usr[ind],
- tsrv[ind], zi, lat[ind]), 2))
- elif (gust == 1):
+ np.power(get_gust(gust[1], Ta[ind], usr[ind],
+ tsrv[ind], gust[2], lat[ind]), 2))
+ elif (gust[0] == 1):
wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
- np.power(get_gust(1, Ta[ind], usr[ind],
- tsrv[ind], zi, lat[ind]), 2))
- elif (gust == 0):
+ np.power(get_gust(gust[1], Ta[ind], usr[ind],
+ tsrv[ind], gust[2], lat[ind]), 2))
+ elif (gust[0] == 0):
wind[ind] = np.copy(spd[ind])
- if (meth == "UA"):
- u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(h_in[0, ind]/10) -
- psim[ind]))
- u10n[u10n < 0] = np.nan
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
- u10n[ind] = ((wind[ind]+usr[ind]/kappa*(np.log(10/h_in[0, ind]) -
- psiu_26(10/monob[ind], meth) +
- psiu_26(h_in[0, ind]/monob[ind], meth)) +
- psiu_26(10/monob[ind], meth)*usr[ind]/kappa) /
- (wind[ind]/spd[ind]))
- u10n[u10n < 0] = np.nan
- elif (meth == "ERA5"):
- u10n[ind] = (spd[ind]-(usr[ind]/kappa)*(np.log(h_in[0, ind] /
- ref_ht)-psim[ind]))
- u10n[u10n < 0] = np.nan
- else:
- u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(h_in[0, ind]/10) -
- psim[ind]))
- u10n[u10n < 0] = np.nan
+ u10n[ind] = wind[ind]-usr[ind]/kappa*(np.log(h_in[0, ind]/10) -
+ psim[ind])
+ u10n[u10n < 0] = np.nan
itera[ind] = np.ones(1)*it
- new = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
- np.copy(usr), np.copy(tsr), np.copy(qsr)])
- d = np.abs(new-old)
- ind = np.where((d[0, :] > tol[0])+(d[1, :] > tol[1]) +
- (d[2, :] > tol[2])+(d[3, :] > tol[3]) +
- (d[4, :] > tol[4])+(d[5, :] > tol[5]))
+ sensible = -rho*cp*usr*tsr
+ latent = -rho*lv*usr*qsr
+ if (gust[0] == 1):
+ tau = rho*np.power(usr, 2)*(spd/wind)
+ elif (gust[0] == 0):
+ tau = rho*np.power(usr, 2)
+ if (tol[0] == 'flux'):
+ new = np.array([np.copy(tau), np.copy(sensible), np.copy(latent)])
+ elif (tol[0] == 'ref'):
+ new = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n)])
+ elif (tol[0] == 'all'):
+ new = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
+ np.copy(tau), np.copy(sensible), np.copy(latent)])
+ d = np.fabs(new-old)
+ if (tol[0] == 'flux'):
+ ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
+ (d[2, :] > tol[3]))
+ elif (tol[0] == 'ref'):
+ ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
+ (d[2, :] > tol[3]))
+ elif (tol[0] == 'all'):
+ ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
+ (d[2, :] > tol[3])+(d[3, :] > tol[4]) +
+ (d[4, :] > tol[5])+(d[5, :] > tol[6]))
if (ind[0].size == 0):
ii = False
else:
@@ -460,22 +518,16 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
# calculate output parameters
rho = (0.34838*P)/(tv10n)
t10n = t10n-(273.16+tlapse*ref_ht)
- sensible = -rho*cp*usr*tsr
- latent = -rho*lv*usr*qsr
- if (gust == 1):
- tau = rho*np.power(usr, 2)*(spd/wind)
- elif (gust == 0):
- tau = rho*np.power(usr, 2)
- zo = ref_ht/np.exp(kappa/cdn**0.5)
- 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))))
- urefs = (spd-(usr/kappa)*(np.log(h_in[0]/h_out[0])-psim +
- psim_calc(h_out[0]/monob, meth)))
- trefs = (Ta-(tsr/kappa)*(np.log(h_in[1]/h_out[1])-psit +
- psit_calc(h_out[0]/monob, meth)))
- trefs = trefs-(273.16+tlapse*h_out[1])
- qrefs = (qair-(qsr/kappa)*(np.log(h_in[2]/h_out[2]) -
- psit+psit_calc(h_out[2]/monob, meth)))
+ # zo = ref_ht/np.exp(kappa/cdn**0.5)
+ # 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))))
+ uref = (spd-usr/kappa*(np.log(h_in[0]/h_out[0])-psim +
+ psim_calc(h_out[0]/monob, meth)))
+ tref = (Ta-tsr/kappa*(np.log(h_in[1]/h_out[1])-psit +
+ psit_calc(h_out[0]/monob, meth)))
+ tref = tref-(273.16+tlapse*h_out[1])
+ qref = (qair-qsr/kappa*(np.log(h_in[2]/h_out[2]) -
+ psit+psit_calc(h_out[2]/monob, meth)))
res = np.zeros((28, len(spd)))
res[0][:] = tau
res[1][:] = sensible
@@ -501,8 +553,13 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
res[21][:] = zo
res[22][:] = zot
res[23][:] = zoq
- res[24][:] = urefs
- res[25][:] = trefs
- res[26][:] = qrefs
+ res[24][:] = uref
+ res[25][:] = tref
+ res[26][:] = qref
res[27][:] = itera
+ if (out == 0):
+ res[:, ind] = np.nan
+ # set missing values where data have non acceptable values
+ res = np.where(spd < 0, np.nan, res)
return res, ind
+
diff --git a/flux_subs.py b/flux_subs.py
index 4a149d71139419a006d2511972e5aa51d300cd3d..a8054434c4855b7fe1d754b864131586a900f609 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -2,11 +2,39 @@ import numpy as np
from VaporPressure import VaporPressure
CtoK = 273.16 # 273.15
-""" Conversion factor for (^\circ\,C) to (^\\circ\\,K) """
+""" Conversion factor for $^\circ\,$C to K """
kappa = 0.4 # NOTE: 0.41
""" von Karman's constant """
# ---------------------------------------------------------------------
+def get_stabco(meth="S80"):
+ """ Gives the coefficients \\alpha, \\beta, \\gamma for stability functions
+
+ Parameters
+ ----------
+ meth : str
+
+ Returns
+ -------
+ coeffs : float
+ """
+ alpha, beta, gamma = 0, 0, 0
+ if (meth == "S80" or meth == "S88" or meth == "LY04" or
+ meth == "UA" or meth == "ERA5" or meth == "C30" or meth == "C35" or
+ meth == "C40"):
+ 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:
+ print("unknown method stabco: "+meth)
+ coeffs = np.zeros(3)
+ coeffs[0] = alpha
+ coeffs[1] = beta
+ coeffs[2] = gamma
+ return coeffs
+# ---------------------------------------------------------------------
def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
@@ -20,6 +48,8 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
air temperature (K)
Tp : float
wave period
+ lat : float
+ latitude
meth : str
Returns
@@ -63,6 +93,8 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
air temperature (K)
Tp : float
wave period
+ lat : float
+ latitude
meth : str
Returns
@@ -92,9 +124,9 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
elif (meth == "C35"):
a = 0.011*np.ones(Ta.shape)
-# a = np.where(u10n > 19, 0.0017*19-0.0050,
-# np.where((u10n > 7) & (u10n <= 18),
-# 0.0017*u10n-0.0050, a))
+ # a = np.where(u10n > 19, 0.0017*19-0.0050,
+ # np.where((u10n > 7) & (u10n <= 18),
+ # 0.0017*u10n-0.0050, a))
a = np.where(u10n > 19, 0.0017*19-0.0050, 0.0017*u10n-0.0050)
zo = 0.11*visc_air(Ta)/usr+a*np.power(usr, 2)/g
elif (meth == "C40"):
@@ -102,13 +134,36 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
a = np.where(u10n > 22, 0.0016*22-0.0035, 0.0016*u10n-0.0035)
zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr # surface roughness
elif (meth == "ERA5"):
- # eq. (3.26) p.40 over sea IFS Documentation cy46r1
+ # eq. (3.26) p.38 over sea IFS Documentation cy46r1
zo = 0.018*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
else:
print("unknown method for cdn_from_roughness "+meth)
cdnn = (kappa/np.log(10/zo))**2
cdn = np.where(np.abs(cdnn-cdn) < tol, cdnn, np.nan)
- return cdnn
+ return cdn
+# ---------------------------------------------------------------------
+
+
+def cd_calc(cdn, height, ref_ht, psim):
+ """ Calculates drag coefficient at reference height
+
+ Parameters
+ ----------
+ cdn : float
+ neutral drag coefficient
+ height : float
+ original sensor height (m)
+ ref_ht : float
+ reference height (m)
+ psim : float
+ momentum stability function
+
+ Returns
+ -------
+ cd : float
+ """
+ cd = (cdn/np.power(1+(np.sqrt(cdn)*(np.log(height/ref_ht)-psim))/kappa, 2))
+ return cd
# ---------------------------------------------------------------------
@@ -120,7 +175,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
zol : float
height over MO length
cdn : float
- neatral drag coefficient
+ neutral drag coefficient
u10n : float
neutral 10m wind speed (m/s)
zo : float
@@ -187,7 +242,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
elif (meth == "ERA5"):
- # eq. (3.26) p.40 over sea IFS Documentation cy46r1
+ # eq. (3.26) p.38 over sea IFS Documentation cy46r1
usr = np.sqrt(cdn*np.power(u10n, 2))
zot = 0.40*visc_air(Ta)/usr
zoq = 0.62*visc_air(Ta)/usr
@@ -199,30 +254,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
# ---------------------------------------------------------------------
-def cd_calc(cdn, height, ref_ht, psim):
- """ Calculates drag coefficient at reference height
-
- Parameters
- ----------
- cdn : float
- neutral drag coefficient
- height : float
- original sensor height (m)
- ref_ht : float
- reference height (m)
- psim : float
- momentum stability function
-
- Returns
- -------
- cd : float
- """
- cd = (cdn/np.power(1+(np.sqrt(cdn)*(np.log(height/ref_ht)-psim))/kappa, 2))
- return cd
-# ---------------------------------------------------------------------
-
-
-def ctcq_calc(cdn, cd, ctn, cqn, h_t, h_q, ref_ht, psit, psiq):
+def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq):
""" Calculates heat and moisture exchange coefficients at reference height
Parameters
@@ -249,12 +281,14 @@ def ctcq_calc(cdn, cd, ctn, cqn, h_t, h_q, ref_ht, psit, psiq):
Returns
-------
ct : float
+ heat exchange coefficient
cq : float
+ moisture exchange coefficient
"""
ct = (ctn*np.sqrt(cd/cdn) /
- (1+ctn*((np.log(h_t/ref_ht)-psit)/(kappa*np.sqrt(cdn)))))
+ (1+ctn*((np.log(ht/ref_ht)-psit)/(kappa*np.sqrt(cdn)))))
cq = (cqn*np.sqrt(cd/cdn) /
- (1+cqn*((np.log(h_q/ref_ht)-psiq)/(kappa*np.sqrt(cdn)))))
+ (1+cqn*((np.log(hq/ref_ht)-psiq)/(kappa*np.sqrt(cdn)))))
return ct, cq
# ---------------------------------------------------------------------
@@ -272,16 +306,13 @@ def psim_calc(zol, meth="S80"):
-------
psim : float
"""
- coeffs = get_stabco(meth)
- alpha, beta, gamma = coeffs[0], coeffs[1], coeffs[2]
if (meth == "ERA5"):
- psim = np.where(zol < 0, psim_conv(zol, alpha, beta, gamma),
- psim_stab_era5(zol, alpha, beta, gamma))
+ psim = psim_era5(zol)
elif (meth == "C30" or meth == "C35" or meth == "C40"):
psim = psiu_26(zol, meth)
else:
- psim = np.where(zol < 0, psim_conv(zol, alpha, beta, gamma),
- psim_stab(zol, alpha, beta, gamma))
+ psim = np.where(zol < 0, psim_conv(zol, meth),
+ psim_stab(zol, meth))
return psim
# ---------------------------------------------------------------------
@@ -294,55 +325,25 @@ def psit_calc(zol, meth="S80"):
zol : float
height over MO length
meth : str
+ parameterisation method
Returns
-------
psit : float
"""
- coeffs = get_stabco(meth)
- alpha, beta, gamma = coeffs[0], coeffs[1], coeffs[2]
if (meth == "ERA5"):
- psit = np.where(zol < 0, psi_conv(zol, alpha, beta, gamma),
- psi_stab_era5(zol, alpha, beta, gamma))
+ psit = np.where(zol < 0, psi_conv(zol, meth),
+ psi_era5(zol))
elif (meth == "C30" or meth == "C35" or meth == "C40"):
psit = psit_26(zol)
else:
- psit = np.where(zol < 0, psi_conv(zol, alpha, beta, gamma),
- psi_stab(zol, alpha, beta, gamma))
+ psit = np.where(zol < 0, psi_conv(zol, meth),
+ psi_stab(zol, meth))
return psit
# ---------------------------------------------------------------------
-def get_stabco(meth="S80"):
- """ Gives the coefficients \\alpha, \\beta, \\gamma for stability functions
-
- Parameters
- ----------
- meth : str
-
- Returns
- -------
- coeffs : float
- """
- if (meth == "S80" or meth == "S88" or meth == "LY04" or
- meth == "UA" or meth == "ERA5" or meth == "C30" or meth == "C35" or
- meth == "C40"):
- 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:
- print("unknown method stabco: "+meth)
- coeffs = np.zeros(3)
- coeffs[0] = alpha
- coeffs[1] = beta
- coeffs[2] = gamma
- return coeffs
-# ---------------------------------------------------------------------
-
-
-def psi_stab_era5(zol, alpha, beta, gamma):
+def psi_era5(zol):
""" Calculates heat stability function for stable conditions
for method ERA5
@@ -350,14 +351,12 @@ def psi_stab_era5(zol, alpha, beta, gamma):
----------
zol : float
height over MO length
- alpha, beta, gamma : float
- constants given by get_stabco
Returns
-------
psit : float
"""
- # eq (3.22) p. 39 IFS Documentation cy46r1
+ # 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
@@ -377,94 +376,103 @@ def psit_26(zol):
psi : float
"""
b, d = 2/3, 0.35
- dzol = np.where(d*zol > 50, 50, d*zol) # stable
- psi = -((1+b*zol)**1.5+b*(zol-14.28)*np.exp(-dzol)+8.525)
- psik = 2*np.log((1+np.sqrt(1-15*zol))/2)
- psic = (1.5*np.log((1+np.power(1-34.15*zol, 1/3) +
- np.power(1-34.15*zol, 2/3))/3)-np.sqrt(3) *
- np.arctan(1+2*np.power(1-34.15*zol, 1/3))/np.sqrt(3) +
- 4*np.arctan(1)/np.sqrt(3))
+ dzol = np.where(d*zol > 50, 50, d*zol)
+ psi = np.where(zol > 0,-(np.power(1+b*zol, 1.5)+b*(zol-14.28) *
+ np.exp(-dzol)+8.525), np.nan)
+ psik = np.where(zol < 0, 2*np.log((1+np.sqrt(1-15*zol))/2), np.nan)
+ psic = np.where(zol < 0, 1.5*np.log((1+np.power(1-34.15*zol, 1/3) +
+ np.power(1-34.15*zol, 2/3))/3)-np.sqrt(3) *
+ np.arctan(1+2*np.power(1-34.15*zol, 1/3))/np.sqrt(3) +
+ 4*np.arctan(1)/np.sqrt(3), np.nan)
f = np.power(zol, 2)/(1+np.power(zol, 2))
psi = np.where(zol < 0, (1-f)*psik+f*psic, psi)
return psi
# ---------------------------------------------------------------------
-def psi_conv(zol, alpha, beta, gamma):
+def psi_conv(zol, meth):
""" Calculates heat stability function for unstable conditions
Parameters
----------
zol : float
height over MO length
- alpha, beta, gamma : float
- constants given by get_stabco
+ 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, alpha, beta, gamma):
+def psi_stab(zol, meth):
""" Calculates heat stability function for stable conditions
Parameters
----------
zol : float
height over MO length
- alpha, beta, gamma : float
- constants given by get_stabco
+ meth : str
+ parameterisation method
Returns
-------
psit : float
"""
+ coeffs = get_stabco(meth)
+ gamma = coeffs[2]
psit = -gamma*zol
return psit
# ---------------------------------------------------------------------
-def psim_stab_era5(zol, alpha, beta, gamma):
- """ Calculates momentum stability function for stable conditions
- for method ERA5
+def psim_era5(zol):
+ """ Calculates momentum stability function for method ERA5
Parameters
----------
zol : float
height over MO length
- alpha, beta, gamma : float
- constants given by get_stabco
Returns
-------
psim : float
"""
- # eq (3.22) p. 39 IFS Documentation cy46r1
+ # eq (3.20, 3.22) p. 37 IFS Documentation cy46r1
+ coeffs = get_stabco("ERA5")
+ alpha, beta = coeffs[0], coeffs[1]
+ xtmp = np.power(1-alpha*zol, beta)
a, b, c, d = 1, 2/3, 5, 0.35
- psim = -b*(zol-c/d)*np.exp(-d*zol)-a*zol-(b*c)/d
+ 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 psim_conv(zol, alpha, beta, gamma):
+def psim_conv(zol, meth):
""" Calculates momentum stability function for unstable conditions
Parameters
----------
zol : float
height over MO length
- alpha, beta, gamma : float
- constants given by get_stabco
+ 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)
@@ -472,20 +480,22 @@ def psim_conv(zol, alpha, beta, gamma):
# ---------------------------------------------------------------------
-def psim_stab(zol, alpha, beta, gamma):
+def psim_stab(zol, meth):
""" Calculates momentum stability function for stable conditions
Parameters
----------
zol : float
height over MO length
- alpha, beta, gamma : float
- constants given by get_stabco
+ meth : str
+ parameterisation method
Returns
-------
psim : float
"""
+ coeffs = get_stabco(meth)
+ gamma = coeffs[2]
psim = -gamma*zol
return psim
# ---------------------------------------------------------------------
@@ -505,28 +515,31 @@ def psiu_26(zol, meth):
"""
if (meth == "C30"):
dzol = np.where(0.35*zol > 50, 50, 0.35*zol) # stable
- psi = -((1+zol)+0.6667*(zol-14.28)*np.exp(-dzol)+8.525)
- k = np.where(zol < 0) # unstable
- x = np.power(1-15*zol[k], 0.25)
- psik = (2*np.log((1+x)/2)+np.log((1+np.power(x, 2))/2)-2*np.arctan(x) +
- 2*np.arctan(1))
- x = np.power(1-10.15*zol[k], 0.3333)
- 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
+ psi = np.where(zol > 0, -((1+zol)+0.6667*(zol-14.28)*np.exp(-dzol) +
+ 8.525), np.nan)
+ x = np.where(zol < 0, np.power(1-15*zol, 0.25), np.nan)
+ psik = np.where(zol < 0, 2*np.log((1+x)/2)+np.log((1+np.power(x, 2)) /
+ 2)-2*np.arctan(x)+2*np.arctan(1), np.nan)
+ x = np.where(zol < 0, np.power(1-10.15*zol, 0.3333), np.nan)
+ psic = np.where(zol < 0, 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), np.nan)
+ f = np.power(zol, 2)/(1+np.power(zol, 2))
+ psi = np.where(zol < 0, (1-f)*psik+f*psic, psi)
elif (meth == "C35" or meth == "C40"):
dzol = np.where(0.35*zol > 50, 50, 0.35*zol) # stable
a, b, c, d = 0.7, 3/4, 5, 0.35
- psi = -(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], 0.25)
- psik = 2*np.log((1+x)/2)+np.log((1+x**2)/2)-2*np.arctan(x)+2*np.arctan(1)
- x = np.power(1-10.15*zol[k], 0.3333)
- 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
+ psi = np.where(zol > 0, -(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d),
+ np.nan)
+ x = np.where(zol < 0, np.power(1-15*zol, 0.25), np.nan)
+ psik = np.where(zol < 0, 2*np.log((1+x)/2)+np.log((1+x**2)/2) -
+ 2*np.arctan(x)+2*np.arctan(1), np.nan)
+ x = np.where(zol < 0, np.power(1-10.15*zol, 0.3333), np.nan)
+ psic = np.where(zol < 0, 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), np.nan)
+ f = np.power(zol, 2)/(1+np.power(zol, 2))
+ psi = np.where(zol < 0, (1-f)*psik+f*psic, psi)
return psi
# ---------------------------------------------------------------------
@@ -638,6 +651,8 @@ def get_heights(h, dim_len):
----------
h : float
input heights (m)
+ dim_len : int
+ length dimension
Returns
-------
@@ -646,17 +661,17 @@ def get_heights(h, dim_len):
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 h.ndim == 1):
+ 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 h.ndim == 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 h.ndim == 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 h.ndim == 2):
+ 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 h.ndim == 2):
+ elif (len(h) == 3 and np.ndim(h) == 2):
hh = np.zeros((3, h.shape[1]))
hh = np.copy(h)
return hh