#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""**Meteorological constants and functions**
Based on earlier versions of the module dating back to 2009 named admeteopy.
© Dietmar Thaler 2009-2019
| E-Mail: info@foehnwall.at
| Web : https://python3.foehnwall.at/
.. note:: This module makes use of the **numpy** library. NumPy is the
fundamental package for scientific computing with Python. Make sure
you have it installed. See http://www.numpy.org/ and
http://sourceforge.net/projects/numpy/
------------------------------------------------------------
*This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.*
*This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.*
*You should have received a copy of the GNU General Public License.
If not, see* <http://www.gnu.org/licenses/gpl.html> .
-----------------------------------------------------------
"""
import csv
import sys
import matplotlib.pyplot as plt
import numpy as np
__version__ = "0.3.1-20190809" # : module version
# __version__ = "0.3.0-20180508" # : module version
# __version__ = "0.1.0-20180211" # : module version
__copyright__ = "D. Thaler"
M0 = 28.9644 # : mean relative molecular mass of dry air (kg/kmol)
R = 8314.32 # : J/(K kmol)
Rl = R / M0 # : Ideal gas constant of dry air in J/kg/K
cp_l = 1004.0 # : cp for dry air [J/kg/K]
cv_l = cp_l - Rl # : valid for an ideal gas
g0 = 9.80665 # : Standard gravity (concerning to ICAO)
Tice = 273.15 # : Freezing point of water in K
gammastand = 0.0065 # : standard temperature gradient (ICAO)
sigma = 5.670e-8 # : W/m^2/K^4
k_Boltzmann = 1.38064852e-23 # : Boltzmann KonstanteJ K-1
c_light = 299792458. # : m/s speed of light
h_Planck = 6.626070040e-34 # : Planck constant J s
_hc_k = h_Planck * c_light / k_Boltzmann # : hc/k
_h_k = h_Planck / k_Boltzmann # : h/k
NAUTICAL_MILE = 1852.0 # m
KT2MS = NAUTICAL_MILE / 3600.
KT2KMH = NAUTICAL_MILE / 1000.
def __makearray__(obj):
"""Converts a scalar or a list etc ... to a numpy array and returns
a numpy array unchanged"""
if np.isscalar(obj):
return np.array([obj])
else:
return np.array(obj)
[docs]def kelvin2celsius(tk):
"""Converts temperature in Kelvin to temperature in Ceslius"""
return tk - Tice
[docs]def c2k(tk):
"""Converts temperature in Celsius to temperature in Kelvin"""
return celsius2kelvin(tk)
[docs]def celsius2kelvin(tc):
"""Converts temperature in Celsius to temperature in Kelvin"""
return tc + Tice
[docs]def k2c(tk):
"""Converts temperature in Kelvin to temperature in Ceslius"""
return kelvin2celsius(tk)
[docs]def g_welmec(phi=45.0, z=0.0):
"""
WELMEC-formula for the caculation of gravity as function
of latitude and height
| phi ... latitude in decimal grade (positive north)
| z ... height in meter
| g_welmec... gravity in m/s^2
"""
p = phi / 180.0 * np.pi
# print phi,p
g = 9.780318 * (1.0 + 0.0053024 * np.sin(p) * np.sin(p) -
0.0000058 * np.sin(2 * p) * np.sin(2 * p)) - \
0.000003085 * z
return g
[docs]def gpotheight(phi=45.0, z=0.0):
"""
Geopotential (ICAO-) height according to the WELMEC-formula
for the caculation of gravity as function of latitude and height.
gp(z) = integral(g(z)dz) from ground with z=0 to height z.
| phi ... latitude in decimal grade (positive north)
| z ... height in meter
| gpotheight ... geopotential in units of g0 (ICAO)
"""
p = phi / 180.0 * np.pi
gp = 9.780318 * (1.0 + 0.0053024 * np.sin(p) * np.sin(p) -
0.0000058 * np.sin(2 * p) * np.sin(2 * p)) * z\
-0.5 * 0.000003085 * z * z
return gp / g0
[docs]def p_iso(p0, T0, z, g=g0):
""" Isothermal atmosphere - pressure reduction
| p0 .. pressure at level 0
| T0 .. temperature in K
| z .. thickness level in [m]
| g .. gravity, defaults to g0
| p_iso .. pressure at level z
"""
return p0 * np.exp(-g * z / Rl / T0)
[docs]def z_iso(T0, p0, p1, g=g0):
"""
Isothermal atmosphere - thickness
| p0 .. pressure at level 0
| p1 .. pressure at level 1
| T0 .. temperature in K
| g .. gravity, defaults to standard gravity g0
| z_iso .. thickness between level 1 an 0 in [m]
"""
return Rl * T0 / g * np.log(p0 / p1)
[docs]def p_poly(p0, T0, z, gamma, g=g0):
"""
Polytropic atmosphere - pressure reduction
| p0 .. pressure at level 0
| T0 .. temperature in K
| z .. thickness level in [m]
| gamma = -dT/dz .. vertical temperature gradient [K/m]
| g .. gravity, defaults to standard gravity g0
| p_poly .. pressure at level z
"""
return p0 * ((1.0 - z / T0 * gamma) ** (g / gamma / Rl))
[docs]def z_poly(T0, p0, p1, gamma, g=g0):
""" Polytropic atmosphere - thickness
| p0 .. pressure at level 0
| p1 .. pressure at level 1
| T0 .. temperature in K
| gamma=-dT/dz .. vertical temperature gradient [K/m]
| g .. gravity, defaults to standard gravity g0
| z_poly .. thickness between level 1 an 0 in [m]
"""
return T0 / gamma * (1 - (p0 / p1) ** (-Rl / g * gamma))
[docs]def z_poly2(T0, p0, T1, p1, g=g0):
"""
Polytropic atmosphere - thickness
| p0 .. pressure at level 0
| p1 .. pressure at level 1
| T0 .. temperature in K in level 0
| T1 .. temperature in K in level 1
| g .. gravity, defaults to standard gravity g0
| z_poly2 .. thickness between level 1 an 0 in [gpm]
"""
return (T0 - T1) * Rl / g * np.log(p1 / p0) / np.log(T1 / T0)
[docs]def p0_wmocimo2008(ps, Hp, ts, es):
"""
Pressure reduction to mean sealevel according to
WMO CIMO Guide, Part I, Chapter 3 (Edition 2008, updated in 2010)
http://www.wmo.int/pages/prog/www/IMOP/CIMO-Guide.html
| ps ... station pressure
| Hp ... station height in geopotential meter (local gravity correction)
| ts ... station temperature in C
| es ... station vapor presser in hPa
"""
a = 0.0065 # K/gpm
Ch = 0.12 # K
gn = 9.80665 # m/s^2
R = 287.05 # J/(K kg)
Ts = ts + 273.15 # K
return ps * np.exp((gn / R * Hp) / (Ts + a * Hp / 2 + es * Ch))
def _p_corr_factor_e_(p):
"""Pressur correction for water vapor correction f(p) according to
ANNEX 4.B. FORMULAE FOR THE COMPUTATION OF MEASURES OF HUMIDITY
p. 163 (Guide to Meteorological Instruments and Methods of Observation,
WMO, Geneva, 2014, 2017). - e_water*f(p), e_ice*f(p)"""
f = 1.0016 + 3.15e-6 * p - 0.074 / p
return f
[docs]def e_water(t=0.0):
"""
Saturation water vapor over a plane liqid water surface
(according http://cires.colorado.edu/~voemel/vp.html
or http://cires.colorado.edu/~voemel/vp.html
Guide to Meteorological Instruments and Methods of Observation,
CIMO Guide, WMO 2008)
| t ... temperature t in C
| e_water ... sat. water vapor in hPa
"""
E0 = 6.112 # hPa
E = E0 * np.exp(17.62 * t / (243.12 + t))
return E
[docs]def e_water_p(t=0.0, p=1013.25):
"""Saturation water vapor over water with pressure correction.
ANNEX 4.B. FORMULAE FOR THE COMPUTATION OF MEASURES OF HUMIDITY
p. 163 (Guide to Meteorological Instruments and Methods of Observation,
WMO, Geneva, 2014, 2017). - e_water*f(p)
| t ... temperature t in C
| p ... pressure in hPa
| e_water_p ... sat. water vapor in hPa"""
return e_water(t) * _p_corr_factor_e_(p)
[docs]def e_ice(t=0.0):
"""
Saturation water vapor over a plane frozen water (ice) surface
(according http://cires.colorado.edu/~voemel/vp.html
or http://cires.colorado.edu/~voemel/vp.html
"Guide to Meteorological Instruments and Methods of Observation",
CIMO Guide, WMO 2008)
| t ... temperature t in C
| e_ice ... sat. water vapor in hPa
"""
E0 = 6.112 # hPa
E = E0 * np.exp(22.46 * t / (272.62 + t))
return E
[docs]def e_ice_p(t=0.0, p=1013.25):
"""Saturation water vapor over ice with pressure correction.
ANNEX 4.B. FORMULAE FOR THE COMPUTATION OF MEASURES OF HUMIDITY
p. 163 (Guide to Meteorological Instruments and Methods of Observation,
WMO, Geneva, 2014, 2017). - e_ice*f(p)
| t ... temperature t in C
| p ... pressure in hPa
| e_ice_p ... sat. water vapor in hPa"""
return e_ice(t) * _p_corr_factor_e_(p)
[docs]def l_evw(t=0.0):
"""
Latent heat of evaporation/condensation of water (vapor - liquid)
valid from -40 to +40 C.
Cubic fit to Table 2.1,p.16, Textbook: R.R.Rogers & M.K. Yau,
A Short Course in Cloud Physics, 3e,(1989), Pergamon press
http://en.wikipedia.org/wiki/Latent_heat (2009-11-09)
V 2009-11-09, (p) dietmar.thaler@gmx.at
| T .. Temperature in C
| l_evw .. Latent heat in J/kg as function of Temperature
"""
L = (-0.0000614342 * t ** 3 + 0.00158927 * t **
2 - 2.36418 * t + 2500.79) * 1000.0
return L
[docs]def q_spechum(e, p):
"""
Specific humidity as function of water-vapor pressure and air pressure
| e ... vapor pressure
| p ... air pressure
| q_spechum ... specific humidity in kg/kg
"""
q = 0.622 * e / (p - 0.378 * e)
return q
[docs]def m_mixingratio(e, p):
"""
mixing ratio as function of water-vapor pressure and air pressure
| e ... vapor pressur
| p ... air pressure
| m_mixingratio ... mixing ratio in kg/kg
"""
m = 0.622 * e / (p - e) # Bolton (1980)
return m
[docs]def sat_mixingratio_water(T, p):
"""
Saturation mixing ratio as function of temperature and air pressure
| T ... Temperature in C
| p ... Pressure in hPa
| sat_mixingratio_water .. Saturation mixing ratio
"""
e = e_water(T)
m = m_mixingratio(e, p)
return m
[docs]def t_virt(T, q):
"""
Virtual Temperatur as function of air temperature
and spezific humidity
| T ... air temperature [K]
| q ... spezific humidity [kg/kg]
| t_virt ... virtual temperature [K]
"""
return T * (1.0 + 0.608 * q)
[docs]def t_virt2(t, td, p):
"""Virtual Temperature as function of temp., dewpoint and
pressure
| t ... temp. in C(!)
| td ... dewpoint temp. in C(!)
| p ... air pressure in hPa
| t_virt2 ... virtual temp. in K(!)
"""
T = t + Tice # K
ew = e_water(td)
q = q_spechum(ew, p)
Tv = t_virt(T, q)
return Tv
[docs]def pot_temp(T, p):
"""
Potential temperatur as function of temperature and pressure
(reduction to ps=1000.0 hpa)
| T .. temperatur of dry air [K]
| p .. pressure of air [hPa !!!]
| pot_temp .. pot. Temp. [K]
"""
kappa = Rl / cp_l
p_standard = 1000.0 # hPa
return T * (p_standard / p) ** kappa
[docs]def e_pot_temp(T, p, m):
"""
Equivalent potential temperature K
| T ... Temp in K
| p ... pressure in hPa
| m ... mixing ration in kg/kg (!)
| e_pot_temp ... equivalent potential temp. in K
"""
Th = pot_temp(T, p)
L = l_evw(T - Tice)
Th_e = Th * np.exp(L * m / cp_l / T)
return Th_e
[docs]def e_pot_temp_bolton(Tk, Tl, p, m):
"""
Equivalent potential temperature K after Bolton(1980):
The Computation of Eqivalent Potential Temperature (MWR Vol.108)
| Tk ... Temp in the starting level of ascend in K
| Tl ... Temp. in the lifting condensation level in K
| p ... pressure in hPa
| m ... mixing ration in kg/kg (!)
| e_pot_temp_bolton ... equivalent potential temp. in K"""
mg = m * 1000.0
A = Tk * (1000.0 / p) ** (0.2854 * (1 - 0.00028 * mg))
B = np.exp((3.376 / Tl - 0.00254) * (mg * (1 + 0.00081 * mg)))
return A * B
[docs]def t_lifting_condensation_level_bolton1(Tk, Td):
"""
Lifting condensation according to Bolton(1980):
"The Computation of Eqivalent Potential Temperature" (MWR Vol.108)
| Tk ... Temp in the starting level of ascend in K
| Td ... Dewpoint Temp. in the starting in the starting level of ascend
| in K
| p ... pressure in hPa
| m ... mixing ration in kg/kg (!)
| t_lifting_condensation_level_bolton1 ... lift. condens.level in K"""
return (1 / (1 / (Td - 56.0) + np.log(Tk / Td) / 800)) + 56.0
[docs]def showalter_temperature_bolton1(ThetaE, p, tswi):
"""
Showalter temperature [C] as a function of the eqivalent-potential
temperature according to Bolton(1980): "The Computation of Eqivalent
Potential Temperature" (MWR Vol.108) and pressure.
Terminates with an error message when there is no
convergence within 100 iterations.
| ThetaE ... equivalent (pseudo-) potential temp.
| p ... pressure in hPa
| tswi ... start value for the Showalter temperature as 0th-approximation
| [C]
| showalter_temperature_bolton1 ... Showalter temperature [C]"""
THETA_EPSILON = 1e-5
NMAX_ITER = 100
n = 0
while n < NMAX_ITER:
e = e_water(tswi)
m = m_mixingratio(e, p)
th_eb = e_pot_temp_bolton(tswi + Tice, tswi + Tice, p, m)
dtheta = ThetaE - th_eb
# print "%6.0f %6.1f %6.1f %7.3f" % (p,ThetaE,th_eb,tswi)
tswi = tswi + dtheta / 2.0
if abs(dtheta) < THETA_EPSILON:
return tswi
n = n + 1
print("\nERROR in function <showalter_temperature_Bolton1>:")
print("No convergence after %d iterations." % NMAX_ITER)
print("Try a better first guess for the SWI-temperature\n")
sys.exit(1)
[docs]def showalter_index_bolton1(
t_lower, t_upper, td_lower, p_lower=850.0, p_upper=500.0):
"""
Showalter index [C] as a function of pressure, temperature and
dewpoint at the lower level and pressure and temperature at
the upper level. Calculation is done according to Bolton(1980):
"The Computation of Eqivalent Potential Temperature" (MWR Vol.108).
Terminates with an error message when there is no
convergence within 100 iterations.
| t_lower ... temperature at p_lower in C
| t_upper ... temperature at p_upper in C
| td_lower ... dewpoint temperature at p_lower in C
| p_lower ... pressure at the lower level in hPa
| p_upper ... pressure at the upper level in hPa
| showalter_index_bolton1 ... Showalter Index (SWI) in C
SWI = Temp at the upper level - Showalter temperature for the upper level
"""
t_lift = t_lifting_condensation_level_bolton1(
t_lower +
Tice,
td_lower +
Tice)
ew_lower = e_water(td_lower)
mw_lower = m_mixingratio(ew_lower, p_lower)
th_e = e_pot_temp_bolton(t_lower + Tice, t_lift, p_lower, mw_lower)
t_swi_1guess = t_upper
t_swi = showalter_temperature_bolton1(th_e, p_upper, t_swi_1guess)
return t_upper - t_swi
[docs]def e_stefan_boltzmann(T, eps=1.0):
"""
Radiation flux density of a black or grey body
(Stefan-Boltzmann law)
| T .... temperature of the body [K]
| eps .. emissivity <= 1, defaults to 1
| e_stefan_boltzmann ... emitted radiation flux density W/m^2/K^4
"""
return eps * sigma * T ** 4
[docs]def planck_black_body_wl(T, lam, eps=1.):
"""Planck thermal radiance as a function of
| T .. temparature of the body [K]
| lam .. wavelenghth [m]
| eps .. spectral emissivity, 0 =< eps <=1"""
b = (2 * h_Planck * c_light ** 2 / lam ** 5 /
(np.exp(h_Planck * c_light / lam / k_Boltzmann / T) - 1.))
return b * eps
def __planck_black_body_fr__(T, nu, eps=1.):
"""Planck thermal radiance as a function of
| T ... temparature of the body [K]
| nu ... frequency [Hz]
| eps ... spectral emissivity, 0 =< eps <=1
"""
c1 = 2. * h_Planck * (nu ** 3) / (c_light ** 2)
c2 = h_Planck * nu / (k_Boltzmann * T)
b = c1 * 1. / (np.exp(c2) - 1.)
# x = _hc_k * nu / k_Boltzmann / T
# f = 15. / np.pi ** 4 * x ** 3 / (np.exp(x) - 1)
# b = sigma / np.pi * T ** 4 * f
return b * eps
[docs]def planck_black_body_fr(T, nu, eps=1.):
"""Planck thermal radiance as a function of
| T ... temparature of the body [K]
| nu ... frequency [Hz]
| eps ... spectral emissivity, 0 =< eps <=1
"""
lnB = (np.log(2 * h_Planck)
+3 * np.log(nu)
-2 * np.log(c_light)
-np.log(np.exp(h_Planck / k_Boltzmann * nu / T)))
return np.exp(lnB) * eps
[docs]def wind_dd2phi(dd):
"""Converts meteorological wind direction dd in degrees into the
mathematical direction phi in degrees."""
phi = np.mod(270. - dd, 360.)
return phi
[docs]def wind_phi2dd(phi):
"""Converts mathematical wind direction phi in degrees into the
meteorological direction ddd in degrees."""
dd = np.mod(270. - phi, 360.)
return dd
[docs]def wind_ddff2uv(dd, ff):
"""Convert wind given in (meteorological) direction 0..360 and
speed to its x- and y-components u and v"""
phi = wind_dd2phi(dd)
phir = np.radians(phi)
u = ff * np.cos(phir)
v = ff * np.sin(phir)
return u, v
[docs]def wind_uv2ddff(u, v):
"""Convert wind given in x- and y-components u and v
to the (meteorological) direction 0..360 and speed."""
phi = np.arctan2(u, v)
phid = np.degrees(phi)
dd = phid + 180.
ff = np.sqrt(u * u + v * v)
return dd, ff
[docs]def wind_kt2kmh(v):
"""Knots to km/h"""
return v * KT2KMH
[docs]def wind_kmh2kt(v):
"""km/h to knots"""
return v / KT2KMH
[docs]def wind_kt2ms(v):
"""knots 2 m/s"""
return v * KT2MS
[docs]def wind_ms2kt(v):
"""m/s to knots"""
return v / KT2MS
[docs]def wind_ms2kmh(v):
"""m/s to km/h"""
return v * 3.6
[docs]def wind_kmh2ms(v):
"""km/h to m/s"""
return v / 3.6
[docs]def wind_equal_hh_field(h0, h1, dh):
"""Returns an equidistant hh-field starting with
h0, ending with h1, stepping dh."""
return np.arange(h0, h1 + dh, dh)
[docs]def wind_mean_uv(hh, uh, vh):
"""Calculates from an height and cartesian wind component field
arithmetic mean component fields.
All field entries are treated equally"""
N = len(hh)
uz = 0.0
vz = 0.0
for k in range(N - 1):
dz = (hh[k + 1] - hh[k])
duz = 0.5 * (uh[k] + uh[k + 1]) * dz
dvz = 0.5 * (vh[k] + vh[k + 1]) * dz
uz = uz + duz
vz = vz + dvz
dh = hh[N - 1] - hh[0]
uq = uz / dh
vq = vz / dh
return uq, vq
[docs]def wind_mean_uv_equidistant(uh, vh):
"""Calculates from an height and cartesian wind component field
the arithmetic mean component fields for an equidistant windfield.
All field entries are treated equally"""
uq = np.mean(uh)
vq = np.mean(vh)
return uq, vq
[docs]def wind_single_w_inter(hx, hh, ww):
""""Linear interpolation between sampling points
For hx lower then the (geometrically)lowest entry is set to the
lowest entry, for hx higher then the (geometrically) highes entry
is set to the highes entry, meaning constant interpolation beyond
the boundaries."""
N = len(ww)
h0 = hh[0]
h1 = hh[N - 1]
wx = None
if hx <= h0: # lower equal than lowest
wx = ww[0]
return wx
if hx >= h1: # higher equal than highest
wx = ww[N - 1]
return wx
# Somewhere in between
for k, _h in enumerate(hh):
if (hx >= hh[k]) and (hx < hh[k + 1]):
wx = (ww[k] + (ww[k + 1] - ww[k]) / (hh[k + 1] -
hh[k]) * (hx - hh[k]))
return wx
[docs]def wind_interpolate_uv(h0, h1, dh, hh, uu, vv):
"""h0 .. lower interpolation boundary
h1 .. upper interpolation boundary,
dh .. interpolation intervall
hh .. height array for measured values
uu, vv .. cartesian wind component arrays at hh
Returns interpolated arrays hn .. heights,
un, vn .. cartesian wind components"""
unew = []
vnew = []
hnew = wind_equal_hh_field(h0, h1, dh)
# k = 0
for h in hnew:
ui = wind_single_w_inter(h, hh, uu)
vi = wind_single_w_inter(h, hh, vv)
unew.append(ui)
vnew.append(vi)
# k = k + 1
hn = np.array(hnew)
un = np.array(unew)
vn = np.array(vnew)
return hn, un, vn
[docs]def wind_print_huv(comment, h, u, v):
"""Text output for a cartesian wind field and conversion
to meteorological dd,ff components"""
print(comment)
print(' k: h u v | dd ff ')
print('-----------------------------------+-------------')
dd, ff = wind_uv2ddff(u, v)
N = len(h)
for k in range(N):
print('{:02d}: {:8.0f} {:8.2f} {:8.2f} | {:4.0f} {:7.2f}'.format(
k, h[k], u[k], v[k], dd[k], ff[k]))
return
[docs]def wind_print_raw(comment, hh, dd, ff):
"""Text output for measured raw wind data"""
uu, vv = wind_ddff2uv(dd, ff)
phi = wind_dd2phi(dd)
print(comment)
print(' k: h dd ff | phi u v')
print('---------------------------+-------------------------')
fmt = '{:02d}: {:6.0f} {:4.0f} {:7.1f} | {:4.0f} {:8.1f} {:8.1f}'
N = len(hh)
for k in range(N):
print(fmt.format(k, hh[k], dd[k], ff[k], phi[k], uu[k], vv[k]))
return
[docs]def wind_plot_uv(hh, uu, vv,
outpic='uv_h.png',
title='Windcomponents as function of height',
show=False):
"""Plots wind data
| hh .. array with height array
| uu .. array with cartesian-x componet of wind
| vv .. array with cartesian y-component of wind
|
| Returns None"""
plt.plot(hh, uu, label='u', c='r', lw=2)
plt.plot(hh, vv, label='v', c='g', lw=2)
plt.grid()
plt.xlabel('Height')
plt.ylabel('Velocity u,v')
plt.title(title)
plt.legend(loc='best', framealpha=0.5)
plt.savefig(outpic)
if show:
plt.show()
plt.clf()
return
def _testadmet():
"""A simple testing routine for module wmetpy"""
print('\n')
print(_testadmet.__doc__)
print("Version " + __version__)
phi = 47.0 + 32.0 / 60.0 # latitude of LOXA
z = 649.0 # Barometer height of LOXA
gg = g_welmec(phi, z)
gp = gpotheight(phi, z) # geopotential height (ICAO)
t = 18.9
T = t + Tice # Temperature C resp. K
td = 2.3
Td = td + Tice # Dew point temperature C resp. K
Tlift = t_lifting_condensation_level_bolton1(t + Tice, td + Tice)
p = 923.3 # Air pressure
e = e_water(td) # actual water vapor
q = q_spechum(e, p) # specific humdity
m = m_mixingratio(e, p) # mixing ratio
Tv = t_virt(T, q)
tv = Tv - Tice # virt. Temp. K
Tv2 = t_virt2(t, td, p)
tv2 = Tv2 - Tice # virt. temp. K
Tvquer = (Tv + (Tv + gammastand * gp)) / 2.0 # mean virt. temperature
tvquer = Tvquer - Tice
p0i = p_iso(p, Tvquer, -gp) # Reduced pressure isotherm
p0p = p_poly(p, T, -gp, gammastand) # -"- polytrop
Theta = pot_temp(T, p) # Potential Temperature
ThetaE = e_pot_temp(T, p, m)
ThetaE_Bolt = e_pot_temp_bolton(T, Tlift, p, m)
t_showalter = showalter_temperature_bolton1(ThetaE_Bolt, 500.0, t)
swi = showalter_index_bolton1(t, -18.0, td, p, 500.0)
epsil = 0.9 # emissivity
# Emitted therm. radiation flux dens.
E_stbo = e_stefan_boltzmann(T, epsil)
print(30 * "=")
print("Rl ... %7.2f" % Rl)
print("cp_l ... %7.2f" % cp_l)
print("cv_l ... %7.2f" % cv_l)
print(30 * "-")
print("phi ... %6.1f" % (phi))
print("z ... %6.2f" % (z))
print(30 * "-")
print("t ... %6.1f" % (t))
print("td(C) ... %6.1f" % (td))
print("Td(K) ... %6.2f" % Td)
print("p ... %5.1f" % (p))
print(30 * "-")
print("G_welmec . %9.7f" % (gg))
print("gp ... %6.2f" % (gp))
print("e ... %6.2f" % (e))
print("q ... %7.5f" % (q))
print("m ... %7.5f" % (m))
print("tv ... %4.1f" % (tv))
print("tv2 ... %4.1f" % (tv2))
print("tvquer ... %4.1f" % (tvquer))
print(30 * "-")
print("p0 iso ... %5.1f" % (p0i))
print("p0 poly .. %5.1f" % (p0p))
print(30 * "-")
print("T ... %5.1f" % (T))
print("p ... %5.1f" % (p))
print("PotTemp .. %5.1f" % (Theta))
print("EquvPoT .. %5.1f" % (ThetaE))
print("T_lift ... %5.1f" % (Tlift))
print("EquvPoT_Bolt %5.1f" % (ThetaE_Bolt))
print("500 hPa Showalter \ntemp. ... %5.1f" % (t_showalter))
print("SWI ... %5.1f" % (swi))
print(30 * "-")
print("T ... %5.1f" % (T))
print("eps ... %5.1f" % (epsil))
print("E_StefBoltz %5.1f" % (E_stbo))
print(75 * "=")
ff = 10.
fmt_ddff = ("(dd, ff)=({0:3d}/{1:4.1f}) --> (u, v)=" +
"({2:+5.1f}/{3:+5.1f}) --> (dd, ff)=({0:3d}/{1:4.1f})")
for dd in range(0, 360 + 30, 30):
u, v = wind_ddff2uv(dd, ff)
d0, f0 = wind_uv2ddff(u, v)
print(fmt_ddff.format(dd, ff, u, v, d0, f0))
print(75 * "=")
if __name__ == "__main__":
_testadmet()