Skip to content

GitLab

Commit b35e0d99 authored by sbiri's avatar sbiri
Browse files

get_init checks input values and sets defaults 

hum_subs contains all subroutines related to humidity calculations

util_subs contains all utility subroutines
parent 36e23fff
Hide whitespace changes
Inline Side-by-side
Showing with 680 additions and 0 deletions
get_init.py 0 → 100644
View file @ b35e0d99
import numpy as np
import sys
def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
"""
Checks initial input values and sets defaults if needed
Parameters
----------
spd : float
relative wind speed in m/s (is assumed as magnitude difference
between wind and surface current vectors)
T : float
air temperature in K
SST : float
sea surface temperature in K
lat : float
latitude (deg), default 45deg
P : float
air pressure (hPa), default 1013hPa
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
cskin : int
0 switch cool skin adjustment off, else 1
default is 1
gust : int
3x1 [x, beta, zi] x=1 to include the effect of gustiness, else 0
beta gustiness parameter, beta=1 for UA, beta=1.2 for COARE
zi PBL height (m) 600 for COARE, 1000 for UA and ERA5, 800 default
default for COARE [1, 1.2, 600]
default for UA, ERA5 [1, 1, 1000]
default else [1, 1.2, 800]
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
tol : float
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, 0.01, 1, 1]
meth : str
"S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
qmeth : str
is the saturation evaporation method to use amongst
"HylandWexler","Hardy","Preining","Wexler","GoffGratch","CIMO",
"MagnusTetens","Buck","Buck2","WMO","WMO2000","Sonntag","Bolton",
"IAPWS","MurphyKoop"]
default is Buck2
Returns
-------
lat : float
latitude
P : float
air pressure (hPa)
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
cskin : int
cool skin adjustment switch
gust : int
gustiness switch
tol : float
tolerance limits
L : int
MO length switch
"""
if ((type(spd) != np.ndarray) or (type(T) != np.ndarray) or
(type(SST) != np.ndarray)):
sys.exit("input type of spd, T and SST should be numpy.ndarray")
# if input values are nan break
if meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
"C40","ERA5"]:
sys.exit("unknown method")
if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler", "CIMO",
"GoffGratch", "MagnusTetens", "Buck", "Buck2", "WMO",
"WMO2000", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
sys.exit("unknown q-method")
if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))):
sys.exit("input wind, T or SST is empty")
if (np.all(lat == None)): # set latitude to 45deg if empty
lat = 45*np.ones(spd.shape)
elif ((np.all(lat != None)) and (np.size(lat) == 1)):
lat = np.ones(spd.shape)*np.copy(lat)
if ((np.all(P == None)) or np.all(np.isnan(P))):
P = np.ones(spd.shape)*1013
elif (((np.all(P != None)) or np.all(~np.isnan(P))) and np.size(P) == 1):
P = np.ones(spd.shape)*np.copy(P)
if (np.all(Rl == None) or np.all(np.isnan(Rl))):
Rl = np.ones(spd.shape)*370 # set to default for COARE3.5
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")):
cskin = 0
elif ((cskin == None) and (meth == "UA" or 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]
elif ((gust == None) and (meth == "UA" or meth == "ERA5")):
gust = [1, 1, 1000]
elif (gust == None):
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]):
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
elif ((L == None) and (meth == "ERA5")):
L = 2
elif ((L == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
L = 3
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
hum_subs.py 0 → 100644
View file @ b35e0d99
import numpy as np
import sys
import logging
from util_subs import (CtoK)
def VaporPressure(temp, P, phase, meth):
"""
Calculate the saturation vapor pressure
For temperatures above 0 deg C the vapor pressure over liquid water
is calculated.
The optional parameter 'liquid' changes the calculation to vapor pressure
over liquid water over the entire temperature range.
The current default fomulas are Hyland and Wexler for liquid and
Goff Gratch for ice.
Parameters
----------
temp : float
Temperature in [deg C]
phase : str
'liquid' : Calculate vapor pressure over liqiud water or
'ice' : Calculate vapor pressure over ice
meth : str
formula to be used
MartiMauersberger : vaporpressure formula from Marti Mauersberger
MagnusTetens : vaporpressure formula from Magnus Tetens
GoffGratch : vaporpressure formula from Goff Gratch
Buck_original : vaporpressure formula from Buck (original)
Buck_manual : vaporpressure formula from the Buck manual
CIMO : vaporpressure formula recommended by CIMO
Vaisala : vaporpressure formula from Vaisala
WMO_Goff : vaporpressure formula from WMO, which should have been Goff
WMO2000 : vaporpressure formula from WMO (2000) containing a typo
Wexler : vaporpressure formula from Wexler (1976)
Sonntag : vaporpressure formula from Sonntag (1994)
Bolton : vaporpressure formula from Bolton (1980)
Fukuta : vaporpressure formula from Fukuta (2003)
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)
Returns
-------
P : float
Saturation vapor pressure in [hPa]
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
temp = temp-273.16
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 == ''):
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.")
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")
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) +
7.0229056e-10*np.power(T, 3) -
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.")
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")
Psat = np.exp(-0.29912729e4*np.power(T, -2) -
0.60170128e4*np.power(T, -1) +
0.1887643854e2-0.28354721e-1*T +
0.17838301e-4*np.power(T, 2) -
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')):
# 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")
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) -
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 == 'CIMO'):
logging.info("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 \
saturation vapor pressure, J. Appl. Meteorol., \
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.")
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")
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.")
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.")
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.")
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)")
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.")
Tc = 647.096 # K : Temperature at the critical point
Pc = 22.064e4 # hPa : Vapor pressure at the critical point
nu = (1-T/Tc)
a1, a2, a3 = -7.85951783, 1.84408259, -11.7866497
a4, a5, a6 = 22.6807411, -15.9618719, 1.80122502
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'):
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.")
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.")
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.")
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.")
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.")
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.")
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")
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.")
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.")
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")
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)")
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.")
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.")
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.")
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.")
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) +
0.65459673e1*np.log(T))/100
Psat[s] = Psat_w[s]
return Psat
#-----------------------------------------------------------------------------
def qsat_sea(T, P, qmeth):
""" Computes surface saturation specific humidity (g/kg)
Parameters
----------
T : float
temperature ($^\\circ$\\,C)
P : float
pressure (mb)
qmeth : str
method to calculate vapor pressure
Returns
-------
qs : float
"""
T = np.asarray(T)
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) as in C35
Parameters
----------
T : float
temperature ($^\circ$\,C)
P : float
pressure (mb)
rh : float
relative humidity (%)
qmeth : str
method to calculate vapor pressure
Returns
-------
q : float
em : float
"""
T = np.asarray(T)
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
Parameters
----------
hum : array
humidity input switch 2x1 [x, values] default is relative humidity
x='rh' : relative humidity in %
x='q' : specific humidity (g/kg)
x='Td' : dew point temperature (K)
T : float
air temperature in K
sst : float
sea surface temperature in K
P : float
air pressure at sea level in hPa
qmeth : str
method to calculate specific humidity from vapor pressure
Returns
-------
qair : float
specific humidity of air
qsea : float
specific humidity over sea surface
"""
if (hum == None):
RH = np.ones(sst.shape)*80
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg)
elif (hum[0] not in ['rh', 'q', 'Td']):
sys.exit("unknown humidity input")
qair, qsea = np.nan, np.nan
elif (hum[0] == 'rh'):
RH = hum[1]
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 (g/kg)
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg)
elif (hum[0] == 'q'):
qair = hum[1]
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
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)))
es = 611.21*np.exp(17.502*((T-273.16)/(T-32.19)))
RH = 100*esd/es
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
util_subs.py 0 → 100644
View file @ b35e0d99
import numpy as np
CtoK = 273.16 # 273.15
""" Conversion factor for $^\circ\,$C to K """
kappa = 0.4 # NOTE: 0.41
""" von Karman's constant """
#------------------------------------------------------------------------------
def get_heights(h, dim_len):
""" Reads input heights for velocity, temperature and humidity
Parameters
----------
h : float
input heights (m)
dim_len : int
length dimension
Returns
-------
hh : array
"""
hh = np.zeros((3, dim_len))
if (type(h) == float or type(h) == int):
hh[0, :], hh[1, :], hh[2, :] = h, h, h
elif (len(h) == 2 and np.ndim(h) == 1):
hh[0, :], hh[1, :], hh[2, :] = h[0], h[1], h[1]
elif (len(h) == 3 and np.ndim(h) == 1):
hh[0, :], hh[1, :], hh[2, :] = h[0], h[1], h[2]
elif (len(h) == 1 and np.ndim(h) == 2):
hh = np.zeros((3, h.shape[1]))
hh[0, :], hh[1, :], hh[2, :] = h[0, :], h[0, :], h[0, :]
elif (len(h) == 2 and np.ndim(h) == 2):
hh = np.zeros((3, h.shape[1]))
hh[0, :], hh[1, :], hh[2, :] = h[0, :], h[1, :], h[1, :]
elif (len(h) == 3 and np.ndim(h) == 2):
hh = np.zeros((3, h.shape[1]))
hh = np.copy(h)
return hh
# ---------------------------------------------------------------------
def gc(lat, lon=None):
""" Computes gravity relative to latitude
Parameters
----------
lat : float
latitude ($^\circ$)
lon : float
longitude ($^\circ$, optional)
Returns
-------
gc : float
gravity constant (m/s^2)
"""
gamma = 9.7803267715
c1 = 0.0052790414
c2 = 0.0000232718
c3 = 0.0000001262
c4 = 0.0000000007
if lon is not None:
lon_m, lat_m = np.meshgrid(lon, lat)
else:
lat_m = lat
phi = lat_m*np.pi/180.
xx = np.sin(phi)
gc = (gamma*(1+c1*np.power(xx, 2)+c2*np.power(xx, 4)+c3*np.power(xx, 6) +
c4*np.power(xx, 8)))
return gc
# ---------------------------------------------------------------------
def visc_air(T):
""" Computes the kinematic viscosity of dry air as a function of air temp.
following Andreas (1989), CRREL Report 89-11.
Parameters
----------
Ta : float
air temperature ($^\circ$\,C)
Returns
-------
visa : float
kinematic viscosity (m^2/s)
"""
T = np.asarray(T)
if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius
T = T-273.16
visa = 1.326e-5*(1+6.542e-3*T+8.301e-6*np.power(T, 2) -
4.84e-9*np.power(T, 3))
return visa
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