Skip to content

GitLab

Commit e3b57391 authored by sbiri's avatar sbiri
Browse files

Update Code/hum_subs.py

parent 7db6b68f
Hide whitespace changes
Inline Side-by-side
Showing with 67 additions and 67 deletions
Code/hum_subs.py
View file @ e3b57391
import numpy as np
import sys
import numpy as np
from util_subs import (CtoK)
def VaporPressure(temp, P, phase, meth):
"""
Calculate the saturation vapor pressure
Calculate the saturation vapor pressure.
For temperatures above 0 deg C the vapor pressure over liquid water
is calculated.
......@@ -34,7 +35,8 @@ def VaporPressure(temp, P, phase, meth):
Wexler : vaporpressure formula from Wexler (1976)
Sonntag : vaporpressure formula from Sonntag (1994)
Bolton : vaporpressure formula from Bolton (1980)
HylandWexler : vaporpressure formula from Hyland and Wexler (1983)
HylandWexler : vaporpressure formula from Hyland and Wexler
(1983)
IAPWS : vaporpressure formula from IAPWS (2002)
Preining : vaporpressure formula from Preining (2002)
MurphyKoop : vaporpressure formula from Murphy and Koop (2005)
......@@ -49,16 +51,14 @@ def VaporPressure(temp, P, phase, meth):
------
Ported to Python and modified by S. Biri from Holger Voemel's original
"""
Psat = np.zeros(temp.size)*np.nan
if (np.nanmin(temp) > 200): # if Ta in Kelvin convert to Celsius
if np.nanmin(temp) > 200: # if Ta in Kelvin convert to Celsius
temp = temp-273.16
T = np.copy(temp)+273.16 # Most formulas use T in [K]
T = np.copy(temp)+273.16 # Most formulas use T in [K]
# Formulas using [C] use the variable temp
# Calculate saturation pressure over liquid water
if (phase == 'liquid'):
if (meth == 'HylandWexler' or meth == ''):
if phase == 'liquid':
if meth in (['HylandWexler', '']):
"""
Source Hyland, R. W. and A. Wexler, Formulations for the
Thermodynamic Properties of the saturated Phases of H2O from
......@@ -67,7 +67,7 @@ def VaporPressure(temp, P, phase, meth):
0.41764768e-4*np.power(T, 2) -
0.14452093e-7*np.power(T, 3) +
0.65459673e1*np.log(T))/100
if (meth == 'Hardy'):
if meth == 'Hardy':
"""
Source Hardy, B., 1998, ITS-90 Formulations for Vapor Pressure,
Frostpoint Temperature, Dewpoint Temperature, and Enhancement
......@@ -79,13 +79,13 @@ def VaporPressure(temp, P, phase, meth):
7.0229056e-10*np.power(T, 3) -
1.8680009e-13*np.power(T, 4) +
2.7150305*np.log(T))/100
if (meth == 'Preining'):
if meth == 'Preining':
"""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'):
if meth == 'Wexler':
"""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,
......@@ -97,7 +97,7 @@ def VaporPressure(temp, P, phase, meth):
0.84150417e-9*np.power(T, 3) +
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')):
if meth in (['GoffGratch', 'MartiMauersberger']):
"""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
......@@ -109,7 +109,7 @@ def VaporPressure(temp, P, phase, meth):
1.3816e-7*(np.power(10, 11.344*(1-T/Ts))-1) +
8.1328e-3*(np.power(10, -3.49149*(Ts/T-1))-1) +
np.log10(ews))
if (meth == 'MagnusTetens'):
if meth == 'MagnusTetens':
"""Source: Murray, F. W., On the computation of \
saturation vapor pressure, J. Appl. Meteorol., \
6, 203-204, 1967."""
......@@ -117,19 +117,19 @@ def VaporPressure(temp, P, phase, meth):
# 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'):
if meth == 'Buck':
"""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'):
if meth == 'Buck2':
"""Bucks vapor pressure formulation based on Tetens formula.
Source: Buck Research, Model CR-1A Hygrometer Operating Manual,
May 2012"""
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'):
if meth == 'WMO':
"""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,
......@@ -141,22 +141,22 @@ def VaporPressure(temp, P, phase, meth):
1.50475e-4*(1-10**(-8.2969*(T/Ts-1))) +
0.42873e-3*(10**(4.76955*(1-Ts/T))-1) + # in eq. 13 is -4.76955; in aerobulk is like this
0.78614)
if (meth == 'WMO2018'):
if meth == 'WMO2018':
"""WMO 2018 edition. Annex 4.B, eq. 4.B.1, 4.B.2, 4.B.5 """
Psat = 6.112*np.exp(17.62*temp/(243.12+temp))*(1.0016+3.15e-6*P -
0.074/P)
if (meth == 'Sonntag'):
if meth == 'Sonntag':
"""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'):
2.433502*np.log(T)) # *(1.0016+P*3.15e-6-0.074/P)
if meth == 'Bolton':
"""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'):
if meth == 'IAPWS':
"""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),
......@@ -172,7 +172,7 @@ def VaporPressure(temp, P, phase, meth):
Psat = (Pc*np.exp(Tc/T*(a1*nu+a2*np.power(nu, 1.5) +
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'):
if meth == 'MurphyKoop':
"""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."""
......@@ -180,15 +180,15 @@ def VaporPressure(temp, P, phase, meth):
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'):
elif phase == 'ice':
"""Default uses Goff Gratch over ice. There is little ambiguity in the
ice saturation curve. Goff Gratch is widely used."""
if (meth == 'MartiMauersberger'):
if meth == 'MartiMauersberger':
"""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'):
if meth == 'HylandWexler':
"""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."""
......@@ -197,14 +197,14 @@ def VaporPressure(temp, P, phase, meth):
0.20747825e-8*np.power(T, 3) -
0.9484024e-12*np.power(T, 4) +
0.41635019e1*np.log(T))/100
if (meth == 'Wexler'):
if meth == 'Wexler':
"""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'):
if meth == 'Hardy':
"""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
......@@ -216,7 +216,7 @@ def VaporPressure(temp, P, phase, meth):
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'):
if meth in (['GoffGratch', '', 'IAPWS']):
"""IAPWS does not provide a vapor pressure formulation over ice use
Goff Gratch instead.
Source : Smithsonian Meteorological Tables, 5th edition, p. 350,
......@@ -225,32 +225,32 @@ def VaporPressure(temp, P, phase, meth):
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'):
if meth == 'MagnusTetens':
"""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'):
if meth == 'Buck':
"""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'):
if meth == 'Buck2':
"""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'):
if meth == 'CIMO':
"""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'):
if meth in ('WMO', 'WMO2000'):
"""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
......@@ -259,20 +259,20 @@ def VaporPressure(temp, P, phase, meth):
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'):
if meth == 'Sonntag':
"""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'):
if meth == 'MurphyKoop':
"""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
s = np.where(temp > 0)
if (s.size[0] >= 1):
# s = np.where(temp > 0)
if np.where(temp > 0).size[0] >= 1:
"""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
......@@ -282,20 +282,20 @@ def VaporPressure(temp, P, phase, meth):
0.41764768e-4*np.power(T, 2) -
0.14452093e-7*np.power(T, 3) +
0.65459673e1*np.log(T))/100
Psat[s] = Psat_w[s]
Psat[np.where(temp > 0)] = Psat_w[np.where(temp > 0)]
return Psat
#-----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
def qsat_sea(T, P, qmeth):
"""
computes surface saturation specific humidity [g/kg]
r"""
Compute surface saturation specific humidity [g/kg].
Parameters
----------
T : float
temperature ($^\\circ$\\,C)
temperature ($^\circ$\,C)
P : float
pressure (mb)
qmeth : str
......@@ -306,18 +306,18 @@ def qsat_sea(T, P, qmeth):
qs : float
"""
T = np.asarray(T)
if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius
if np.nanmin(T) > 200: # if Ta in Kelvin convert to Celsius
T = T-CtoK
ex = VaporPressure(T, P, 'liquid', qmeth)
es = 0.98*ex # reduction at sea surface
qs = 622*es/(P-0.378*es)
return qs
#-----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
def qsat_air(T, P, rh, qmeth):
"""
computes saturation specific humidity [g/kg]
r"""
Compute saturation specific humidity [g/kg].
Parameters
----------
......@@ -336,18 +336,18 @@ def qsat_air(T, P, rh, qmeth):
em : float
"""
T = np.asarray(T)
if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius
if np.nanmin(T) > 200: # if Ta in Kelvin convert to Celsius
T = T-CtoK
es = VaporPressure(T, P, 'liquid', qmeth)
em = 0.01*rh*es
q = 622*em/(P-0.378*em)
return q
#------------------------------------------------------------------------------
# -----------------------------------------------------------------------------
def get_hum(hum, T, sst, P, qmeth):
"""
Get specific humidity output
Get specific humidity output.
Parameters
----------
......@@ -378,16 +378,16 @@ def get_hum(hum, T, sst, P, qmeth):
qair, qsea = np.nan, np.nan
elif ((hum[0] == 'rh') or (hum[0] == 'no')):
RH = hum[1]
if (np.all(RH < 1)):
sys.exit("input relative humidity units should be \%")
if np.all(RH < 1):
sys.exit("input relative humidity units should be %")
qair, qsea = np.nan, np.nan
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (kg/kg)
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (kg/kg)
elif (hum[0] == 'q'):
elif hum[0] == 'q':
qair = hum[1]
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (kg/kg)
elif (hum[0] == 'Td'):
Td = hum[1] # dew point temperature (K)
elif hum[0] == 'Td':
Td = hum[1] # dew point temperature (K)
Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
T = np.where(T < 200, np.copy(T)+CtoK, np.copy(T))
esd = 611.21*np.exp(17.502*((Td-273.16)/(Td-32.19)))
......@@ -396,12 +396,12 @@ def get_hum(hum, T, sst, P, qmeth):
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (kg/kg)
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (kg/kg)
return qair, qsea
#------------------------------------------------------------------------------
# -----------------------------------------------------------------------------
def gamma(opt, sst, t, q, cp):
"""
Computes the adiabatic lapse-rate
Compute the adiabatic lapse-rate.
Parameters
----------
......@@ -424,26 +424,26 @@ def gamma(opt, sst, t, q, cp):
lapse rate [K/m]
"""
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
if np.nanmin(sst) < 200: # if sst in Celsius convert to Kelvin
sst = sst+CtoK
if (np.nanmin(t) < 200): # if sst in Celsius convert to Kelvin
if np.nanmin(t) < 200: # if t in Celsius convert to Kelvin
t = t+CtoK
if (opt == "moist"):
if opt == "moist":
t = np.maximum(t, 180)
q = np.maximum(q, 1e-6)
w = q/(1-q) # mixing ratio w = q/(1-q)
w = q/(1-q) # mixing ratio w = q/(1-q)
iRT = 1/(287.05*t)
# latent heat of vaporization of water as a function of temperature
lv = (2.501-0.00237*(sst-CtoK))*1e6
gamma = 9.8*(1+lv*w*iRT)/(1005+np.power(lv, 2)*w*(287.05/461.495) *
iRT/t)
elif (opt == "dry_c"):
elif opt == "dry_c":
gamma = 0.0098*np.ones(t.shape)
elif (opt == "dry"):
elif opt == "dry":
gamma = 9.81/cp
elif (opt == "dry_v"):
w = q/(1-q) # mixing ratio
f_v = 1-0.85*w #(1+w)/(1+w*)
elif opt == "dry_v":
w = q/(1-q) # mixing ratio
f_v = 1-0.85*w # (1+w)/(1+w*)
gamma = f_v*9.81/cp
return gamma
#------------------------------------------------------------------------------
# -----------------------------------------------------------------------------
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