diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index 19183943da1d667df7a962574347ea02edc455e4..038ed20517f2487b01b4538c9404d30b9125b0e7 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -71,10 +71,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
set 1 to keep points
L : int
Monin-Obukhov length definition options
- 0 : default for S80, S88, LP82, YT96 and LY04
- 1 : following UA (Zeng et al., 1998), default for UA
- 2 : following ERA5 (IFS Documentation cy46r1), default for ERA5
- 3 : COARE3.5 (Edson et al., 2013), default for C30, C35 and C40
+ "S80" : default for S80, S88, LP82, YT96 and LY04
+ "ERA5" : following ERA5 (IFS Documentation cy46r1), default for ERA5
Returns
-------
res : array that contains
@@ -113,6 +111,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
"""
logging.basicConfig(filename='flux_calc.log',
format='%(asctime)s %(message)s',level=logging.INFO)
+ logging.captureWarnings(True)
+ # check input values and set defaults where appropriate
lat, P, Rl, Rs, cskin, gust, tol, L = get_init(spd, T, SST, lat, P, Rl, Rs,
cskin, gust, L, tol, meth,
qmeth)
@@ -123,7 +123,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
' Rs: %s | gust: %s | cskin: %s | L : %s', meth,
np.nanmedian(lat), np.nanmedian(P), np.nanmedian(Rl),
np.nanmedian(Rs), gust, cskin, L)
- ####
+ # set up/calculate temperatures and specific humidities
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
@@ -135,9 +135,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
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")
- # first guesses
dt = Ta - sst
dq = qair - qsea
+ # first guesses
t10n, q10n = np.copy(Ta), np.copy(qair)
tv10n = t10n*(1 + 0.61*q10n)
# Zeng et al. 1998
@@ -162,6 +162,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
Rnl = 0.97*(5.67e-8*np.power(sst-0.3*cskin+CtoK, 4)-Rl)
dter = np.ones(T.shape)*0.3
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
+ # gustiness adjustment
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)))
@@ -169,6 +170,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
elif (gust[0] == 0):
wind = np.copy(spd)
+ # stars and roughness lengths
usr = np.sqrt(cd*np.power(wind, 2))
zo = 0.0001*np.ones(spd.shape)
zot, zoq = 0.0001*np.ones(spd.shape), 0.0001*np.ones(spd.shape)
@@ -177,11 +179,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
psit_calc(h_in[1]/monob, meth))
qsr = (dq+dqer*cskin)*kappa/(np.log(h_in[2]/zoq) -
psit_calc(h_in[2]/monob, meth))
+ # set-up to feed into iteration loop
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)
+ # iteration loop
while np.any(ii):
it += 1
if it > n:
@@ -288,14 +292,14 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
d = np.fabs(new-old)
if (tol[0] == 'flux'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
- (d[2, :] > tol[3]))
+ (d[2, :] > tol[3]))
elif (tol[0] == 'ref'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
- (d[2, :] > tol[3]))
+ (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]))
+ (d[2, :] > tol[3])+(d[3, :] > tol[4]) +
+ (d[4, :] > tol[5])+(d[5, :] > tol[6]))
if (ind[0].size == 0):
ii = False
else:
diff --git a/flux_subs.py b/flux_subs.py
index 1e17d5be14b092284b45059326d46f776079e45c..8e1c57c24dfca53730c491ac07c44f0b708aa240 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -651,10 +651,8 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
----------
L : int
Monin-Obukhov length definition options
- 0 : default for S80, S88, LP82, YT96 and LY04
- 1 : following UA (Zeng et al., 1998), default for UA
- 2 : following ERA5 (IFS Documentation cy46r1), default for ERA5
- 3 : COARE3.5 (Edson et al., 2013), default for C30, C35 and C40
+ "S80" : default for S80, S88, LP82, YT96 and LY04
+ "ERA5" : following ERA5 (IFS Documentation cy46r1), default for ERA5
lat : float
latitude
usr : float
@@ -702,19 +700,11 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
"""
g = gc(lat)
- if (L == 0):
+ if (L == "S80"):
tsrv = tsr+0.61*t10n*qsr
monob = ((tv10n*np.power(usr, 2))/(g*kappa*tsrv))
monob = np.where(np.fabs(monob) < 1, np.where(monob < 0, -1, 1), monob)
- elif (L == 1):
- Rb = g*h_in[0]*dtv/(tv*np.power(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))
- monob = h_in[0]/zol
- tsrv = tsr*(1.+0.61*qair)+0.61*th*qsr
- # monob = ((tv*np.power(usr, 2))/(kappa*g*tsrv))
- elif (L == 2):
+ elif (L == "ERA5"):
tsrv = tsr+0.61*t10n*qsr
Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_in[0])+0.61*dq)) /
np.power(wind, 2))
@@ -727,11 +717,6 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
monob = h_in[0]/zol
- elif (L == 3):
- tsrv = tsr+0.61*(T+CtoK)*qsr
- zol = (kappa*g*h_in[0]/(T+CtoK)*(tsr+0.61*(T+CtoK)*qsr) /
- np.power(usr, 2))
- monob = h_in[0]/zol
return tsrv, monob
#------------------------------------------------------------------------------
@@ -784,48 +769,50 @@ def get_strs(h_in, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, ct, cq,
"""
if (meth == "UA"):
- usr = np.where(h_in[0]/monob < -1.574, kappa*wind /
+ usr = np.where(h_in[0]/monob <= -1.574, kappa*wind /
(np.log(-1.574*monob/zo)-psim_calc(-1.574, meth) +
psim_calc(zo/monob, meth) +
1.14*(np.power(-h_in[0]/monob, 1/3) -
np.power(1.574, 1/3))),
- np.where((h_in[0]/monob > -1.574) & (h_in[0]/monob < 0),
- kappa*wind/(np.log(h_in[0]/zo) -
- psim_calc(h_in[0]/monob, meth) +
- psim_calc(zo/monob, meth)),
- np.where((h_in[0]/monob > 0) &
- (h_in[0]/monob < 1),
- kappa*wind/(np.log(h_in[0]/zo) +
- 5*h_in[0]/monob-5*zo/monob),
- kappa*wind/(np.log(monob/zo)+5-5*zo/monob +
- 5*np.log(h_in[0]/monob)+h_in[0]/monob-1))))
+ np.where(h_in[0]/monob < 0, kappa*wind /
+ (np.log(h_in[0]/zo) -
+ psim_calc(h_in[0]/monob, meth) +
+ psim_calc(zo/monob, meth)),
+ np.where(h_in[0]/monob <= 1, kappa*wind /
+ (np.log(h_in[0]/zo) +
+ 5*h_in[0]/monob-5*zo/monob),
+ kappa*wind/(np.log(monob/zo)+5 -
+ 5*zo/monob +
+ 5*np.log(h_in[0]/monob) +
+ h_in[0]/monob-1))))
# Zeng et al. 1998 (7-10)
tsr = np.where(h_in[1]/monob < -0.465, kappa*(dt+dter*cskin) /
(np.log((-0.465*monob)/zot) -
psit_calc(-0.465, meth)+0.8*(np.power(0.465, -1/3) -
np.power(-h_in[1]/monob, -1/3))),
- np.where((h_in[1]/monob > -0.465) & (h_in[1]/monob < 0),
- kappa*(dt+dter*cskin)/(np.log(h_in[1]/zot) -
- psit_calc(h_in[1]/monob, meth) +
- psit_calc(zot/monob, meth)),
- np.where((h_in[1]/monob > 0) & (h_in[1]/monob < 1),
- kappa*(dt+dter*cskin)/(np.log(h_in[1]/zot) +
- 5*h_in[1]/monob-5*zot/monob),
- kappa*(dt+dter*cskin)/(np.log(monob/zot)+5 -
- 5*zot/monob+5*np.log(h_in[1]/monob) +
- h_in[1]/monob-1))))
+ np.where(h_in[1]/monob < 0, kappa*(dt+dter*cskin) /
+ (np.log(h_in[1]/zot) -
+ psit_calc(h_in[1]/monob, meth) +
+ psit_calc(zot/monob, meth)),
+ np.where(h_in[1]/monob <= 1,
+ kappa*(dt+dter*cskin) /
+ (np.log(h_in[1]/zot) +
+ 5*h_in[1]/monob-5*zot/monob),
+ kappa*(dt+dter*cskin) /
+ (np.log(monob/zot)+5 -
+ 5*zot/monob+5*np.log(h_in[1]/monob) +
+ h_in[1]/monob-1))))
# Zeng et al. 1998 (11-14)
qsr = np.where(h_in[2]/monob < -0.465, kappa*(dq+dqer*cskin) /
(np.log((-0.465*monob)/zoq) -
psit_calc(-0.465, meth)+psit_calc(zoq/monob, meth) +
0.8*(np.power(0.465, -1/3) -
np.power(-h_in[2]/monob, -1/3))),
- np.where((h_in[2]/monob > -0.465) & (h_in[2]/monob < 0),
+ np.where(h_in[2]/monob < 0,
kappa*(dq+dqer*cskin)/(np.log(h_in[1]/zot) -
psit_calc(h_in[2]/monob, meth) +
psit_calc(zoq/monob, meth)),
- np.where((h_in[2]/monob > 0) &
- (h_in[2]/monob<1),
+ np.where(h_in[2]/monob <= 1,
kappa*(dq+dqer*cskin) /
(np.log(h_in[1]/zoq)+5*h_in[2]/monob -
5*zoq/monob),
diff --git a/get_init.py b/get_init.py
index 26699e32c3e99f3dae92401b4217369572634173..f991f0caa2977b91dc197cf578f978b205c6d4a9 100644
--- a/get_init.py
+++ b/get_init.py
@@ -99,11 +99,11 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
if (np.all(Rs == None) or np.all(np.isnan(Rs))):
Rs = np.ones(spd.shape)*150 # set to default for COARE3.5
if ((cskin == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
- or meth == "YT96" or meth == "LY04")):
+ or meth == "YT96" or meth == "UA" or
+ meth == "LY04")):
cskin = 0
- elif ((cskin == None) and (meth == "UA" or meth == "C30"
- or meth == "C35" or meth == "C40"
- or meth == "ERA5")):
+ elif ((cskin == None) and (meth == "C30" or meth == "C35" or meth == "C40"
+ or meth == "ERA5")):
cskin = 1
if ((gust == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
gust = [1, 1.2, 600]
@@ -113,19 +113,17 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
gust = [1, 1.2, 800]
elif (np.size(gust) < 3):
sys.exit("gust input must be a 3x1 array")
- if (L not in [None, 0, 1, 2, 3]):
+ if (L not in [None, "S80", "ERA5"]):
sys.exit("L input must be either None, 0, 1, 2 or 3")
if ((L == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
- or meth == "YT96" or meth == "LY04")):
- L = 0
- elif ((L == None) and (meth == "UA")):
- L = 1
+ or meth == "YT96" or meth == "LY04" or
+ meth == "UA" or meth == "C30" or meth == "C35"
+ or meth == "C40")):
+ L = "S80"
elif ((L == None) and (meth == "ERA5")):
- L = 2
- elif ((L == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
- L = 3
+ L = "ERA5"
if (tol == None):
tol = ['flux', 0.01, 1, 1]
elif (tol[0] not in ['flux', 'ref', 'all']):
sys.exit("unknown tolerance input")
- return lat, P, Rl, Rs, cskin, gust, tol, L
\ No newline at end of file
+ return lat, P, Rl, Rs, cskin, gust, tol, L
diff --git a/hum_subs.py b/hum_subs.py
index cb6a1b546e75ac7a925c87310b7fe66118d9d9f5..8528b80c1a8eda0384a7bc910fe8e921c44313df 100644
--- a/hum_subs.py
+++ b/hum_subs.py
@@ -1,6 +1,5 @@
import numpy as np
import sys
-import logging
from util_subs import (CtoK)
def VaporPressure(temp, P, phase, meth):
@@ -52,9 +51,7 @@ def VaporPressure(temp, P, phase, meth):
Author
------
Ported to Python and modified by S. Biri from Holger Voemel's original
-"""
- logging.basicConfig(filename='VaporPressure.log',
- format='%(asctime)s %(message)s',level=logging.info)
+ """
Psat = np.zeros(temp.size)*np.nan
if (np.nanmin(temp) > 200): # if Ta in Kelvin convert to Celsius
@@ -64,21 +61,20 @@ def VaporPressure(temp, P, phase, meth):
# Calculate saturation pressure over liquid water
if (phase == 'liquid'):
if (meth == 'HylandWexler' or meth == ''):
- logging.info("Source Hyland, R. W. and A. Wexler, Formulations \
- for the Thermodynamic Properties of the saturated \
- Phases of H2O from 173.15K to 473.15K, \
- ASHRAE Trans, 89(2A), 500-519, 1983.")
+ """
+ Source Hyland, R. W. and A. Wexler, Formulations for the
+ Thermodynamic Properties of the saturated Phases of H2O from
+ 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983."""
Psat = np.exp(-0.58002206e4/T+0.13914993e1-0.48640239e-1*T +
0.41764768e-4*np.power(T, 2) -
0.14452093e-7*np.power(T, 3) +
0.65459673e1*np.log(T))/100
if (meth == 'Hardy'):
- logging.info("Source Hardy, B., 1998, ITS-90 Formulations \
- for Vapor Pressure, Frostpoint Temperature, \
- Dewpoint Temperature, and Enhancement Factors \
- in the Range -100 to +100°C, The Proceedings of \
- the Third International Symposium on Humidity & \
- Moisture, London, England")
+ """
+ Source Hardy, B., 1998, ITS-90 Formulations for Vapor Pressure,
+ Frostpoint Temperature, Dewpoint Temperature, and Enhancement
+ Factors in the Range -100 to +100°C, The Proceedings of the Third
+ International Symposium on Humidity & Moisture, London, England"""
Psat = np.exp(-2.8365744e3/np.power(T, 2)-6.028076559e3/T +
1.954263612e1-2.737830188e-2*T +
1.6261698e-5*np.power(T, 2) +
@@ -86,18 +82,16 @@ def VaporPressure(temp, P, phase, meth):
1.8680009e-13*np.power(T, 4) +
2.7150305*np.log(T))/100
if (meth == 'Preining'):
- logging.info("Source : Vehkamaeki, H., M. Kulmala, I. Napari, \
- K. E. J. Lehtinen, C.Timmreck, M. Noppel, and \
- A. Laaksonen (2002), J. Geophys. Res., 107, \
- doi:10.1029/2002JD002184.")
+ """Source : Vehkamaeki, H., M. Kulmala, I. Napari, K. E. J.
+ Lehtinen, C.Timmreck, M. Noppel, and A. Laaksonen (2002),
+ J. Geophys. Res., 107, doi:10.1029/2002JD002184."""
Psat = np.exp(-7235.424651/T+77.34491296+5.7113e-3*T -
8.2*np.log(T))/100
if (meth == 'Wexler'):
- logging.info("Wexler, A., Vapor Pressure Formulation for Water \
- in Range 0 to 100 C. A Revision, Journal of Research \
- of the National Bureau of Standards - A. Physics and \
- Chemistry, September - December 1976, Vol. 80A, \
- Nos.5 and 6, 775-785")
+ """Wexler, A., Vapor Pressure Formulation for Water in Range 0 to
+ 100 C. A Revision, Journal of Research of the National Bureau of
+ Standards - A. Physics and Chemistry, September - December 1976,
+ Vol. 80A, Nos.5 and 6, 775-785"""
Psat = np.exp(-0.29912729e4*np.power(T, -2) -
0.60170128e4*np.power(T, -1) +
0.1887643854e2-0.28354721e-1*T +
@@ -106,14 +100,11 @@ def VaporPressure(temp, P, phase, meth):
0.44412543e-12*np.power(T, 4) + # This line was corrected from '-' to '+' following the original citation. (HV 20140819). The change makes only negligible difference
2.858487*np.log(T))/100
if ((meth == 'GoffGratch') or (meth == 'MartiMauersberger')):
- # logging.info("Marti and Mauersberger don't have a vapor pressure \
- # curve over liquid. Using Goff Gratch instead; \
- # Goff Gratch formulation Source : Smithsonian \
- # Meteorological Tables, 5th edition, p. 350, 1984 \
- # From original source: Goff & Gratch (1946), p. 107) \
- # Goff Gratch formulation; Source : Smithsonian \
- # Meteorological Tables, 5th edition, p. 350, 1984 \
- # From original source: Goff and Gratch (1946), p. 107")
+ """Marti and Mauersberger don't have a vapor pressure curve over
+ liquid. Using Goff Gratch instead; Goff Gratch formulation Source :
+ Smithsonian Meteorological Tables, 5th edition, p. 350, 1984
+ From original source: Goff & Gratch (1946), p. 107)
+ """
Ts = 373.16 # steam point temperature in K
ews = 1013.246 # saturation pressure at steam point temperature
Psat = np.power(10, -7.90298*(Ts/T-1)+5.02808*np.log10(Ts/T) -
@@ -121,80 +112,72 @@ def VaporPressure(temp, P, phase, meth):
8.1328e-3*(np.power(10, -3.49149*(Ts/T-1))-1) +
np.log10(ews))
if (meth == 'CIMO'):
- logging.info("Source: Annex 4B, Guide to Meteorological \
- Instruments and Methods of Observation, \
- WMO Publication No 8, 7th edition, Geneva, 2008. \
- (CIMO Guide)")
+ """Source: Annex 4B, Guide to Meteorological Instruments and
+ Methods of Observation, WMO Publication No 8, 7th edition, Geneva,
+ 2008. (CIMO Guide)"""
Psat = (6.112*np.exp(17.62*temp/(243.12+temp)) *
(1.0016+3.15e-6*P-0.074/P))
if (meth == 'MagnusTetens'):
- logging.info("Source: Murray, F. W., On the computation of \
+ """Source: Murray, F. W., On the computation of \
saturation vapor pressure, J. Appl. Meteorol., \
- 6, 203-204, 1967.")
+ 6, 203-204, 1967."""
Psat = np.power(10, 7.5*(temp)/(temp+237.5)+0.7858)
# Murray quotes this as the original formula and
Psat = 6.1078*np.exp(17.269388*temp/(temp+237.3))
# this as the mathematical aquivalent in the form of base e.
if (meth == 'Buck'):
- logging.info("Bucks vapor pressure formulation based on \
- Tetens formula. Source: Buck, A. L., \
- New equations for computing vapor pressure and \
- enhancement factor, J. Appl. Meteorol., 20, \
- 1527-1532, 1981.")
+ """Bucks vapor pressure formulation based on Tetens formula.
+ Source: Buck, A. L., New equations for computing vapor pressure and
+ enhancement factor, J. Appl. Meteorol., 20, 1527-1532, 1981."""
Psat = (6.1121*np.exp(17.502*temp/(240.97+temp)) *
(1.0007+(3.46e-6*P)))
if (meth == 'Buck2'):
- logging.info("Bucks vapor pressure formulation based on \
- Tetens formula. Source: Buck Research, \
- Model CR-1A Hygrometer Operating Manual, Sep 2001")
+ """Bucks vapor pressure formulation based on Tetens formula.
+ Source: Buck Research, Model CR-1A Hygrometer Operating Manual,
+ Sep 2001"""
Psat = (6.1121*np.exp((18.678-(temp)/234.5)*(temp)/(257.14+temp)) *
(1+1e-4*(7.2+P*(0.0320)+5.9e-6*np.power(T, 2))))
if (meth == 'WMO'):
- logging.info("Intended WMO formulation, originally published by \
- Goff (1957) incorrectly referenced by WMO technical \
- regulations, WMO-NO 49, Vol I, General \
- Meteorological Standards and Recommended Practices, \
- App. A, Corrigendum Aug 2000. and incorrectly \
- referenced by WMO technical regulations, \
- WMO-NO 49, Vol I, General Meteorological Standards \
- and Recommended Practices, App. A, 1988.")
+ """Intended WMO formulation, originally published by Goff (1957)
+ incorrectly referenced by WMO technical regulations, WMO-NO 49,
+ Vol I, General Meteorological Standards and Recommended Practices,
+ App. A, Corrigendum Aug 2000. and incorrectly referenced by WMO
+ technical regulations, WMO-NO 49, Vol I, General Meteorological
+ Standards and Recommended Practices, App. A, 1988."""
Ts = 273.16 # triple point temperature in K
Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
0.42873e-3*(np.power(10, -4.76955*(1-Ts/T))-1) +
0.78614)
if (meth == 'WMO2000'):
- logging.info("WMO formulation, which is very similar to Goff \
- Gratch. Source : WMO technical regulations, \
- WMO-NO 49, Vol I, General Meteorological Standards \
- and Recommended Practices, App. A, \
- Corrigendum Aug 2000.")
+ """WMO formulation, which is very similar to Goff & Gratch.
+ Source : WMO technical regulations, WMO-NO 49, Vol I, General
+ Meteorological Standards and Recommended Practices, App. A,
+ Corrigendum Aug 2000."""
Ts = 273.16 # triple point temperature in K
Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
0.42873e-3*(np.power(10, -4.76955*(1.-Ts/T))-1) +
0.78614)
if (meth == 'Sonntag'):
- logging.info("Source: Sonntag, D., Advancements in the field of \
- hygrometry, Meteorol. Z., N. F., 3, 51-66, 1994.")
+ """Source: Sonntag, D., Advancements in the field of hygrometry,
+ Meteorol. Z., N. F., 3, 51-66, 1994."""
Psat = np.exp(-6096.9385*np.power(T, -1)+16.635794 -
2.711193e-2*T+1.673952e-5*np.power(T, 2) +
2.433502*np.log(T))#*(1.0016+P*3.15e-6-0.074/P)
if (meth == 'Bolton'):
- logging.info("Source: Bolton, D., The computation of equivalent \
- potential temperature, Monthly Weather Report, 108, \
- 1046-1053, 1980. equation (10)")
+ """Source: Bolton, D., The computation of equivalent potential
+ temperature, Monthly Weather Report, 108, 1046-1053, 1980.
+ equation (10)"""
Psat = 6.112*np.exp(17.67*temp/(temp+243.5))
if (meth == 'IAPWS'):
- logging.info("Source: Wagner W. and A. Pruss (2002), The IAPWS \
- formulation 1995 for the thermodynamic properties \
- of ordinary water substance for general and \
- scientific use, J. Phys. Chem. Ref. Data, 31(2), \
- 387-535. This is the 'official' formulation from \
- the International Association for the Properties of \
- Water and Steam The valid range of this formulation \
- is 273.16 <= T <= 647.096 K and is based on the \
- ITS90 temperature scale.")
+ """Source: Wagner W. and A. Pruss (2002), The IAPWS formulation
+ 1995 for the thermodynamic properties of ordinary water substance
+ for general and scientific use, J. Phys. Chem. Ref. Data, 31(2),
+ 387-535. This is the 'official' formulation from the International
+ Association for the Properties of Water and Steam The valid range
+ of this formulation is 273.16 <= T <= 647.096 K and is based on the
+ ITS90 temperature scale."""
Tc = 647.096 # K : Temperature at the critical point
Pc = 22.064e4 # hPa : Vapor pressure at the critical point
nu = (1-T/Tc)
@@ -204,132 +187,111 @@ def VaporPressure(temp, P, phase, meth):
a3*np.power(nu, 3)+a4*np.power(nu, 3.5) +
a5*np.power(nu, 4)+ a6*np.power(nu, 7.5))))
if (meth == 'MurphyKoop'):
- logging.info("Source : Murphy and Koop, Review of the vapour \
- pressure of ice and supercooled water for \
- atmospheric applications, Q. J. R. Meteorol. \
- Soc (2005), 131, pp. 1539-1565.")
+ """Source : Murphy and Koop, Review of the vapour pressure of ice
+ and supercooled water for atmospheric applications, Q. J. R.
+ Meteorol. Soc (2005), 131, pp. 1539-1565."""
Psat = np.exp(54.842763-6763.22/T-4.210*np.log(T)+0.000367*T +
np.tanh(0.0415*(T-218.8))*(53.878-1331.22/T -
9.44523*np.log(T)+0.014025*T))/100
# Calculate saturation pressure over ice ----------------------------------
elif (phase == 'ice'):
- logging.info("Default uses Goff Gratch over ice. There is little \
- ambiguity in the ice saturation curve. Goff Gratch \
- is widely used.")
+ """Default uses Goff Gratch over ice. There is little ambiguity in the
+ ice saturation curve. Goff Gratch is widely used."""
if (meth == 'MartiMauersberger'):
- logging.info("Source : Marti, J. and K Mauersberger, A survey and \
- new measurements of ice vapor pressure at \
- temperatures between 170 and 250 K, GRL 20, \
- 363-366, 1993.")
+ """Source : Marti, J. and K Mauersberger, A survey and new
+ measurements of ice vapor pressure at temperatures between 170 and
+ 250 K, GRL 20, 363-366, 1993."""
Psat = np.power(10, -2663.5/T+12.537)/100
if (meth == 'HylandWexler'):
- logging.info("Source Hyland, R. W. and A. Wexler, Formulations \
- for the Thermodynamic Properties of the saturated \
- Phases of H2O from 173.15K to 473.15K, ASHRAE Trans,\
- 89(2A), 500-519, 1983.")
+ """Source Hyland, R. W. and A. Wexler, Formulations for the
+ Thermodynamic Properties of the saturated Phases of H2O from
+ 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983."""
Psat = np.exp(-0.56745359e4/T+0.63925247e1-0.96778430e-2*T +
0.62215701e-6*np.power(T, 2) +
0.20747825e-8*np.power(T, 3) -
0.9484024e-12*np.power(T, 4) +
0.41635019e1*np.log(T))/100
if (meth == 'Wexler'):
- logging.info("Wexler, A., Vapor pressure formulation for ice, \
- Journal of Research of the National Bureau of \
- Standards-A. 81A, 5-20, 1977.")
+ """Wexler, A., Vapor pressure formulation for ice, Journal of
+ Research of the National Bureau of Standards-A. 81A, 5-20, 1977."""
Psat = np.exp(-0.58653696e4*np.power(T, -1)+0.22241033e2 +
0.13749042e-1*T-0.34031775e-4*np.power(T, 2) +
0.26967687e-7*np.power(T, 3) +
0.6918651*np.log(T))/100
if (meth == 'Hardy'):
- logging.info("Source Hardy, B., 1998, ITS-90 Formulations for \
- Vapor Pressure, Frostpoint Temperature, Dewpoint \
- Temperature, and Enhancement Factors in the Range \
- -100 to +100°C, The Proceedings of the Third \
- International Symposium on Humidity & Moisture, \
- London, England. These coefficients are updated to \
- ITS90 based on the work by Bob Hardy at Thunder \
- Scientific: http://www.thunderscientific.com/\
- tech_info/reflibrary/its90formulas.pdf \
- The difference to the older ITS68 coefficients used \
- by Wexler is academic.")
+ """Source Hardy, B., 1998, ITS-90 Formulations for Vapor Pressure,
+ Frostpoint Temperature, Dewpoint Temperature, and Enhancement
+ Factors in the Range -100 to +100°C, The Proceedings of the Third
+ International Symposium on Humidity & Moisture, London, England.
+ These coefficients are updated to ITS90 based on the work by
+ Bob Hardy at Thunder Scientific: http://www.thunderscientific.com/
+ tech_info/reflibrary/its90formulas.pdf """
Psat = np.exp(-0.58666426e4*np.power(T, -1)+0.2232870244e2 +
0.139387003e-1*T-0.34262402e-4*np.power(T, 2) +
0.27040955e-7*np.power(T, 3) +
0.67063522e-1*np.log(T))/100
if (meth == 'GoffGratch' or meth == '' or meth == 'IAPWS'):
- logging.info("IAPWS does not provide a vapor pressure formulation \
- over ice use Goff Gratch instead.\
- Source : Smithsonian Meteorological Tables, \
- 5th edition, p. 350, 1984")
+ """IAPWS does not provide a vapor pressure formulation over ice use
+ Goff Gratch instead.
+ Source : Smithsonian Meteorological Tables, 5th edition, p. 350,
+ 1984"""
ei0 = 6.1071 # mbar
T0 = 273.16 # triple point in K
Psat = np.power(10, -9.09718*(T0/T-1)-3.56654*np.log10(T0/T) +
0.876793*(1-T/T0)+np.log10(ei0))
if (meth == 'MagnusTetens'):
- logging.info("Source: Murray, F. W., On the computation of \
- saturation vapor pressure, J. Appl. Meteorol., 6, \
- 203-204, 1967.")
+ """Source: Murray, F. W., On the computation of saturation vapor
+ pressure, J. Appl. Meteorol., 6, 203-204, 1967."""
Psat = np.power(10, 9.5*temp/(265.5+temp)+0.7858)
# Murray quotes this as the original formula and
Psat = 6.1078*np.exp(21.8745584*(T-273.16)/(T-7.66))
# this as the mathematical aquivalent in the form of base e.
if (meth == 'Buck'):
- logging.info("Bucks vapor pressure formulation based on Tetens \
- formula. Source: Buck, A. L., New equations for \
- computing vapor pressure and enhancement factor, \
- J. Appl. Meteorol., 20, 1527-1532, 1981.")
+ """Bucks vapor pressure formulation based on Tetens formula.
+ Source: Buck, A. L., New equations for computing vapor pressure and
+ enhancement factor, J. Appl. Meteorol., 20, 1527-1532, 1981."""
Psat = (6.1115*np.exp(22.452*temp/(272.55+temp)) *
(1.0003+(4.18e-6*P)))
if (meth == 'Buck2'):
- logging.info("Bucks vapor pressure formulation based on Tetens \
- formula. Source: Buck Research, Model CR-1A \
- Hygrometer Operating Manual, Sep 2001")
+ """Bucks vapor pressure formulation based on Tetens formula.
+ Source: Buck Research, Model CR-1A Hygrometer Operating Manual,
+ Sep 2001"""
Psat = (6.1115*np.exp((23.036-temp/333.7)*temp/(279.82+temp)) *
(1+1e-4*(2.2+P*(0.0383+6.4e-6*np.power(T, 2)))))
if (meth == 'CIMO'):
- logging.info("Source: Annex 4B, Guide to Meteorological \
- Instruments and Methods of Observation, \
- WMO Publication No 8, 7th edition, Geneva, 2008. \
- (CIMO Guide)")
+ """Source: Annex 4B, Guide to Meteorological Instruments and
+ Methods of Observation, WMO Publication No 8, 7th edition, Geneva,
+ 2008. (CIMO Guide)"""
Psat = (6.112*np.exp(22.46*temp/(272.62+temp)) *
(1.0016+3.15e-6*P-0.074/P))
if (meth == 'WMO' or meth == 'WMO2000'):
- logging.info("There is no typo issue in the WMO formulations for \
- ice. WMO formulation, which is very similar to Goff \
- Gratch. Source : WMO technical regulations, \
- WMO-NO 49, Vol I, General Meteorological Standards \
- and Recommended Practices, Aug 2000, App. A.")
+ """There is no typo issue in the WMO formulations for ice.
+ WMO formulation, which is very similar to Goff & Gratch.
+ Source : WMO technical regulations, WMO-NO 49, Vol I, General
+ Meteorological Standards and Recommended Practices, Aug 2000,
+ App. A."""
T0 = 273.16 # triple point temperature in K
Psat = np.power(10, -9.09685*(T0/T-1)-3.56654*np.log10(T0/T) +
0.87682*(1-T/T0)+0.78614)
if (meth == 'Sonntag'):
- logging.info("Source: Sonntag, D., Advancements in the field of \
- hygrometry, Meteorol. Z., N. F., 3, 51-66, 1994.")
+ """Source: Sonntag, D., Advancements in the field of hygrometry,
+ Meteorol. Z., N. F., 3, 51-66, 1994."""
Psat = np.exp(-6024.5282*np.power(T, -1)+24.721994+1.0613868e-2*T -
1.3198825e-5*np.power(T, 2)-0.49382577*np.log(T))
if (meth == 'MurphyKoop'):
- logging.info("Source : Murphy and Koop, Review of the vapour \
- pressure of ice and supercooled water for \
- atmospheric applications, Q. J. R. Meteorol. Soc \
- (2005), 131, pp. 1539-1565.")
+ """Source : Murphy and Koop, Review of the vapour pressure of ice
+ and supercooled water for atmospheric applications, Q. J. R.
+ Meteorol. Soc (2005), 131, pp. 1539-1565."""
Psat = np.exp(9.550426-5723.265/T+3.53068*np.log(T) -
0.00728332*T)/100
- if (meth == 'McIDAS'):
- logging.info("Source : Unknown, Received from Xin Jin \
- <xjin@ssec.wisc.edu>")
- A0, A1, A2 = 0.7859063157, 0.3579242320, -0.1292820828
- A3, A4, A5 = 0.5937519208, 0.4482949133, 0.2176664827
- E = A0+temp*(A1+temp*(A2+temp*(A3+temp*(A4+temp*A5))))
- Psat = np.power(10, E)
s = np.where(temp > 0)
if (s.size[0] >= 1):
- logging.info("Independent of the formula used for ice, use \
- Hyland Wexler (water) for temperatures above freezing\
- (see above). Source Hyland, R. W. and A. Wexler, \
- Formulations for the Thermodynamic Properties of the \
- saturated Phases of H2O from 173.15K to 473.15K, \
- ASHRAE Trans, 89(2A), 500-519, 1983.")
+ """Independent of the formula used for ice, use Hyland Wexler
+ (water) for temperatures above freezing (see above).
+ Source Hyland, R. W. and A. Wexler, Formulations for the
+ Thermodynamic Properties of the saturated Phases of H2O from
+ 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983."""
Psat_w = np.exp(-0.58002206e4/T+0.13914993e1-0.48640239e-1*T +
0.41764768e-4*np.power(T, 2) -
0.14452093e-7*np.power(T, 3) +
@@ -450,4 +412,4 @@ def get_hum(hum, T, sst, P, qmeth):
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)
return qair, qsea
-#------------------------------------------------------------------------------
\ No newline at end of file
+#------------------------------------------------------------------------------