#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) Météo France (2014-)
# This software is governed by the CeCILL-C license under French law.
# http://www.cecill.info
"""
This module contains functions to:
- to convert vertical coordinates systems from the one to another;
- compute air specific gas constant R according to specific humidity
and hydrometeors.
- to plot vertical profiles.
"""
from __future__ import print_function, absolute_import, unicode_literals, division
import numpy
from bronx.meteo import constants
from .config import default_Ptop
# Constants
Cpd = constants.Cpd
Rd = constants.Rd
Rv = constants.Rv
g0 = constants.g0
# FORMULAS #
############
[docs]def hybridP2fluxpressure(A, B, Psurf):
"""
Computes and returns the pressure at flux levels defined by hybrid-pressure
coefficients. Coefficients are assumed to define flux levels.
:param A: table of A coefficients
:param B: table of B coefficients
:param Psurf: surface pressure in Pa.
A and B must not contain the first coefficient (A=B=0.)
"""
if not len(A) == len(B):
raise ValueError("A, B must have the same size.")
if not isinstance(Psurf, numpy.ndarray):
if isinstance(Psurf, int) or isinstance(Psurf, float):
Psurf = numpy.array([Psurf])
else:
Psurf = numpy.array(Psurf)
pi_tilde = numpy.zeros([len(A)] + list(Psurf.shape))
# computation
for k in range(len(A)):
pi_tilde[k] = A[k] + B[k] * Psurf
if not numpy.all(pi_tilde[1:] >= pi_tilde[:-1]):
raise ValueError('pi_tilde array must be in ascending order.')
return pi_tilde
[docs]def hybridP2masspressure(A, B, Psurf, vertical_mean):
"""
Computes and returns the pressure at mass levels defined by hybrid-pressure
coefficients. Coefficients are assumed to define flux levels.
:param A: table of A coefficients
:param B: table of B coefficients
:param Psurf: surface pressure in Pa.
:param vertical_mean: kind of vertical averaging
A and B must not contain the first coefficient (A=B=0.)
"""
pi_tilde = hybridP2fluxpressure(A, B, Psurf)
pi = flux2masspressures(pi_tilde, vertical_mean)
if not numpy.all(pi[1:] >= pi[:-1]):
raise ValueError('pi array must be in ascending order.')
return pi
[docs]def hybridH2fluxheight(A, B, Zsurf, conv2height=False):
"""
Computes and returns the altitude at flux levels defined by hybrid-height
coefficients. Coefficients are assumed to define flux levels.
:param A: table of A coefficients
:param B: table of B coefficients
:param Zsurf: surface altitude in m.
:param conv2height: True to compute height instead of altitude
A and B must include the underground level
"""
if not len(A) == len(B):
raise ValueError("A, B must have the same size.")
if not isinstance(Zsurf, numpy.ndarray):
if isinstance(Zsurf, int) or isinstance(Zsurf, float):
Zsurf = numpy.array([Zsurf])
else:
Zsurf = numpy.array(Zsurf)
H_tilde = numpy.zeros([len(A)] + list(Zsurf.shape))
# computation
for k in range(len(A)):
H_tilde[k] = A[k] + B[k] * Zsurf
if conv2height:
H_tilde -= Zsurf
if not numpy.all(H_tilde[1:] >= H_tilde[:-1]):
raise ValueError('H_tilde array must be in ascending order.')
return H_tilde
[docs]def hybridH2massheight(A, B, Zsurf, conv2height=False):
"""
Computes and returns the altitude at mass levels defined by hybrid-height
coefficients. Coefficients are assumed to define flux levels.
:param A: table of A coefficients
:param B: table of B coefficients
:param Zsurf: surface altitude in m.
:param conv2height: True to compute height instead of altitude
A and B must include the underground level
"""
if not len(A) == len(B):
raise ValueError("A, B must have the same size.")
H_tilde = hybridH2fluxheight(A, B, Zsurf, conv2height=conv2height)
H = flux2massheights(H_tilde)
if not numpy.all(H[1:] >= H[:-1]):
raise ValueError('H array must be in ascending order.')
return H
[docs]def flux2masspressures(pi_tilde, vertical_mean, Ptop=default_Ptop):
"""Converts pressures at flux levels to mass levels."""
if not numpy.all(pi_tilde[1:] >= pi_tilde[:-1]):
raise ValueError('pi_tilde array must be in ascending order.')
L = len(pi_tilde)
if not isinstance(pi_tilde, numpy.ndarray):
pi_tilde = numpy.array(pi_tilde)
pi = numpy.zeros(pi_tilde.shape)
for k in range(1, L + 1):
ik = k - 1 # python arranging
if vertical_mean == 'geometric':
if k == 1:
pi[ik] = (pi_tilde[ik] - Ptop) / (1. + Cpd / Rd)
else:
pi[ik] = numpy.sqrt(pi_tilde[ik] * pi_tilde[ik - 1])
elif vertical_mean == 'arithmetic':
if k == 1:
pi[ik] = (pi_tilde[ik] + Ptop) / 2.
else:
pi[ik] = (pi_tilde[ik] + pi_tilde[ik - 1]) / 2.
else:
raise NotImplementedError("vertical_mean not among" +
" ('geometric', 'arithmetic').")
return pi
[docs]def mass2fluxpressures(pi, vertical_mean, Ptop=default_Ptop):
"""Converts pressures at mass levels to flux levels."""
if not numpy.all(pi[1:] >= pi[:-1]):
raise ValueError('pi array must be in ascending order.')
L = len(pi)
if not isinstance(pi, numpy.ndarray):
pi = numpy.array(pi)
pi_tilde = numpy.zeros(pi.shape)
for k in range(1, L + 1):
ik = k - 1 # python arranging
if vertical_mean == 'geometric':
if k == 1:
pi_tilde[ik] = (pi[ik] + Ptop) * (1. + Cpd / Rd)
else:
pi_tilde[ik] = pi[ik] ** 2 / pi_tilde[ik - 1]
elif vertical_mean == 'arithmetic':
if k == 1:
pi_tilde[ik] = (pi[ik] + Ptop) / 2
else:
pi_tilde[ik] = (pi[ik] + pi[ik - 1]) / 2.
else:
raise NotImplementedError("vertical_mean not among" +
" ('geometric', 'arithmetic').")
return pi_tilde
[docs]def flux2massheights(H_tilde):
"""Converts altitudes at flux levels to mass levels."""
if not numpy.all(H_tilde[1:] >= H_tilde[:-1]):
raise ValueError('H_tilde array must be in ascending order.')
if not isinstance(H_tilde, numpy.ndarray):
H_tilde = numpy.array(H_tilde)
H = numpy.zeros(H_tilde.shape)
H[:-1] = (H_tilde[:-1] + H_tilde[1:]) / 2.
H[-1] = 2 * H_tilde[-1] - H_tilde[-2]
return H[1:-1]
[docs]def mass2fluxheights(H):
"""Converts altitudes at mass levels to flux levels."""
if not numpy.all(H[1:] >= H[:-1]):
raise ValueError('H array must be in ascending order.')
# L = len(H)
if not isinstance(H, numpy.ndarray):
H = numpy.array(H)
H_tilde = numpy.zeros(H.shape)
raise NotImplementedError('mass2fluxheights is not yet implemented.')
# for k in range(1, L+1):
# ik = k-1 # python arranging
# if k == 1:
# pi_tilde[ik] = (pi[ik] + Ptop) * (1. + Cpd/Rd)
# else:
# pi_tilde[ik] = pi[ik]**2 / pi_tilde[ik-1]
return H_tilde
[docs]def pressure2altitude(R, T, vertical_mean,
pi=None, pi_tilde=None, Pdep=None, Phi_surf=None):
"""
Converts a pressure profile to mass-levels altitude (= height above ground
if *Phi_surf == 0*).
The pressure can be given at mass levels (*pi*) or flux levels (*pi_tilde*),
or both (assuming they are consistent), with unit: Pa.
:param pi: mass levels pressure
:param pi_tilde: flux levels pressure
:param R: table of R = specific gas constant (J/kg/K)
:param T: table of Temperature (K)
:param vertical_mean: kind of vertical averaging
:param Pdep: table of Pressure departure: P-Pi, where P is the hydrostatic
pressure and Pi the Non-Hydrostatic pressure (Pa).
0. if Hydrostatic (default)
:param Phi_surf: surface geopotential (J/kg)
"""
# get full pressures
if pi is None:
pi = flux2masspressures(pi_tilde, vertical_mean)
elif pi_tilde is None:
pi_tilde = mass2fluxpressures(pi, vertical_mean)
else:
if not isinstance(pi, numpy.ndarray):
pi = numpy.array(pi)
if not isinstance(pi_tilde, numpy.ndarray):
pi_tilde = numpy.array(pi_tilde)
L = len(pi)
if not len(pi) == len(pi_tilde) == len(R) == len(T):
raise ValueError("R, T, pi, pi_tilde must have the same size.")
if not isinstance(R, numpy.ndarray):
R = numpy.array(R)
if not isinstance(T, numpy.ndarray):
T = numpy.array(T)
if Phi_surf is None:
myPhi_surf = numpy.zeros(T[0].shape)
else:
if not isinstance(Phi_surf, numpy.ndarray):
myPhi_surf = numpy.array(Phi_surf)
else:
myPhi_surf = Phi_surf
if Pdep is None:
myPdep = numpy.zeros(R.shape)
else:
if not isinstance(Pdep, numpy.ndarray):
myPdep = numpy.array(Pdep)
else:
myPdep = Pdep
if not pi.shape == pi_tilde.shape:
raise ValueError("pi and pi_tilde must have the same shape.")
if not R.shape == T.shape == myPdep.shape:
raise ValueError("R, T, and Pdep must have the same shape.")
if (not T[0].shape == myPhi_surf.shape) and (not T[0].size == myPhi_surf.size == 1):
raise ValueError("Phi_surf shape must be compatible with other shapes.")
# Pre-computation
delta = numpy.zeros(pi.shape)
alpha = numpy.zeros(pi.shape)
for k in range(1, L + 1):
ik = k - 1 # python arranging
if k == 1:
delta[ik] = 1. + Cpd / Rd
alpha[ik] = 1.
else:
delta[ik] = (pi_tilde[ik] - pi_tilde[ik - 1]) / pi[ik]
alpha[ik] = 1 - pi[ik] / pi_tilde[ik]
# Geopotential
Phi = numpy.zeros(T.shape)
partialsum = numpy.zeros(T.shape)
for k in reversed(range(1, L + 1)):
ik = k - 1 # python arranging
if k == L:
partialsum[ik] = 0.
else:
partialsum[ik] = (partialsum[ik + 1] +
R[ik + 1] * T[ik + 1] * delta[ik + 1] /
(1. + myPdep[ik + 1] / pi[ik + 1]))
Phi[ik] = myPhi_surf + partialsum[ik] + R[ik] * T[ik] * alpha[ik] / \
(1. + myPdep[ik] / pi[ik])
# Altitude
z = Phi / g0
return z
[docs]def hybridP2altitude(A, B, R, T, Psurf, vertical_mean,
Pdep=None, Phi_surf=None, Ptop=default_Ptop):
"""
Computes the altitude of mass levels defined by hybrid-pressure
coefficients. (= height above ground if *Phi_surf == 0*).
:param A: table of A coefficients
:param B: table of B coefficients
:param Psurf: surface pressure in Pa.
:param R: table of R = specific gas constant (J/kg/K)
:param T: table of Temperature (K)
:param vertical_mean: kind of vertical averaging
:param Pdep: table of Pressure departure: P-Pi, where P is the hydrostatic
pressure and Pi the Non-Hydrostatic pressure (Pa).
0. if Hydrostatic (default)
:param Phi_surf: surface geopotential (J/kg)
A and B must not contain the first coefficient (A=B=0.)
"""
if not len(A) == len(B) == len(R) == len(T):
raise ValueError("A, B, R, T must have the same size.")
pi_tilde = hybridP2fluxpressure(A, B, Psurf)
pi = flux2masspressures(pi_tilde, vertical_mean, Ptop=Ptop)
z = pressure2altitude(R, T, vertical_mean,
pi=pi, pi_tilde=pi_tilde, Pdep=Pdep,
Phi_surf=Phi_surf)
return z