From 5976a8b6589d3947d153e236a4f09e9cc47962c6 Mon Sep 17 00:00:00 2001
From: sbiri <sbiri@noc.ac.uk>
Date: Tue, 16 Mar 2021 07:16:09 +0000
Subject: [PATCH] C40 option removed
---
AirSeaFluxCode.py | 23 +++++++++++------------
flux_subs.py | 44 +++++++++++++++-----------------------------
get_init.py | 22 ++++++++++------------
toy_ASFC.py | 10 +++++-----
4 files changed, 41 insertions(+), 58 deletions(-)
diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index 377ad34..b7ef1d0 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -58,7 +58,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
default for UA, ecmwf [1, 1, 1000]
default else [1, 1.2, 800]
meth : str
- "S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35", "C40",
+ "S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
"ecmwf", "Beljaars"
qmeth : str
is the saturation evaporation method to use amongst
@@ -86,17 +86,17 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
Returns
-------
res : array that contains
- 1. momentum flux (N/m^2)
- 2. sensible heat (W/m^2)
- 3. latent heat (W/m^2)
- 4. Monin-Obhukov length (mb)
+ 1. momentum flux (N/m^2)
+ 2. sensible heat (W/m^2)
+ 3. latent heat (W/m^2)
+ 4. Monin-Obhukov length (m)
5. drag coefficient (cd)
6. neutral drag coefficient (cdn)
- 7. heat exhange coefficient (ct)
- 8. neutral heat exhange coefficient (ctn)
+ 7. heat exchange coefficient (ct)
+ 8. neutral heat exchange coefficient (ctn)
9. moisture exhange coefficient (cq)
- 10. neutral moisture exhange coefficient (cqn)
- 11. star virtual temperature (tsrv)
+ 10. neutral moisture exchange coefficient (cqn)
+ 11. star virtual temperatcure (tsrv)
12. star temperature (tsr)
13. star specific humidity (qsr)
14. star wind speed (usr)
@@ -110,7 +110,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
22. surface roughness length (zo)
23. heat roughness length (zot)
24. moisture roughness length (zoq)
- 25. velocity at reference height (uref)
+ 25. wind speed at reference height (uref)
26. temperature at reference height (tref)
27. specific humidity at reference height (qref)
28. number of iterations until convergence
@@ -349,8 +349,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
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")):
+ elif (gust[0] == 1 and (meth == "C30" or meth == "C35")):
wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
np.power(get_gust(gust[1], Ta[ind], usr[ind],
tsrv[ind], gust[2], lat[ind]), 2))
diff --git a/flux_subs.py b/flux_subs.py
index bede98c..f8e079b 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -6,7 +6,7 @@ from util_subs import (CtoK, kappa, gc, visc_air)
def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
"""
- Calculates 10m neutral drag coefficient
+ Calculates neutral drag coefficient
Parameters
----------
@@ -33,7 +33,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
np.where((u10n < 11) & (u10n >= 4), 1.2*0.001,
(0.49+0.065*u10n)*0.001))
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 == "Beljaars"):
cdn = cdn_from_roughness(u10n, Ta, None, lat, meth)
elif (meth == "YT96"):
# for u<3 YT96 convert usr in eq. 21 to cdn
@@ -56,7 +56,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
"""
- Calculates 10m neutral drag coefficient from roughness length
+ Calculates neutral drag coefficient from roughness length
Parameters
----------
@@ -99,10 +99,6 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
a = 0.011*np.ones(Ta.shape)
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"):
- a = 0.011*np.ones(Ta.shape)
- 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 == "ecmwf" or meth == "Beljaars")):
# eq. (3.26) p.38 over sea IFS Documentation cy46r1
zo = 0.018*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
@@ -140,7 +136,7 @@ def cd_calc(cdn, hin, hout, psim):
def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
"""
- Calculates 10m neutral heat and moisture exchange coefficients
+ Calculates neutral heat and moisture exchange coefficients
Parameters
----------
@@ -198,16 +194,6 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
zot=zoq # temperature roughness
cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
- elif (meth == "C40"):
- usr = np.sqrt(cdn*np.power(u10n, 2))
- rr = zo*usr/visc_air(Ta)
- zot = np.where(1.0e-4/np.power(rr, 0.55) > 2.4e-4/np.power(rr, 1.2),
- 2.4e-4/np.power(rr, 1.2),
- 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),
- 1.1e-4/np.power(rr,0.9), 2.0e-5/np.power(rr,0.22))
- cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
- ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
elif (meth == "ecmwf" or meth == "Beljaars"):
# eq. (3.26) p.38 over sea IFS Documentation cy46r1
usr = np.sqrt(cdn*np.power(u10n, 2))
@@ -236,11 +222,11 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, hout, psit, psiq):
cqn : float
neutral moisture exchange coefficient
ht : float
- original temperature height [m]
+ original temperature sensor height [m]
hq : float
- original moisture height [m]
+ original moisture sensor height [m]
hout : float
- reference height [m]
+ reference height [m]
psit : float
heat stability function
psiq : float
@@ -276,7 +262,7 @@ def get_stabco(meth="S80"):
alpha, beta, gamma = 0, 0, 0
if (meth == "S80" or meth == "S88" or meth == "LY04" or
meth == "UA" or meth == "ecmwf" or meth == "C30" or
- meth == "C35" or meth == "C40" or meth == "Beljaars"):
+ meth == "C35" or meth == "Beljaars"):
alpha, beta, gamma = 16, 0.25, 5 # Smith 1980, from Dyer (1974)
elif (meth == "LP82"):
alpha, beta, gamma = 16, 0.25, 7
@@ -308,7 +294,7 @@ def psim_calc(zol, meth="S80"):
"""
if (meth == "ecmwf"):
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)
elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
psim = np.where(zol < 0, psim_conv(zol, meth), psi_Bel(zol))
@@ -337,7 +323,7 @@ def psit_calc(zol, meth="S80"):
if (meth == "ecmwf"):
psit = np.where(zol < 0, psi_conv(zol, meth),
psi_ecmwf(zol))
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
+ elif (meth == "C30" or meth == "C35"):
psit = psit_26(zol)
elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
psit = np.where(zol < 0, psi_conv(zol, meth), psi_Bel(zol))
@@ -514,7 +500,7 @@ def psiu_26(zol, meth):
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"):
+ elif (meth == "C35"):
dzol = np.where(0.35*zol > 50, 50, 0.35*zol) # stable
a, b, c, d = 0.7, 3/4, 5, 0.35
psi = np.where(zol > 0, -(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d),
@@ -949,7 +935,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
Monin-Obukhov length from previous iteration step (m)
meth : str
bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
- "UA", "LY04", "C30", "C35", "C40", "ecmwf", "Beljaars"
+ "UA", "LY04", "C30", "C35", "ecmwf", "Beljaars"
Returns
-------
@@ -988,7 +974,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
(np.log((hin[1]+zo)/zot) -
psit_calc((hin[1]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
- monob = hin[1]/zol
+ monob = hin[1]/zol
return tsrv, monob, Rb
#------------------------------------------------------------------------------
@@ -1032,7 +1018,7 @@ def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, dtwl, ct, cq,
warm layer correction switch
meth : str
bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA",
- "LY04", "C30", "C35", "C40", "ecmwf", "Beljaars"
+ "LY04", "C30", "C35", "ecmwf", "Beljaars"
Returns
-------
@@ -1096,7 +1082,7 @@ def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, dtwl, ct, cq,
(np.log(monob/zoq)+5-5*zoq/monob +
5*np.log(hin[2]/monob) +
hin[2]/monob-1))))
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
+ elif (meth == "C30" or meth == "C35"):
usr = (wind*kappa/(np.log(hin[0]/zo)-psiu_26(hin[0]/monob, meth)))
tsr = ((dt+dter*cskin-dtwl*wl)*(kappa/(np.log(hin[1]/zot) -
psit_26(hin[1]/monob))))
diff --git a/get_init.py b/get_init.py
index 8f2b49f..73b18ff 100644
--- a/get_init.py
+++ b/get_init.py
@@ -1,8 +1,8 @@
import numpy as np
import sys
-def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, meth,
- qmeth):
+def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
+ meth, qmeth):
"""
Checks initial input values and sets defaults if needed
@@ -17,7 +17,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
sea surface temperature in K
lat : float
latitude (deg), default 45deg
- hum : float
+ hum : array
relative humidity, if None is set to 80%
P : float
air pressure (hPa), default 1013hPa
@@ -40,14 +40,14 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
default else [1, 1.2, 800]
L : int
Monin-Obukhov length definition options
- tol : float
+ tol : array
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, 1e-3, 0.1, 0.1]
meth : str
- "S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ecmwf",
+ "S80","S88","LP82","YT96","UA","LY04","C30","C35","ecmwf",
"Beljaars"
qmeth : str
is the saturation evaporation method to use amongst
@@ -90,7 +90,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
sys.exit("input dtype of spd, T and SST should be float")
# if input values are nan break
if meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
- "C40", "ecmwf", "Beljaars"]:
+ "ecmwf", "Beljaars"]:
sys.exit("unknown method")
if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler",
"GoffGratch", "WMO", "MagnusTetens", "Buck", "Buck2",
@@ -119,11 +119,10 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
or meth == "YT96" or meth == "UA" or
meth == "LY04")):
cskin = 0
- elif ((cskin == None) and (meth == "C30" or meth == "C35" or meth == "C40"
+ elif ((cskin == None) and (meth == "C30" or meth == "C35"
or meth == "ecmwf" or meth == "Beljaars")):
cskin = 1
- if ((skin == None) and (meth == "C30" or meth == "C35"
- or meth == "C40")):
+ if ((skin == None) and (meth == "C30" or meth == "C35")):
skin = "C35"
elif ((skin == None) and (meth == "ecmwf")):
skin = "ecmwf"
@@ -131,8 +130,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
skin = "Beljaars"
if (wl == None):
wl = 0
- if (np.all(gust == None) and (meth == "C30" or meth == "C35" or
- meth == "C40")):
+ if (np.all(gust == None) and (meth == "C30" or meth == "C35")):
gust = [1, 1.2, 600]
elif (np.all(gust == None) and (meth == "UA" or meth == "ecmwf" or
meth == "Beljaars")):
@@ -148,7 +146,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
if ((L == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96" or meth == "LY04" or
meth == "UA" or meth == "C30" or meth == "C35"
- or meth == "C40" or meth == "Beljaars")):
+ or meth == "Beljaars")):
L = "S80"
elif ((L == None) and (meth == "ecmwf")):
L = "ecmwf"
diff --git a/toy_ASFC.py b/toy_ASFC.py
index abe3b37..31318e9 100644
--- a/toy_ASFC.py
+++ b/toy_ASFC.py
@@ -493,7 +493,7 @@ start_time = time.perf_counter()
inF = input("Give input file name (data_all.csv or era5_r360x180.nc): \n")
meth = input("Give prefered method: \n")
while meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
- "C40", "ecmwf","Beljaars"]:
+ "ecmwf","Beljaars"]:
print("method unknown")
meth = input("Give prefered method: \n")
else:
@@ -519,8 +519,8 @@ if (cskinIn == ''):
meth == "LY04")):
cskinIn = 0
ext = ext+'noskin_'
- elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
- or meth == "ecmwf" or meth == "Beljaars")):
+ elif ((cskinIn == None) and (meth == "C30" or meth == "C35"
+ or meth == "ecmwf" or meth == "Beljaars")):
cskinIn = 1
ext = ext+'skin_'
else:
@@ -589,10 +589,10 @@ if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96" or meth == "UA" or
meth == "LY04")):
cskinIn = 0
-elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
+elif ((cskinIn == None) and (meth == "C30" or meth == "C35"
or meth == "ecmwf" or meth == "Beljaars")):
cskinIn = 1
-if (np.all(gustIn == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
+if (np.all(gustIn == None) and (meth == "C30" or meth == "C35")):
gustIn = [1, 1.2, 600]
elif (np.all(gustIn == None) and (meth == "UA" or meth == "ecmwf")):
gustIn = [1, 1, 1000]
--
GitLab