Skip to content

GitLab

Commit edcb5c52 authored by sbiri's avatar sbiri
Browse files

adjusted get_height subroutine to read in multiple (3xn) height input. Added... 

adjusted get_height subroutine to read in multiple (3xn) height input. Added psiq as output. Included all switches jcool, gust, qmeth.
parent 69316791
Hide whitespace changes
Inline Side-by-side
Showing with 215 additions and 237 deletions
AirSeaFluxCode.py
View file @ edcb5c52
......@@ -3,18 +3,23 @@ import sys
import logging
from flux_subs import (kappa, CtoK, get_heights, cdn_calc, cd_calc, get_skin,
psim_calc, psit_calc, ctcq_calc, ctcqn_calc, get_gust,
gc, q_calc, qsea_calc, qsat26sea, qsat26air,
visc_air, psit_26, psiu_26)
gc, qsat_sea, qsat_air, visc_air, psit_26, psiu_26)
def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
Rl=None, Rs=None, jcool=1, meth="S88", n=10):
def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
gust, meth, qmeth, n):
""" 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)
......@@ -29,8 +34,7 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
P : float
air pressure
hin : float
sensor heights in m (array of 1->3 values: 1 -> u=t=q; /
2 -> u,t=q; 3 -> u,t,q ) default 10m
sensor heights in m (array 3x1 or 3xn) default 10m ([10, 10, 10])
hout : float
output height, default is 10m
zi : int
......@@ -41,6 +45,8 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
downward shortwave radiation (W/m^2)
jcool : int
0 if sst is true ocean skin temperature, else 1
gust : int
1 to include the effect of gustiness, else 0
n : int
number of iterations
......@@ -62,18 +68,19 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
13. star humidity (qsr)
14. star velocity (usr)
15. momentum stability function (psim)
16. heat stability funciton (psit)
17. 10m neutral velocity (u10n)
18. 10m neutral temperature (t10n)
19. 10m neutral virtual temperature (tv10n)
20. 10m neutral specific humidity (q10n)
21. surface roughness length (zo)
22. heat roughness length (zot)
23. moisture roughness length (zoq)
24. velocity at reference height (urefs)
25. temperature at reference height (trefs)
26. specific humidity at reference height (qrefs)
27. number of iterations until convergence
16. heat stability function (psit)
17. moisture stability function (psiq)
18. 10m neutral velocity (u10n)
19. 10m neutral temperature (t10n)
20. 10m neutral virtual temperature (tv10n)
21. 10m neutral specific humidity (q10n)
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)
28. number of iterations until convergence
ind : int
the indices in the matrix for the points that did not converge
after the maximum number of iterations
......@@ -81,15 +88,12 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
"""
logging.basicConfig(filename='flux_calc.log',
format='%(asctime)s %(message)s',level=logging.INFO)
ref_ht, tlapse = 10, 0.0098 # reference height, lapse rate
hh_in = get_heights(hin) # heights of input measurements/fields
hh_out = get_heights(hout) # desired height of output variables
if np.all(np.isnan(lat)): # set latitude to 45deg if empty
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
ctn, ct, cqn, cq = (np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan,
np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan)
cdn = np.zeros(spd.shape)*np.nan
# 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")
......@@ -121,34 +125,28 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
elif ((np.all(np.isnan(zi))) and (meth == "ERA5" or meth == "UA")):
zi = 1000
####
th = np.where(np.nanmax(T) < 200, (np.copy(T)+CtoK) *
th = np.where(T < 200, (np.copy(T)+CtoK) *
np.power(1000/P,287.1/1004.67),
np.copy(T)*np.power(1000/P,287.1/1004.67)) # potential T
Ta = np.where(np.nanmax(T) < 200, np.copy(T)+CtoK+tlapse*hh_in[1],
np.copy(T)+tlapse*hh_in[1]) # convert to Kelvin if needed
sst = np.where(np.nanmax(SST) < 200, np.copy(SST)+CtoK, np.copy(SST))
if (meth == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
meth == "ERA5"):
qsea = qsat26sea(sst, P)/1000 # surface water specific humidity (g/kg)
Q, _ = qsat26air(T, P, RH) # specific humidity of air (g/kg)
qair = Q/1000
del Q
logging.info('method %s | qsea:%s, qair:%s', meth,
np.ma.median(qsea[~np.isnan(qsea)]),
np.ma.median(qair[~np.isnan(qair)]))
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:
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 = qsea_calc(sst, P)
qair = q_calc(Ta, RH, P)
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.ma.median(qsea[~np.isnan(qsea)]),
np.ma.median(qair[~np.isnan(qair)]))
np.nanmedian(qsea[~np.isnan(qsea)]),
np.nanmedian(qair[~np.isnan(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.ma.median(sst[~np.isnan(sst)]),
np.ma.median(Ta[~np.isnan(Ta)]),
np.ma.median(P[~np.isnan(P)]),
np.ma.median(RH[~np.isnan(RH)]))
'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
......@@ -159,70 +157,66 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
dtv=dt*(1.+0.61*qair)+0.61*th*dq
# ------------
rho = P*100/(287.1*tv10n)
# rho=P*100/(287.1*sst*(1+0.61*qsea)) # Zeng et al. 1998
lv = (2.501-0.00237*SST)*1e6
cp = 1004.67*(1 + 0.00084*qsea)
# cp = 1004.67*np.ones(Ta.shape) # Zeng et al. 1998, C3.0, C3.5
u10n = np.copy(spd)
monob = -100*np.ones(u10n.shape) # Monin-Obukhov length
monob = -100*np.ones(spd.shape) # Monin-Obukhov length
cdn = cdn_calc(u10n, Ta, None, lat, meth)
psim, psit, psiq = (np.zeros(u10n.shape), np.zeros(u10n.shape),
np.zeros(u10n.shape))
cd = cd_calc(cdn, hh_in[0], ref_ht, psim)
tsr, tsrv = np.zeros(u10n.shape), np.zeros(u10n.shape)
qsr = np.zeros(u10n.shape)
# if jcool == 1:
ctn, ct, cqn, cq = (np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan,
np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan)
psim, psit, psiq = (np.zeros(spd.shape), np.zeros(spd.shape),
np.zeros(spd.shape))
cd = cd_calc(cdn, h_in[0], ref_ht, psim)
tsr, tsrv = np.zeros(spd.shape), np.zeros(spd.shape)
qsr = np.zeros(spd.shape)
# jcool parameters
tkt = 0.001*np.ones(T.shape)
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 (meth == "UA"):
if (gust == 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):
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
elif (gust == 0):
wind = np.copy(spd)
if (meth == "UA"):
usr = 0.06
for i in range(5):
zo = 0.013*np.power(usr,2)/g+0.11*visc_air(Ta)/usr
usr=kappa*wind/np.log(hh_in[0]/zo)
Rb = g*hh_in[0]*dtv/(tv*wind**2)
zol = np.where(Rb >= 0, Rb*np.log(hh_in[0]/zo) /
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(hh_in[0]/zo))
monob = hh_in[0]/zol
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.ma.median(wind[~np.isnan(wind)]),
np.ma.median(usr[~np.isnan(usr)]),
np.ma.median(zo[~np.isnan(zo)]),
np.ma.median(zot[~np.isnan(zot)]),
np.ma.median(zoq[~np.isnan(zoq)]),
np.ma.median(Rb[~np.isnan(Rb)]),
np.ma.median(monob[~np.isnan(monob)]))
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(zoq), np.nanmedian(Rb),
np.nanmedian(monob))
elif (meth == "ERA5"):
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
usr = np.sqrt(cd*np.power(wind, 2))
Rb = ((g*hh_in[0]*((2*dt)/(Ta+sst-g*hh_in[0]) +
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
zol = (Rb*((np.log((hh_in[0]+zo)/zo)-psim_calc((hh_in[0]+zo) /
zol = (Rb*((np.log((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) /
monob, meth)+psim_calc(zo/monob, meth)) /
(np.log((hh_in[0]+zo)/zot) -
psit_calc((hh_in[0]+zo)/monob, meth) +
(np.log((h_in[0]+zo)/zot) -
psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
monob = hh_in[0]/zol
monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
np.ma.median(wind[~np.isnan(wind)]),
np.ma.median(usr[~np.isnan(usr)]),
np.ma.median(zo[~np.isnan(zo)]),
np.ma.median(zot[~np.isnan(zot)]),
np.ma.median(Rb[~np.isnan(Rb)]),
np.ma.median(monob[~np.isnan(monob)]))
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"):
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
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),
......@@ -231,38 +225,30 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
rr = zo*usr/visc_air(T)
zoq = np.minimum(5e-5/np.power(rr, 0.6), 1.15e-4)
zot = zoq
Rb = g*hh_in[0]*dtv/((T+CtoK)*np.power(wind, 2))
zol = (Rb*((np.log((hh_in[0]+zo)/zo)-psim_calc((hh_in[0]+zo) /
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((hh_in[0]+zo)/zot) -
psit_calc((hh_in[0]+zo)/monob, meth) +
(np.log((h_in[0]+zo)/zot) -
psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
monob = hh_in[0]/zol
monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
np.ma.median(wind[~np.isnan(wind)]),
np.ma.median(usr[~np.isnan(usr)]),
np.ma.median(zo[~np.isnan(zo)]),
np.ma.median(zot[~np.isnan(zot)]),
np.ma.median(Rb[~np.isnan(Rb)]),
np.ma.median(monob[~np.isnan(monob)]))
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
else:
wind = np.copy(spd)
zo, zot = 0.0001*np.ones(u10n.shape), 0.0001*np.ones(u10n.shape)
zo, zot = 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, Rb:%s, monob:%s', meth,
np.ma.median(wind[~np.isnan(wind)]),
np.ma.median(usr[~np.isnan(usr)]),
np.ma.median(zo[~np.isnan(zo)]),
np.ma.median(zot[~np.isnan(zot)]),
np.ma.median(monob[~np.isnan(monob)]))
'zo:%s, zot:%s, monob:%s', meth,
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(monob))
tsr = (dt+dter*jcool)*kappa/(np.log(hin[1]/zot) -
psit_calc(hh_in[1]/monob, meth))
psit_calc(h_in[1]/monob, meth))
qsr = (dq+dqer*jcool)*kappa/(np.log(hin[2]/zot) -
psit_calc(hh_in[2]/monob, meth))
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]) #[1e-06, 0.01, 5e-07, 1e-06, 0.001, 5e-07]
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
while np.any(ii):
......@@ -276,180 +262,183 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
break
logging.info('break %s at iteration %s cdn<0', meth, it)
zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cdn[ind]))
psim[ind] = psim_calc(hh_in[0]/monob[ind], meth)
cd[ind] = cd_calc(cdn[ind], hh_in[0], ref_ht, psim[ind])
ctn[ind], cqn[ind] = ctcqn_calc(hh_in[1]/monob[ind], cdn[ind],
psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
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)
psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
psiq[ind] = psit_calc(hh_in[2]/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)
ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
hh_in[1], hh_in[2], ref_ht,
h_in[1, ind], h_in[2, ind], ref_ht,
psit[ind], psiq[ind])
if (meth == "UA"):
usr[ind] = np.where(hh_in[0]/monob[ind] < -1.574, kappa*wind[ind] /
usr[ind] = np.where(h_in[0, ind]/monob[ind] < -1.574,
kappa*wind[ind] /
(np.log(-1.574*monob[ind]/zo[ind]) -
psim_calc(-1.574, meth) +
psim_calc(zo[ind]/monob[ind], meth) +
1.14*(np.power(-hh_in[0]/monob[ind], 1/3) -
1.14*(np.power(-h_in[0, ind]/monob[ind], 1/3) -
np.power(1.574, 1/3))),
np.where((hh_in[0]/monob[ind] > -1.574) &
(hh_in[0]/monob[ind] < 0),
kappa*wind[ind]/(np.log(hh_in[0]/zo[ind]) -
psim_calc(hh_in[0]/monob[ind], meth) +
psim_calc(zo[ind]/monob[ind], meth)),
np.where((hh_in[0]/monob[ind] > 0) &
(hh_in[0]/monob[ind] < 1),
kappa*wind[ind]/(np.log(hh_in[0]/zo[ind]) +
5*hh_in[0]/monob[ind]-5*zo[ind]/monob[ind]),
kappa*wind[ind]/(np.log(monob[ind]/zo[ind]) +
np.where((h_in[0, ind]/monob[ind] > -1.574) &
(h_in[0, ind]/monob[ind] < 0),
kappa*wind[ind]/(np.log(h_in[0, ind]/zo[ind]) -
psim_calc(h_in[0, ind]/monob[ind], meth) +
psim_calc(zo[ind, ind]/monob[ind], meth)),
np.where((h_in[0, ind]/monob[ind] > 0) &
(h_in[0, ind]/monob[ind] < 1),
kappa*wind[ind]/(np.log(h_in[0, ind]/zo[ind]) +
5*h_in[0, ind]/monob[ind]-5*zo[ind] /
monob[ind]), kappa*wind[ind] /
(np.log(monob[ind]/zo[ind]) +
5-5*zo[ind]/monob[ind] +
5*np.log(hh_in[0]/monob[ind]) +
hh_in[0]/monob[ind]-1))))
5*np.log(h_in[0, ind]/monob[ind]) +
h_in[0, ind]/monob[ind]-1))))
# Zeng et al. 1998 (7-10)
tsr[ind] = np.where(hh_in[1]/monob[ind] < -0.465, kappa*(dt[ind] +
tsr[ind] = np.where(h_in[1, ind]/monob[ind] < -0.465,
kappa*(dt[ind] +
dter[ind]*jcool) /
(np.log((-0.465*monob[ind])/zot[ind]) -
psit_calc(-0.465, meth)+0.8 *
(np.power(0.465, -1/3) -
np.power(-hh_in[1]/monob[ind], -1/3))),
np.where((hh_in[1]/monob[ind] > -0.465) &
(hh_in[1]/monob[ind] < 0),
np.power(-h_in[1, ind]/monob[ind], -1/3))),
np.where((h_in[1, ind]/monob[ind] > -0.465) &
(h_in[1, ind]/monob[ind] < 0),
kappa*(dt[ind]+dter[ind]*jcool) /
(np.log(hh_in[1]/zot[ind]) -
psit_calc(hh_in[1]/monob[ind], meth) +
(np.log(h_in[1, ind]/zot[ind]) -
psit_calc(h_in[1, ind]/monob[ind], meth) +
psit_calc(zot[ind]/monob[ind], meth)),
np.where((hh_in[1]/monob[ind] > 0) &
(hh_in[1]/monob[ind] < 1),
np.where((h_in[1, ind]/monob[ind] > 0) &
(h_in[1, ind]/monob[ind] < 1),
kappa*(dt[ind]+dter[ind]*jcool) /
(np.log(hh_in[1]/zot[ind]) +
5*hh_in[1]/monob[ind]-5*zot[ind]/monob[ind]),
(np.log(h_in[1, ind]/zot[ind]) +
5*h_in[1, ind]/monob[ind]-5*zot[ind] /
monob[ind]),
kappa*(dt[ind]+dter[ind]*jcool) /
(np.log(monob[ind]/zot[ind])+5 -
5*zot[ind]/monob[ind] +
5*np.log(hh_in[1]/monob[ind]) +
hh_in[1]/monob[ind]-1))))
5*np.log(h_in[1, ind]/monob[ind]) +
h_in[1, ind]/monob[ind]-1))))
# Zeng et al. 1998 (11-14)
qsr[ind] = np.where(hh_in[2]/monob[ind] < -0.465, kappa*(dq[ind] +
qsr[ind] = np.where(h_in[2, ind]/monob[ind] < -0.465,
kappa*(dq[ind] +
dqer[ind]*jcool) /
(np.log((-0.465*monob[ind])/zoq[ind]) -
psit_calc(-0.465, meth) +
psit_calc(zoq[ind]/monob[ind], meth) +
0.8*(np.power(0.465, -1/3) -
np.power(-hh_in[2]/monob[ind], -1/3))),
np.where((hh_in[2]/monob[ind] > -0.465) &
(hh_in[2]/monob[ind] < 0),
np.power(-h_in[2, ind]/monob[ind], -1/3))),
np.where((h_in[2, ind]/monob[ind] > -0.465) &
(h_in[2, ind]/monob[ind] < 0),
kappa*(dq[ind]+dqer[ind]*jcool) /
(np.log(hh_in[1]/zot[ind]) -
psit_calc(hh_in[2]/monob[ind], meth) +
(np.log(h_in[1, ind]/zot[ind]) -
psit_calc(h_in[2, ind]/monob[ind], meth) +
psit_calc(zoq[ind]/monob[ind], meth)),
np.where((hh_in[2]/monob[ind] > 0) &
(hh_in[2]/monob[ind]<1), kappa*(dq[ind] +
np.where((h_in[2, ind]/monob[ind] > 0) &
(h_in[2, ind]/monob[ind]<1), kappa*(dq[ind] +
dqer[ind]*jcool) /
(np.log(hh_in[1]/zoq[ind]) +
5*hh_in[2]/monob[ind]-5*zoq[ind]/monob[ind]),
(np.log(h_in[1, ind]/zoq[ind]) +
5*h_in[2, ind]/monob[ind]-5*zoq[ind]/monob[ind]),
kappa*(dq[ind]+dqer[ind]*jcool) /
(np.log(monob[ind]/zoq[ind])+5 -
5*zoq[ind]/monob[ind] +
5*np.log(hh_in[2]/monob[ind]) +
hh_in[2]/monob[ind]-1))))
5*np.log(h_in[2, ind]/monob[ind]) +
h_in[2, ind]/monob[ind]-1))))
elif (meth == "C30" or meth == "C35" or meth == "C40"):
usr[ind] = (wind[ind]*kappa/(np.log(hh_in[0]/zo[ind]) -
psiu_26(hh_in[0]/monob[ind], meth)))
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.ma.median(dter[~np.isnan(dter)]),
np.ma.median(Rnl[~np.isnan(Rnl)]),
np.ma.median(usr[~np.isnan(usr)]),
np.ma.median(tsr[~np.isnan(tsr)]),
np.ma.median(qsr[~np.isnan(qsr)]))
qsr[ind] = ((dq[ind]+dqer[ind]*jcool)*(kappa /
(np.log(hin[2]/zoq[ind])-psit_26(hin[2]/monob[ind]))))
tsr[ind] = ((dt[ind]+dter[ind]*jcool)*(kappa /
(np.log(hin[1]/zot[ind])-psit_26(hin[1]/monob[ind]))))
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] /
zot[ind])-psit_26(hin[1, ind]/monob[ind]))))
else:
usr[ind] = np.sqrt(cd[ind]*np.power(wind[ind], 2))
# = (wind[ind]*kappa/(np.log(hh_in[0]/zo[ind]) -
# psim_calc(hh_in[0]/monob[ind], meth)))
usr[ind] = (wind[ind]*kappa/(np.log(h_in[0, ind]/zo[ind]) -
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]
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(spd[ind].shape)
dqer[ind] = np.zeros(spd[ind].shape)
tkt[ind] = np.zeros(spd[ind].shape)
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])
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(hh_in[1]/ref_ht)-psit[ind])/kappa))
q10n[ind] = (qair[ind] -
(qsr[ind]*(np.log(hh_in[2]/ref_ht)-psiq[ind])/kappa))
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))
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]
monob[ind] = (tv[ind]*np.power(usr[ind], 2))/(kappa*g[ind]*tsrv[ind])
elif (meth == "C30" or meth == "C35" or meth == "C40"):
tsrv[ind] = tsr[ind]+0.61*(T[ind]+CtoK)*qsr[ind]
zol[ind] = (kappa*g[ind]*hh_in[0]/(T[ind]+CtoK)*(tsr[ind] +
zol[ind] = (kappa*g[ind]*h_in[0, ind]/(T[ind]+CtoK)*(tsr[ind] +
0.61*(T[ind]+CtoK)*qsr[ind])/np.power(usr[ind], 2))
monob[ind] = hh_in[0]/zol[ind]
monob[ind] = h_in[0, ind]/zol[ind]
elif (meth == "ERA5"):
tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
Rb[ind] = ((g[ind]*hh_in[0]*((2*dt[ind])/(Ta[ind]+sst[ind]-g[ind] *
hh_in[0])+0.61*dq[ind]))/np.power(wind[ind], 2))
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))
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]
zol[ind] = (Rb[ind]*((np.log((hh_in[0]+zo[ind])/zo[ind]) -
psim_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
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)) /
(np.log((hh_in[0]+zo[ind])/zot[ind]) -
psit_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
(np.log((h_in[0, ind]+zo[ind])/zot[ind]) -
psit_calc((h_in[0, ind]+zo[ind])/monob[ind], meth) +
psit_calc(zot[ind]/monob[ind], meth))))
monob[ind] = hh_in[0]/zol[ind]
monob[ind] = h_in[0, ind]/zol[ind]
else:
tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
monob[ind] = (tv10n[ind]*usr[ind]**2)/(g[ind]*kappa*tsrv[ind])
psim[ind] = psim_calc(hh_in[0]/monob[ind], meth)
psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
psiq[ind] = psit_calc(hh_in[2]/monob[ind], meth)
if (meth == "UA"):
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"):
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)
u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
psim[ind]))
u10n[u10n < 0] = np.nan
elif (meth == "C30" or meth == "C35" or meth == "C40"):
elif (gust == 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))
u10n[ind] = ((wind[ind] + usr[ind]/kappa*(np.log(10/hh_in[0])-
elif (gust == 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):
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(hh_in[0]/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
# u10n[u10n < 0] = np.nan
elif (meth == "ERA5"):
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))
u10n[ind] = (spd[ind]+(usr[ind]/kappa)*(np.log(hh_in[0] /
u10n[ind] = (spd[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind] /
ref_ht)-psim[ind]))
u10n[u10n < 0] = np.nan
# u10n[u10n < 0] = np.nan
else:
u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind]/10) -
psim[ind]))
u10n[u10n < 0] = np.nan
# 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)])
......@@ -461,41 +450,29 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
ii = False
else:
ii = True
# if ((it == 3) or (it == 6) or (it == 10)):
# print('Method {}, # {} tol u10n: {} | t10n: {} | q10n: {} | '
# 'u*: {} | t*: {} | q*: {}'.format(meth, it,
# np.ma.max(d[0, ~np.isnan(d[0,:])]),
# np.ma.max(d[1, ~np.isnan(d[1,:])]),
# np.ma.max(d[2, ~np.isnan(d[2,:])]),
# np.ma.max(d[3, ~np.isnan(d[3,:])]),
# np.ma.max(d[4, ~np.isnan(d[4,:])]),
# np.ma.max(d[5, ~np.isnan(d[5,:])])),
# file=open('tol_mid.txt','a'))
logging.info('method %s | # of iterations:%s', meth, it)
logging.info('method %s | # of points that did not converge :%s', meth,
ind[0].size)
# calculate output parameters
rho = (0.34838*P)/(tv10n)
# rho = P*100./(287.1*(T+CtoK)*(1+0.61*qair)) # C35
t10n = t10n-(273.16+tlapse*ref_ht)
sensible = -rho*cp*usr*tsr
latent = -rho*lv*usr*qsr
if (meth == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
meth == "ERA5"):
if (gust == 1):
tau = rho*np.power(usr, 2)*(spd/wind)
else:
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(hh_in[0]/hh_out[0])-psim +
psim_calc(hh_out[0]/monob, meth)))
trefs = (Ta-(tsr/kappa)*(np.log(hh_in[1]/hh_out[1])-psit +
psit_calc(hh_out[0]/monob, meth)))
trefs = trefs-(273.16+tlapse*hh_out[1])
qrefs = (qair-(qsr/kappa)*(np.log(hh_in[2]/hh_out[2]) -
psit+psit_calc(hh_out[2]/monob, meth)))
res = np.zeros((27, len(spd)))
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)))
res = np.zeros((28, len(spd)))
res[0][:] = tau
res[1][:] = sensible
res[2][:] = latent
......@@ -512,15 +489,16 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
res[13][:] = usr
res[14][:] = psim
res[15][:] = psit
res[16][:] = u10n
res[17][:] = t10n
res[18][:] = tv10n
res[19][:] = q10n
res[20][:] = zo
res[21][:] = zot
res[22][:] = zoq
res[23][:] = urefs
res[24][:] = trefs
res[25][:] = qrefs
res[26][:] = itera
res[16][:] = psiq
res[17][:] = u10n
res[18][:] = t10n
res[19][:] = tv10n
res[20][:] = q10n
res[21][:] = zo
res[22][:] = zot
res[23][:] = zoq
res[24][:] = urefs
res[25][:] = trefs
res[26][:] = qrefs
res[27][:] = itera
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