Skip to content

GitLab

Commit d24e4c69 authored by sbiri's avatar sbiri
Browse files

Update flux_subs.py cdn_calc 

- changed wind speed limits for S80, LP82,
- added cool skin functions following ecmwf and Beljaars, 
- added warm layer correction function
parent beaf5d64
Hide whitespace changes
Inline Side-by-side
Showing with 386 additions and 146 deletions
flux_subs.py
View file @ d24e4c69
...@@ -5,18 +5,19 @@ from util_subs import (CtoK, kappa, gc, visc_air) ...@@ -5,18 +5,19 @@ from util_subs import (CtoK, kappa, gc, visc_air)
def cdn_calc(u10n, Ta, Tp, lat, meth="S80"): def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
""" Calculates neutral drag coefficient """
Calculates neutral drag coefficient
Parameters Parameters
---------- ----------
u10n : float u10n : float
neutral 10m wind speed (m/s) neutral 10m wind speed [m/s]
Ta : float Ta : float
air temperature (K) air temperature [K]
Tp : float Tp : float
wave period wave period
lat : float lat : float
latitude latitude [degE]
meth : str meth : str
Returns Returns
...@@ -28,17 +29,17 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"): ...@@ -28,17 +29,17 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
cdn = np.where(u10n <= 3, (0.61+0.567/u10n)*0.001, cdn = np.where(u10n <= 3, (0.61+0.567/u10n)*0.001,
(0.61+0.063*u10n)*0.001) (0.61+0.063*u10n)*0.001)
elif (meth == "LP82"): elif (meth == "LP82"):
cdn = np.where((u10n < 11) & (u10n >= 4), 1.2*0.001, cdn = np.where(u10n < 4, 1.14*0.001,
np.where((u10n <= 25) & (u10n >= 11), np.where((u10n < 11) & (u10n >= 4), 1.2*0.001,
(0.49+0.065*u10n)*0.001, 1.14*0.001)) (0.49+0.065*u10n)*0.001))
elif (meth == "S88" or meth == "UA" or meth == "ERA5" or meth == "C30" or elif (meth == "S88" or meth == "UA" or meth == "ecmwf" or meth == "C30" or
meth == "C35" or meth == "C40" or meth == "Beljaars"): meth == "C35" or meth == "C40" or meth == "Beljaars"):
cdn = cdn_from_roughness(u10n, Ta, None, lat, meth) cdn = cdn_from_roughness(u10n, Ta, None, lat, meth)
elif (meth == "YT96"): elif (meth == "YT96"):
# for u<3 same as S80 # for u<3 same as S80
cdn = np.where((u10n < 6) & (u10n >= 3), cdn = np.where((u10n < 6) & (u10n >= 3),
(0.29+3.1/u10n+7.7/u10n**2)*0.001, (0.29+3.1/u10n+7.7/u10n**2)*0.001,
np.where((u10n <= 26) & (u10n >= 6), np.where((u10n >= 6),
(0.60 + 0.070*u10n)*0.001, (0.61+0.567/u10n)*0.001)) (0.60 + 0.070*u10n)*0.001, (0.61+0.567/u10n)*0.001))
elif (meth == "LY04"): elif (meth == "LY04"):
cdn = np.where(u10n >= 0.5, cdn = np.where(u10n >= 0.5,
...@@ -51,17 +52,18 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"): ...@@ -51,17 +52,18 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"): def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
""" Calculates neutral drag coefficient from roughness length """
Calculates neutral drag coefficient from roughness length
Parameters Parameters
---------- ----------
u10n : float u10n : float
neutral 10m wind speed (m/s) neutral 10m wind speed [m/s]
Ta : float Ta : float
air temperature (K) air temperature [K]
Tp : float Tp : float
wave period wave period
lat : float lat : float [degE]
latitude latitude
meth : str meth : str
...@@ -92,16 +94,13 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"): ...@@ -92,16 +94,13 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
elif (meth == "C35"): elif (meth == "C35"):
a = 0.011*np.ones(Ta.shape) 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, 0.0017*u10n-0.0050) 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 zo = 0.11*visc_air(Ta)/usr+a*np.power(usr, 2)/g
elif (meth == "C40"): elif (meth == "C40"):
a = 0.011*np.ones(Ta.shape) a = 0.011*np.ones(Ta.shape)
a = np.where(u10n > 22, 0.0016*22-0.0035, 0.0016*u10n-0.0035) 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 zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr # surface roughness
elif ((meth == "ERA5" or meth == "Beljaars")): elif ((meth == "ecmwf" or meth == "Beljaars")):
# eq. (3.26) p.38 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 zo = 0.018*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
else: else:
...@@ -113,16 +112,17 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"): ...@@ -113,16 +112,17 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
def cd_calc(cdn, height, ref_ht, psim): def cd_calc(cdn, height, ref_ht, psim):
""" Calculates drag coefficient at reference height """
Calculates drag coefficient at reference height
Parameters Parameters
---------- ----------
cdn : float cdn : float
neutral drag coefficient neutral drag coefficient
height : float height : float
original sensor height (m) original sensor height [m]
ref_ht : float ref_ht : float
reference height (m) reference height [m]
psim : float psim : float
momentum stability function momentum stability function
...@@ -136,7 +136,8 @@ def cd_calc(cdn, height, ref_ht, psim): ...@@ -136,7 +136,8 @@ def cd_calc(cdn, height, ref_ht, psim):
def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
""" Calculates neutral heat and moisture exchange coefficients """
Calculates neutral heat and moisture exchange coefficients
Parameters Parameters
---------- ----------
...@@ -145,11 +146,11 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): ...@@ -145,11 +146,11 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
cdn : float cdn : float
neutral drag coefficient neutral drag coefficient
u10n : float u10n : float
neutral 10m wind speed (m/s) neutral 10m wind speed [m/s]
zo : float zo : float
surface roughness (m) surface roughness [m]
Ta : float Ta : float
air temperature (K) air temperature [K]
meth : str meth : str
Returns Returns
...@@ -163,10 +164,8 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): ...@@ -163,10 +164,8 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
cqn = np.ones(Ta.shape)*1.20*0.001 # from S88 cqn = np.ones(Ta.shape)*1.20*0.001 # from S88
ctn = np.ones(Ta.shape)*1.00*0.001 ctn = np.ones(Ta.shape)*1.00*0.001
elif (meth == "LP82"): elif (meth == "LP82"):
cqn = np.where((zol <= 0) & (u10n > 4) & (u10n < 14), 1.15*0.001, cqn = np.where((zol <= 0), 1.15*0.001, 1*0.001)
1*0.001) ctn = np.where((zol <= 0), 1.13*0.001, 0.66*0.001)
ctn = np.where((zol <= 0) & (u10n > 4) & (u10n < 25), 1.13*0.001,
0.66*0.001)
elif (meth == "LY04"): elif (meth == "LY04"):
cqn = 34.6*0.001*np.sqrt(cdn) cqn = 34.6*0.001*np.sqrt(cdn)
ctn = np.where(zol <= 0, 32.7*0.001*np.sqrt(cdn), 18*0.001*np.sqrt(cdn)) ctn = np.where(zol <= 0, 32.7*0.001*np.sqrt(cdn), 18*0.001*np.sqrt(cdn))
...@@ -183,6 +182,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): ...@@ -183,6 +182,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
elif (meth == "C30"): elif (meth == "C30"):
usr = np.sqrt(cdn*np.power(u10n, 2)) usr = np.sqrt(cdn*np.power(u10n, 2))
rr = zo*usr/visc_air(Ta) rr = zo*usr/visc_air(Ta)
# moisture roughness determined by the CLIMODE, GASEX and CBLAST data
zoq = np.where(5e-5/np.power(rr, 0.6) > 1.15e-4, 1.15e-4, zoq = np.where(5e-5/np.power(rr, 0.6) > 1.15e-4, 1.15e-4,
5e-5/np.power(rr, 0.6)) # moisture roughness 5e-5/np.power(rr, 0.6)) # moisture roughness
zot=zoq # temperature roughness zot=zoq # temperature roughness
...@@ -204,12 +204,9 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): ...@@ -204,12 +204,9 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
1.0e-4/np.power(rr, 0.55)) # temperature roughness 1.0e-4/np.power(rr, 0.55)) # temperature roughness
zoq = np.where(2.0e-5/np.power(rr,0.22) > 1.1e-4/np.power(rr,0.9), zoq = np.where(2.0e-5/np.power(rr,0.22) > 1.1e-4/np.power(rr,0.9),
1.1e-4/np.power(rr,0.9), 2.0e-5/np.power(rr,0.22)) 1.1e-4/np.power(rr,0.9), 2.0e-5/np.power(rr,0.22))
# moisture roughness determined by the CLIMODE, GASEX and CBLAST data
# zoq = np.where(5e-5/np.power(rr, 0.6) > 1.15e-4, 1.15e-4,
# 5e-5/np.power(rr, 0.6)) # moisture roughness as in C30
cqn = kappa**2/np.log(10/zo)/np.log(10/zoq) cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
ctn = kappa**2/np.log(10/zo)/np.log(10/zot) ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
elif (meth == "ERA5" or meth == "Beljaars"): elif (meth == "ecmwf" or meth == "Beljaars"):
# eq. (3.26) p.38 over sea IFS Documentation cy46r1 # eq. (3.26) p.38 over sea IFS Documentation cy46r1
usr = np.sqrt(cdn*np.power(u10n, 2)) usr = np.sqrt(cdn*np.power(u10n, 2))
zot = 0.40*visc_air(Ta)/usr zot = 0.40*visc_air(Ta)/usr
...@@ -223,7 +220,8 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"): ...@@ -223,7 +220,8 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, 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 """
Calculates heat and moisture exchange coefficients at reference height
Parameters Parameters
---------- ----------
...@@ -235,12 +233,12 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq): ...@@ -235,12 +233,12 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq):
neutral heat exchange coefficient neutral heat exchange coefficient
cqn : float cqn : float
neutral moisture exchange coefficient neutral moisture exchange coefficient
h_t : float ht : float
original temperature sensor height (m) original temperature sensor height [m]
h_q : float hq : float
original moisture sensor height (m) original moisture sensor height [m]
ref_ht : float ref_ht : float
reference height (m) reference height [m]
psit : float psit : float
heat stability function heat stability function
psiq : float psiq : float
...@@ -262,7 +260,8 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq): ...@@ -262,7 +260,8 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq):
def get_stabco(meth="S80"): def get_stabco(meth="S80"):
""" Gives the coefficients \\alpha, \\beta, \\gamma for stability functions """
Gives the coefficients \\alpha, \\beta, \\gamma for stability functions
Parameters Parameters
---------- ----------
...@@ -274,7 +273,7 @@ def get_stabco(meth="S80"): ...@@ -274,7 +273,7 @@ def get_stabco(meth="S80"):
""" """
alpha, beta, gamma = 0, 0, 0 alpha, beta, gamma = 0, 0, 0
if (meth == "S80" or meth == "S88" or meth == "LY04" or if (meth == "S80" or meth == "S88" or meth == "LY04" or
meth == "UA" or meth == "ERA5" or meth == "C30" or meth == "UA" or meth == "ecmwf" or meth == "C30" or
meth == "C35" or meth == "C40" or meth == "Beljaars"): meth == "C35" or meth == "C40" or meth == "Beljaars"):
alpha, beta, gamma = 16, 0.25, 5 # Smith 1980, from Dyer (1974) alpha, beta, gamma = 16, 0.25, 5 # Smith 1980, from Dyer (1974)
elif (meth == "LP82"): elif (meth == "LP82"):
...@@ -292,7 +291,8 @@ def get_stabco(meth="S80"): ...@@ -292,7 +291,8 @@ def get_stabco(meth="S80"):
def psim_calc(zol, meth="S80"): def psim_calc(zol, meth="S80"):
""" Calculates momentum stability function """
Calculates momentum stability function
Parameters Parameters
---------- ----------
...@@ -304,8 +304,8 @@ def psim_calc(zol, meth="S80"): ...@@ -304,8 +304,8 @@ def psim_calc(zol, meth="S80"):
------- -------
psim : float psim : float
""" """
if (meth == "ERA5"): if (meth == "ecmwf"):
psim = psim_era5(zol) psim = psim_ecmwf(zol)
elif (meth == "C30" or meth == "C35" or meth == "C40"): elif (meth == "C30" or meth == "C35" or meth == "C40"):
psim = psiu_26(zol, meth) psim = psiu_26(zol, meth)
elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17 elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
...@@ -318,7 +318,8 @@ def psim_calc(zol, meth="S80"): ...@@ -318,7 +318,8 @@ def psim_calc(zol, meth="S80"):
def psit_calc(zol, meth="S80"): def psit_calc(zol, meth="S80"):
""" Calculates heat stability function """
Calculates heat stability function
Parameters Parameters
---------- ----------
...@@ -331,9 +332,9 @@ def psit_calc(zol, meth="S80"): ...@@ -331,9 +332,9 @@ def psit_calc(zol, meth="S80"):
------- -------
psit : float psit : float
""" """
if (meth == "ERA5"): if (meth == "ecmwf"):
psit = np.where(zol < 0, psi_conv(zol, meth), psit = np.where(zol < 0, psi_conv(zol, meth),
psi_era5(zol)) psi_ecmwf(zol))
elif (meth == "C30" or meth == "C35" or meth == "C40"): elif (meth == "C30" or meth == "C35" or meth == "C40"):
psit = psit_26(zol) psit = psit_26(zol)
elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17 elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
...@@ -346,7 +347,8 @@ def psit_calc(zol, meth="S80"): ...@@ -346,7 +347,8 @@ def psit_calc(zol, meth="S80"):
def psi_Bel(zol): def psi_Bel(zol):
""" Calculates momentum/heat stability function """
Calculates momentum/heat stability function
Parameters Parameters
---------- ----------
...@@ -365,9 +367,10 @@ def psi_Bel(zol): ...@@ -365,9 +367,10 @@ def psi_Bel(zol):
# --------------------------------------------------------------------- # ---------------------------------------------------------------------
def psi_era5(zol): def psi_ecmwf(zol):
""" Calculates heat stability function for stable conditions """
for method ERA5 Calculates heat stability function for stable conditions
for method ecmwf
Parameters Parameters
---------- ----------
...@@ -386,7 +389,8 @@ def psi_era5(zol): ...@@ -386,7 +389,8 @@ def psi_era5(zol):
def psit_26(zol): def psit_26(zol):
""" Computes temperature structure function as in C35 """
Computes temperature structure function as in C35
Parameters Parameters
---------- ----------
...@@ -413,7 +417,8 @@ def psit_26(zol): ...@@ -413,7 +417,8 @@ def psit_26(zol):
def psi_conv(zol, meth): def psi_conv(zol, meth):
""" Calculates heat stability function for unstable conditions """
Calculates heat stability function for unstable conditions
Parameters Parameters
---------- ----------
...@@ -435,7 +440,8 @@ def psi_conv(zol, meth): ...@@ -435,7 +440,8 @@ def psi_conv(zol, meth):
def psi_stab(zol, meth): def psi_stab(zol, meth):
""" Calculates heat stability function for stable conditions """
Calculates heat stability function for stable conditions
Parameters Parameters
---------- ----------
...@@ -455,8 +461,9 @@ def psi_stab(zol, meth): ...@@ -455,8 +461,9 @@ def psi_stab(zol, meth):
# --------------------------------------------------------------------- # ---------------------------------------------------------------------
def psim_era5(zol): def psim_ecmwf(zol):
""" Calculates momentum stability function for method ERA5 """
Calculates momentum stability function for method ecmwf
Parameters Parameters
---------- ----------
...@@ -468,7 +475,7 @@ def psim_era5(zol): ...@@ -468,7 +475,7 @@ def psim_era5(zol):
psim : float psim : float
""" """
# eq (3.20, 3.22) p. 37 IFS Documentation cy46r1 # eq (3.20, 3.22) p. 37 IFS Documentation cy46r1
coeffs = get_stabco("ERA5") coeffs = get_stabco("ecmwf")
alpha, beta = coeffs[0], coeffs[1] alpha, beta = coeffs[0], coeffs[1]
xtmp = np.power(1-alpha*zol, beta) xtmp = np.power(1-alpha*zol, beta)
a, b, c, d = 1, 2/3, 5, 0.35 a, b, c, d = 1, 2/3, 5, 0.35
...@@ -480,7 +487,8 @@ def psim_era5(zol): ...@@ -480,7 +487,8 @@ def psim_era5(zol):
def psiu_26(zol, meth): def psiu_26(zol, meth):
""" Computes velocity structure function C35 """
Computes velocity structure function C35
Parameters Parameters
---------- ----------
...@@ -524,7 +532,8 @@ def psiu_26(zol, meth): ...@@ -524,7 +532,8 @@ def psiu_26(zol, meth):
def psim_conv(zol, meth): def psim_conv(zol, meth):
""" Calculates momentum stability function for unstable conditions """
Calculates momentum stability function for unstable conditions
Parameters Parameters
---------- ----------
...@@ -547,7 +556,8 @@ def psim_conv(zol, meth): ...@@ -547,7 +556,8 @@ def psim_conv(zol, meth):
def psim_stab(zol, meth): def psim_stab(zol, meth):
""" Calculates momentum stability function for stable conditions """
Calculates momentum stability function for stable conditions
Parameters Parameters
---------- ----------
...@@ -567,46 +577,46 @@ def psim_stab(zol, meth): ...@@ -567,46 +577,46 @@ def psim_stab(zol, meth):
# --------------------------------------------------------------------- # ---------------------------------------------------------------------
def get_skin(sst, qsea, rho, Rl, Rs, Rnl, cp, lv, tkt, usr, tsr, qsr, lat): def cs_C35(sst, qsea, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, lat):
""" Computes cool skin """
Computes cool skin
following COARE3.5 (Fairall et al. 1996, Edson et al. 2013)
Parameters Parameters
---------- ----------
sst : float sst : float
sea surface temperature ($^\circ$\,C) sea surface temperature [K]
qsea : float qsea : float
specific humidity over sea (g/kg) specific humidity over sea [g/kg]
rho : float rho : float
density of air (kg/m^3) density of air [kg/m^3]
Rl : float
downward longwave radiation (W/m^2)
Rs : float Rs : float
downward shortwave radiation (W/m^2) downward shortwave radiation [Wm-2]
Rnl : float Rnl : float
upwelling IR radiation (W/m^2) upwelling IR radiation [Wm-2]
cp : float cp : float
specific heat of air at constant pressure specific heat of air at constant pressure [J/K/kg]
lv : float lv : float
latent heat of vaporization latent heat of vaporization [J/kg]
tkt : float delta : float
cool skin thickness cool skin thickness [m]
usr : float usr : float
friction velocity friction velocity [m/s]
tsr : float tsr : float
star temperature star temperature [K]
qsr : float qsr : float
star humidity star humidity [g/kg]
lat : float lat : float
latitude latitude [degE]
Returns Returns
------- -------
dter : float dter : float
cool-skin temperature depression cool skin correction [K]
dqer : float dqer : float
cool-skin humidity depression humidity corrction [g/kg]
tkt : float delta : float
cool skin thickness cool skin thickness [m]
""" """
# coded following Saunders (1967) with lambda = 6 # coded following Saunders (1967) with lambda = 6
g = gc(lat, None) g = gc(lat, None)
...@@ -615,50 +625,280 @@ def get_skin(sst, qsea, rho, Rl, Rs, Rnl, cp, lv, tkt, usr, tsr, qsr, lat): ...@@ -615,50 +625,280 @@ def get_skin(sst, qsea, rho, Rl, Rs, Rnl, cp, lv, tkt, usr, tsr, qsr, lat):
# ************ cool skin constants ******* # ************ cool skin constants *******
# density of water, specific heat capacity of water, water viscosity, # density of water, specific heat capacity of water, water viscosity,
# thermal conductivity of water # thermal conductivity of water
# rhow, cpw, visw, tcw = 1025, 4190, 1e-6, 0.6
rhow, cpw, visw, tcw = 1022, 4000, 1e-6, 0.6 rhow, cpw, visw, tcw = 1022, 4000, 1e-6, 0.6
Al = 2.1e-5*np.power(sst+3.2, 0.79) aw = 2.1e-5*np.power(np.maximum(sst+3.2, 0), 0.79)
be = 0.026
bigc = 16*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(rho, 2)) bigc = 16*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(rho, 2))
wetc = 0.622*lv*qsea/(287.1*np.power(sst+273.16, 2)) wetc = 0.622*lv*qsea/(287.1*np.power(sst+273.16, 2))
Rns = 0.945*Rs # albedo correction Rns = 0.945*Rs # albedo correction
hsb = -rho*cp*usr*tsr shf = -rho*cp*usr*tsr
hlb = -rho*lv*usr*qsr lhf = -rho*lv*usr*qsr
qout = Rnl+hsb+hlb Qnsol = Rnl+shf+lhf
dels = Rns*(0.065+11*tkt-6.6e-5/tkt*(1-np.exp(-tkt/8.0e-4))) fs = 0.065+11*delta-6.6e-5/delta*(1-np.exp(-delta/8.0e-4))
qcol = qout-dels qcol = Qnsol-Rns*fs
alq = Al*qcol+be*hlb*cpw/lv alq = aw*qcol+0.026*lhf*cpw/lv
xlamx = 6*np.ones(sst.shape) xlamx = 6*np.ones(sst.shape)
xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6) xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6)
tkt = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), delta = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr),
np.where(xlamx*visw/(np.sqrt(rho/rhow)*usr) > 0.01, 0.01, np.where(xlamx*visw/(np.sqrt(rho/rhow)*usr) > 0.01, 0.01,
xlamx*visw/(np.sqrt(rho/rhow)*usr))) xlamx*visw/(np.sqrt(rho/rhow)*usr)))
dter = qcol*tkt/tcw dter = qcol*delta/tcw
dqer = wetc*dter dqer = wetc*dter
return dter, dqer, tkt return dter, dqer, delta
# ---------------------------------------------------------------------
def delta(aw, Qnsol, usr, lat):
"""
Computes the thickness (m) of the viscous skin layer.
Based on Fairall et al., 1996 and cited in IFS Documentation Cy46r1
eq. 8.155 p. 164
Parameters
----------
aw : float
thermal expansion coefficient of sea-water [1/K]
Qnsol : float
part of the net heat flux actually absorbed in the warm layer [W/m^2]
usr : float
friction velocity in the air (u*) [m/s]
lat : float
latitude [degE]
Returns
-------
delta : float
the thickness (m) of the viscous skin layer
"""
rhow, visw, tcw = 1025, 1e-6, 0.6
# u* in the water
usr_w = np.maximum(usr, 1e-4)*np.sqrt(1.2/rhow) # rhoa=1.2
rcst_cs = 16*gc(lat)*np.power(visw, 3)/np.power(tcw, 2)
lm = 6/np.power(1+np.power(np.maximum(Qnsol*aw*rcst_cs /
np.power(usr_w, 4), 0), 3/4),
1/3)
ztmp = visw/usr_w
delta = np.where(Qnsol > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
return delta
# ---------------------------------------------------------------------
def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, lat):
"""
cool skin adjustment based on IFS Documentation cy46r1
Parameters
----------
rho : float
density of air [kg/m^3]
Rs : float
downward solar radiation [Wm-2]
Rnl : float
net thermal radiaion [Wm-2]
cp : float
specific heat of air at constant pressure [J/K/kg]
lv : float
latent heat of vaporization [J/kg]
usr : float
friction velocity [m/s]
tsr : float
star temperature [K]
qsr : float
star humidity [g/kg]
sst : float
sea surface temperature [K]
lat : float
latitude [degE]
Returns
-------
dtc : float
cool skin temperature correction [K]
"""
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+CtoK
aw = np.maximum(1e-5, 1e-5*(sst-CtoK))
Rns = 0.945*Rs # (1-0.066)*Rs net solar radiation (albedo correction)
shf = -rho*cp*usr*tsr
lhf = -rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl # eq. 8.152
d = delta(aw, Qnsol, usr, lat)
for jc in range(10): # because implicit in terms of delta...
# # fraction of the solar radiation absorbed in layer delta eq. 8.153
# and Eq.(5) Zeng & Beljaars, 2005
fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/8e-4))
Q = Qnsol-fs*Rns
d = delta(aw, Q, usr, lat)
dtc = Q*d/0.6 # (rhow*cw*kw)eq. 8.151
return dtc
# --------------------------------------------------------------------- # ---------------------------------------------------------------------
def wl_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, skt, dtc, lat):
"""
warm layer correction following IFS Documentation cy46r1
and aerobulk (Brodeau et al., 2016)
Parameters
----------
rho : float
density of air [kg/m^3]
Rs : float
downward solar radiation [Wm-2]
Rnl : TYPE
net thermal radiation [Wm-2]
cp : float
specific heat of air at constant pressure [J/K/kg]
lv : float
latent heat of vaporization [J/kg]
usr : float
friction velocity [m/s]
tsr : float
star temperature [K]
qsr : float
star humidity [g/kg]
sst : float
bulk sst [K]
skt : float
skin sst from previous step [K]
dtc : float
cool skin correction [K]
lat : float
latitude [degE]
Returns
-------
dtwl : float
warm layer correction [K]
"""
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+CtoK
rhow, cpw, visw, rd0 = 1025, 4190, 1e-6, 3
Rns = 0.945*Rs
## Previous value of dT / warm-layer, adapted to depth:
aw = 2.1e-5*np.power(np.maximum(sst-CtoK+3.2, 0), 0.79)
# thermal expansion coefficient of sea-water (SST accurate enough#)
# *** Rd = Fraction of solar radiation absorbed in warm layer (-)
a1, a2, a3 = 0.28, 0.27, 0.45
b1, b2, b3 = -71.5, -2.8, -0.06 # [m-1]
Rd = 1-(a1*np.exp(b1*rd0)+a2*np.exp(b2*rd0)+a3*np.exp(b3*rd0))
shf = -rho*cp*usr*tsr
lhf = -rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl
usrw = np.maximum(usr, 1e-4 )*np.sqrt(1.2/rhow) # u* in the water
zc3 = rd0*kappa*gc(lat)/np.power(1.2/rhow, 3/2)
zc4 = (0.3+1)*kappa/rd0
zc5 = (0.3+1)/(0.3*rd0)
for jwl in range(10): # itteration to solve implicitely eq. for warm layer
dsst = skt-sst+dtc
## Buoyancy flux and stability parameter (zdl = -z/L) in water
ZSRD = (Qnsol-Rns*Rd)/(rhow*cpw)
ztmp = np.maximum(dsst, 0)
zdl = np.where(ZSRD > 0, 2*(np.power(usrw, 2) *
np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))+ZSRD,
np.power(usrw, 2)*np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))
usr = np.maximum(usr, 1e-4 )
zdL = zc3*aw*zdl/np.power(usr, 3)
## Stability function Phi_t(-z/L) (zdL is -z/L) :
zphi = np.where(zdL > 0, (1+(5*zdL+4*np.power(zdL, 2)) /
(1+3*zdL+0.25*np.power(zdL, 2)) +
2/np.sqrt(1-16*(-np.abs(zdL)))),
1/np.sqrt(1-16*(-np.abs(zdL))))
zz = zc4*(usrw)/zphi
zz = np.maximum(np.abs(zz) , 1e-4)*np.sign(zz)
dtwl = np.maximum(0, (zc5*ZSRD)/zz)
return dtwl
#----------------------------------------------------------------------
def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, lat, Qs):
"""
cool skin adjustment based on Beljaars (1997)
air-sea interaction in the ECMWF model
Parameters
----------
rho : float
density of air [kg/m^3]
Rs : float
downward solar radiation [Wm-2]
Rnl : float
net thermal radiaion [Wm-2]
cp : float
specific heat of air at constant pressure [J/K/kg]
lv : float
latent heat of vaporization [J/kg]
usr : float
friction velocity [m/s]
tsr : float
star temperature [K]
qsr : float
star humidity [g/kg]
lat : float
latitude [degE]
Qs : float
radiation balance
Returns
-------
Qs : float
radiation balance
dtc : float
cool skin temperature correction [K]
"""
g = gc(lat, None)
tcw = 0.6 # thermal conductivity of water (at 20C) [W/m/K]
visw = 1e-6 # kinetic viscosity of water [m^2/s]
rhow = 1025 # Density of sea-water [kg/m^3]
cpw = 4190 # specific heat capacity of water
aw = 3e-4 # thermal expansion coefficient [K-1]
Rns = 0.945*Rs # net solar radiation (albedo correction)
shf = -rho*cp*usr*tsr
lhf = -rho*lv*usr*qsr
Q = Rnl+shf+lhf+Qs
xt = (16*Q*g*aw*cpw*np.power(rhow*visw, 3))/(np.power(usr, 4) *
np.power(rho*tcw, 2))
xt1 = 1+np.power(xt, 3/4)
# Saunders const eq. 22
ls = np.where(Q > 0, 6/np.power(xt1, 1/3), 6)
delta = np.where(Q > 0 , (ls*visw)/(np.sqrt(rho/rhow)*usr),
np.where((ls*visw)/(np.sqrt(rho/rhow)*usr) > 0.01, 0.01,
(ls*visw)/(np.sqrt(rho/rhow)*usr))) # eq. 21
# fraction of the solar radiation absorbed in layer delta
fc = 0.065+11*delta-6.6e-5*(1-np.exp(-delta/0.0008))/delta
Qs = fc*Rns
Q = Rnl+shf+lhf+Qs
dtc = Q*delta/tcw
return Qs, dtc
#----------------------------------------------------------------------
def get_gust(beta, Ta, usr, tsrv, zi, lat): def get_gust(beta, Ta, usr, tsrv, zi, lat):
""" Computes gustiness """
Computes gustiness
Parameters Parameters
---------- ----------
beta : float beta : float
constant constant
Ta : float Ta : float
air temperature (K) air temperature [K]
usr : float usr : float
friction velocity (m/s) friction velocity [m/s]
tsrv : float tsrv : float
star virtual temperature of air (K) star virtual temperature of air [K]
zi : int zi : int
scale height of the boundary layer depth (m) scale height of the boundary layer depth [m]
lat : float lat : float
latitude latitude
Returns Returns
------- -------
ug : float ug : float [m/s]
""" """
if (np.nanmax(Ta) < 200): # convert to K if in Celsius if (np.nanmax(Ta) < 200): # convert to K if in Celsius
Ta = Ta+273.16 Ta = Ta+273.16
...@@ -670,8 +910,8 @@ def get_gust(beta, Ta, usr, tsrv, zi, lat): ...@@ -670,8 +910,8 @@ def get_gust(beta, Ta, usr, tsrv, zi, lat):
# --------------------------------------------------------------------- # ---------------------------------------------------------------------
def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, def get_L(L, lat, usr, tsr, qsr, t10n, h_in, Ta, sst, qair, qsea, q10n,
dt, dtv, dq, zo, wind, monob, meth): wind, monob, meth):
""" """
calculates Monin-Obukhov length and virtual star temperature calculates Monin-Obukhov length and virtual star temperature
...@@ -679,8 +919,8 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -679,8 +919,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
"S80" : default for S80, S88, LP82, YT96 and LY04 "S80" : default for S80, S88, LP82, YT96 and LY04
"ERA5" : following ERA5 (IFS Documentation cy46r1), default for ERA5 "ecmwf" : following ecmwf (IFS Documentation cy46r1), default for ecmwf
lat : float lat : float
latitude latitude
usr : float usr : float
...@@ -691,33 +931,25 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -691,33 +931,25 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
star specific humidity (g/kg) star specific humidity (g/kg)
t10n : float t10n : float
neutral temperature at 10m (K) neutral temperature at 10m (K)
tv10n : float
neutral virtual temperature at 10m (K)
qair : float
air specific humidity (g/kg)
h_in : float h_in : float
sensor heights (m) sensor heights (m)
T : float
air temperature (K)
Ta : float Ta : float
air temperature (K) air temperature (K)
th : float
potential temperature (K)
tv : float
virtual temperature (K)
sst : float sst : float
sea surface temperature (K) sea surface temperature (K)
dt : float qair : float
temperature difference (K) air specific humidity (g/kg)
dq : float qsea : float
specific humidity difference (g/kg) specific humidity at sea surface(g/kg)
q10n : float
neutral specific humidity at 10m (g/kg)
wind : float wind : float
wind speed (m/s) wind speed (m/s)
monob : float monob : float
Monin-Obukhov length from previous iteration step (m) Monin-Obukhov length from previous iteration step (m)
meth : str meth : str
bulk parameterisation method option: "S80", "S88", "LP82", "YT96", bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
"UA", "LY04", "C30", "C35", "C40", "ERA5" "UA", "LY04", "C30", "C35", "C40", "ecmwf", "Beljaars"
Returns Returns
------- -------
...@@ -729,13 +961,17 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -729,13 +961,17 @@ 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 == "S80"): if (L == "S80"):
tsrv = tsr+0.61*t10n*qsr # as in NCAR, LY04
monob = ((tv10n*np.power(usr, 2))/(g*kappa*tsrv)) tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
monob = np.where(np.fabs(monob) < 1, np.where(monob < 0, -1, 1), monob) temp = (g*kappa*tsrv /
elif (L == "ERA5"): np.maximum(np.power(usr, 2)*Ta*(1+0.6077*qair), 1e-9))
tsrv = tsr+0.61*t10n*qsr temp = np.minimum(np.abs(temp), 200)*np.sign(temp)
Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_in[0])+0.61*dq)) / temp = np.minimum(np.abs(temp), 10/h_in[0])*np.sign(temp)
np.power(wind, 2)) monob = 1/np.copy(temp)
elif (L == "ecmwf"):
tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
Rb = ((g*h_in[0]/(Ta*(1+0.61*qair))) *
((t10n*(1+0.61*q10n)-sst*(1+0.61*qsea))/np.power(wind, 2)))
zo = (0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g) zo = (0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g)
zot = 0.40*visc_air(Ta)/usr zot = 0.40*visc_air(Ta)/usr
zol = (Rb*(np.power(np.log((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) / zol = (Rb*(np.power(np.log((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) /
...@@ -744,56 +980,60 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst, ...@@ -744,56 +980,60 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
(np.log((h_in[0]+zo)/zot) - (np.log((h_in[0]+zo)/zot) -
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
return tsrv, monob return tsrv, monob
#------------------------------------------------------------------------------ #------------------------------------------------------------------------------
def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq, def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, dtwl, ct, cq,
cskin, meth): cskin, wl, meth):
""" """
calculates star wind speed, temperature and specific humidity calculates star wind speed, temperature and specific humidity
Parameters Parameters
---------- ----------
h_in : float h_in : float
sensor heights (m) sensor heights [m]
monob : float monob : float
M-O length (m) M-O length [m]
wind : float wind : float
wind speed (m/s) wind speed [m/s]
zo : float zo : float
momentum roughness length (m) momentum roughness length [m]
zot : float zot : float
temperature roughness length (m) temperature roughness length [m]
zoq : float zoq : float
moisture roughness length (m) moisture roughness length [m]
dt : float dt : float
temperature difference (K) temperature difference [K]
dq : float dq : float
specific humidity difference (g/kg) specific humidity difference [g/kg]
dter : float dter : float
cskin temperature adjustment (K) cskin temperature adjustment [K]
dqer : float dqer : float
cskin q adjustment (g/kg) cskin q adjustment [g/kg]
dtwl : float
warm layer temperature adjustment [K]
ct : float ct : float
temperature exchange coefficient temperature exchange coefficient
cq : float cq : float
moisture exchange coefficient moisture exchange coefficient
cskin : int cskin : int
cool skin adjustment switch cool skin adjustment switch
wl : int
warm layer correction switch
meth : str meth : str
bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA", bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA",
"LY04", "C30", "C35", "C40", "ERA5" "LY04", "C30", "C35", "C40", "ecmwf", "Beljaars"
Returns Returns
------- -------
usr : float usr : float
friction wind speed (m/s) friction wind speed [m/s]
tsr : float tsr : float
star temperature (K) star temperature [K]
qsr : float qsr : float
star specific humidity (g/kg) star specific humidity [g/kg]
""" """
if (meth == "UA"): if (meth == "UA"):
...@@ -814,19 +1054,19 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq, ...@@ -814,19 +1054,19 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq,
5*np.log(h_in[0]/monob) + 5*np.log(h_in[0]/monob) +
h_in[0]/monob-1)))) 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-dtwl*wl) /
(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, kappa*(dt+dter*cskin) / np.where(h_in[1]/monob < 0, kappa*(dt+dter*cskin-dtwl*wl) /
(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 <= 1, np.where(h_in[1]/monob <= 1,
kappa*(dt+dter*cskin) / kappa*(dt+dter*cskin-dtwl*wl) /
(np.log(h_in[1]/zot) + (np.log(h_in[1]/zot) +
5*h_in[1]/monob-5*zot/monob), 5*h_in[1]/monob-5*zot/monob),
kappa*(dt+dter*cskin) / kappa*(dt+dter*cskin-dtwl*wl) /
(np.log(monob/zot)+5 - (np.log(monob/zot)+5 -
5*zot/monob+5*np.log(h_in[1]/monob) + 5*zot/monob+5*np.log(h_in[1]/monob) +
h_in[1]/monob-1)))) h_in[1]/monob-1))))
...@@ -850,13 +1090,13 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq, ...@@ -850,13 +1090,13 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq,
h_in[2]/monob-1)))) h_in[2]/monob-1))))
elif (meth == "C30" or meth == "C35" or meth == "C40"): elif (meth == "C30" or meth == "C35" or meth == "C40"):
usr = (wind*kappa/(np.log(h_in[0]/zo)-psiu_26(h_in[0]/monob, meth))) usr = (wind*kappa/(np.log(h_in[0]/zo)-psiu_26(h_in[0]/monob, meth)))
tsr = ((dt+dter*cskin)*(kappa/(np.log(h_in[1]/zot) - tsr = ((dt+dter*cskin-dtwl*wl)*(kappa/(np.log(h_in[1]/zot) -
psit_26(h_in[1]/monob)))) psit_26(h_in[1]/monob))))
qsr = ((dq+dqer*cskin)*(kappa/(np.log(h_in[2]/zoq) - qsr = ((dq+dqer*cskin)*(kappa/(np.log(h_in[2]/zoq) -
psit_26(h_in[2]/monob)))) psit_26(h_in[2]/monob))))
else: else:
usr = (wind*kappa/(np.log(h_in[0]/zo)-psim_calc(h_in[0]/monob, meth))) usr = (wind*kappa/(np.log(h_in[0]/zo)-psim_calc(h_in[0]/monob, meth)))
tsr = ct*wind*(dt+dter*cskin)/usr tsr = ct*wind*(dt+dter*cskin-dtwl*wl)/usr
qsr = cq*wind*(dq+dqer*cskin)/usr qsr = cq*wind*(dq+dqer*cskin)/usr
return usr, tsr, qsr return usr, tsr, qsr
# --------------------------------------------------------------------- # ---------------------------------------------------------------------
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