Source code for exo_k.atm_evolution.condensation_gcm

# -*- coding: utf-8 -*-
"""
@author: jeremy leconte
"""
import numpy as np
import numba
import exo_k.util.cst as cst


[docs] class Condensing_species(object): def __init__(self, Latent_heat_vaporization = 0., cp_vap = 0., Mvap = 0., T_ref = 0., Psat_ref = 0., cp_cond = None): """ Parameters ---------- Latent_heat_vaporization: float specific Latent heat vaporization (J/kg) cp_vap: float Specific heat capacity of the vapor (J/kg/K) Mvap: float Molar mass of vapor (kg/mol) T_ref: float Reference temperature Psat_ref: float Saturation vapor pressure at the reference temperature (Pa) cp_cond: float (optional) Specific heat capacity of the condensate (J/kg/K). Assumed equal to cp_vap if not provided. """ self.Latent_heat_vaporization = Latent_heat_vaporization self.cp_vap = cp_vap self.T_ref = T_ref self.Psat_ref = Psat_ref self.Mvap = Mvap self.Rvap = cst.RGP / self.Mvap if cp_cond is None: self.cp_cond = cp_vap else: self.cp_cond = cp_cond self.delta_cp = self.cp_vap - self.cp_cond self.delta_cp_R = self.delta_cp / self.Rvap self.LovR = self.Latent_heat_vaporization / self.Rvap self.c1 = self.LovR / self.T_ref - self.delta_cp_R self.c2 = self.delta_cp_R * self.T_ref - self.LovR
[docs] def Lvap(self, T): """Latent heat at temperature T Parameters ---------- T: array Temperature in layers (K) """ return Lvap_T(T, self.Latent_heat_vaporization, self.T_ref, self.delta_cp)
[docs] def Psat(self, T): """Saturation vapor pressure for the condensing species Parameters ---------- T: array Temperature in layers (K) """ return Psat_T(T, self.T_ref, self.Psat_ref, self.c1, self.c2, self.delta_cp_R)
[docs] def Tsat(self, P): """Boiling temperature for the condensing species NOT EXACT IF DELTA CP != 0. Parameters ---------- P: array Pressure in layers (Pa) """ return Tsat_P(P, self.Psat_ref, self.c1, self.c2)
[docs] def qsat(self, psat, p, epsilon): """Saturation vapor mass mixing ratio for the condensing species Parameters ---------- psat: array saturation vapor pressure p: array pressure at the layer center epsilon: float or array Ratio of the molar mass of the vapor over the background molar mass """ return Qsat(psat, p, epsilon)
[docs] def dPsat_dT(self, T): """Saturation vapor pressure derivative for the condensing species Parameters ---------- T: array Temperature in layers (K) """ return dPsat_dT(T, self.Latent_heat_vaporization, self.T_ref, self.Psat_ref, self.Rvap, self.delta_cp, self.delta_cp_R, self.c1, self.c2)
[docs] def dlnPsat_dlnT(self, T): """Saturation vapor pressure for the condensing species and its derivative Parameters ---------- T: array Temperature in layers (K) """ return dlnPsat_dlnT(T, self.Latent_heat_vaporization, self.T_ref, self.delta_cp, self.Rvap)
[docs] def moist_adiabat(self, T, P, cp, Mgas): """Computes the threshold thermal gradient (d ln T / d ln P) for a moist atmosphere. Parameters ---------- T: array Temperature in layers (K) P: array pressure at the layer center cp: float specific heat capacity at constant pressure of the background gas Mgas: float or array Molar mass of the background atmosphere Returns ------- array Moist adiabat lapse rate Lvap: array Latent heat at temperature T psat: array Saturation vapor pressure for the condensing species (Pa) qsat: array Saturation vapor mass mixing ratio for the condensing species dqsat_dt: array Derivative of qsat with respect to temperature at fixed pressure q_crit: array Critical mass mixing ratio for the inhibition of moist convection (Eq. 17 of Leconte et al. 2017) """ return moist_adiabat(T, P, cp, Mgas, self.cp_vap, self.Mvap, self.Rvap, self.Latent_heat_vaporization, self.T_ref, self.Psat_ref, self.delta_cp, self.delta_cp_R, self.c1, self.c2)
[docs] def compute_condensation_parameters(self, T, P, Mgas): """Computes necessary quantities to compute large scale condensation. Parameters ---------- T: array Temperature in layers (K) P: array pressure at the layer center Mgas: float or array Molar mass of the background atmosphere """ return compute_condensation_parameters(T, P, Mgas, self.Mvap, self.Rvap, self.Latent_heat_vaporization, self.T_ref, self.Psat_ref, self.delta_cp, self.delta_cp_R, self.c1, self.c2)
def __repr__(self): """Method to output parameters """ output=""" Lvap : {lvap} cp (vap) : {cp} T0 (K): {t0} Psat(T0) (Pa): {p0}""".format(lvap=self.Latent_heat_vaporization, cp=self.cp_vap, t0=self.T_ref, p0=self.Psat_ref) return output
[docs] class Condensation_Thermodynamical_Parameters(object): def __init__(self, Latent_heat_vaporization = 0., cp_vap = 0., Mvap = 0., T_ref = 0., Psat_ref = 0., cp_cond = None): """ Parameters ---------- Latent_heat_vaporization: float specific Latent heat vaporization (J/kg) cp_vap: float Specific heat capacity of the vapor (J/kg/K) Mvap: float Molar mass of vapor (kg/mol) T_ref: float Reference temperature Psat_ref: float Saturation vapor pressure at the reference temperature (Pa) cp_cond: float (optional) Specific heat capacity of the condensate (J/kg/K). Assumed equal to cp_vap if not provided. """ Rvap = cst.RGP / Mvap if cp_cond is None: cp_cond = cp_vap else: cp_cond = cp_cond delta_cp = cp_vap - cp_cond delta_cp_R = delta_cp / Rvap LovR = Latent_heat_vaporization / Rvap c1 = LovR / T_ref - delta_cp_R c2 = delta_cp_R * T_ref - LovR self.th_params = np.array([cp_vap, Mvap, Rvap, Latent_heat_vaporization, T_ref, Psat_ref, delta_cp, delta_cp_R, c1, c2])
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def Lvap_T(T, Latent_heat_vaporization, T_ref, delta_cp): """Latent heat at temperature T Parameters ---------- T: array Temperature in layers (K) """ return Latent_heat_vaporization * np.ones_like(T)
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def Psat_T(T, T_ref, Psat_ref, c1, c2, delta_cp_R): """GCM version for test """ Pref_solid_liquid=611.14 Trefvaporization=35.86 Trefsublimation=7.66 Tmin=8. r3vaporization=17.269 r3sublimation=21.875 T_h2O_ice_liq=273.16 psat = np.where(T>T_h2O_ice_liq, Pref_solid_liquid*np.exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization)),# liquid / vapour Pref_solid_liquid*np.exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation)) # solid / vapour ) return psat
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def Tsat_P(P, Psat_ref, c1, c2): """GCM version for test """ Pref_solid_liquid=611.14 Trefvaporization=35.86 Trefsublimation=7.66 Tmin=8. r3vaporization=17.269 r3sublimation=21.875 T_h2O_ice_liq=273.16 Tsat = np.where(P<Pref_solid_liquid, (T_h2O_ice_liq*r3sublimation- Trefsublimation*np.log(P/Pref_solid_liquid))/(r3sublimation-np.log(P/Pref_solid_liquid)), (T_h2O_ice_liq*r3vaporization- Trefvaporization*np.log(P/Pref_solid_liquid))/(r3vaporization-np.log(P/Pref_solid_liquid)) ) return Tsat
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def Qsat(psat, p, epsilon): """Saturation vapor mass mixing ratio for the condensing species Parameters ---------- psat: array saturation vapor pressure p: array pressure at the layer center epsilon: float or array Ratio of the molar mass of the vapor over the background molar mass """ psat_tmp = np.core.umath.minimum(psat,p) fac = p + (epsilon -1.) * psat_tmp qsat = epsilon * psat_tmp / fac return qsat
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def dPsat_dT(T, Latent_heat_vaporization, T_ref, Psat_ref, Rvap, delta_cp, delta_cp_R, c1, c2): """Saturation vapor pressure derivative for the condensing species Parameters ---------- T: array Temperature in layers (K) """ Lvap = Lvap_T(T, Latent_heat_vaporization, T_ref, delta_cp) psat = Psat_T(T, T_ref, Psat_ref, c1, c2, delta_cp_R) Trefvaporization=35.86 Trefsublimation=7.66 r3vaporization=17.269 r3sublimation=21.875 T_h2O_ice_liq=273.16 dpsat = np.where(T>T_h2O_ice_liq, psat*r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2,# liquid / vapour psat*r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2 # solid / vapour ) return dpsat, psat, Lvap
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def dlnPsat_dlnT(T, Latent_heat_vaporization, T_ref, delta_cp, Rvap): """Saturation vapor pressure for the condensing species and its derivative Parameters ---------- T: array Temperature in layers (K) """ Lvap = Lvap_T(T, Latent_heat_vaporization, T_ref, delta_cp) Trefvaporization=35.86 Trefsublimation=7.66 r3vaporization=17.269 r3sublimation=21.875 T_h2O_ice_liq=273.16 dlnpsat_dt = np.where(T>T_h2O_ice_liq, r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2,# liquid / vapour r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2 # solid / vapour ) return dlnpsat_dt*T, Lvap
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def moist_adiabat(T, P, cp, Mgas, cp_vap, Mvap, Rvap, Latent_heat_vaporization, T_ref, Psat_ref, delta_cp, delta_cp_R, c1, c2): """Computes the threshold thermal gradient (d ln T / d ln P) for a moist atmosphere. Parameters ---------- T: array Temperature in layers (K) P: array pressure at the layer center cp: float specific heat capacity at constant pressure of the background gas Mgas: float or array Molar mass of the background atmosphere Returns ------- array Moist adiabat lapse rate Lvap: array Latent heat at temperature T psat: array Saturation vapor pressure for the condensing species (Pa) qsat: array Saturation vapor mass mixing ratio for the condensing species dqsat_dt: array Derivative of qsat with respect to temperature at fixed pressure q_crit: array Critical mass mixing ratio for the inhibition of moist convection (Eq. 17 of Leconte et al. 2017) """ epsilon = Mvap / Mgas psat = Psat_T(T, T_ref, Psat_ref, c1, c2, delta_cp_R) qsat = Qsat(psat, P, epsilon) dlpsat_dlT, Lvap = dlnPsat_dlnT(T, Latent_heat_vaporization, T_ref, delta_cp, Rvap) dqsat_dt = qsat**2 * P * dlpsat_dlT / (epsilon*psat*T) qa=1.-qsat qLvt = qsat * Lvap / T # fac = qsat * cp_vap + qa * cp + qLvt * dlpsat_dlT fac = cp + qLvt * dlpsat_dlT q_crit = epsilon / ((epsilon - 1.) * dlpsat_dlT) return ( cst.RGP / Mgas + qLvt ) / fac, Lvap, psat, qsat, dqsat_dt, q_crit # is missing p/p_a terms
# return ( cst.RGP / Mgas + qLvt ) / fac*0.8, Lvap, psat, qsat, dqsat_dt, q_crit # is missing p/p_a terms # factor 0.8 to mimic 3D
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def compute_condensation_parameters(T, P, Mgas, Mvap, Rvap, Latent_heat_vaporization, T_ref, Psat_ref, delta_cp, delta_cp_R, c1, c2): """Computes necessary quantities to compute large scale condensation. Parameters ---------- T: array Temperature in layers (K) P: array pressure at the layer center Mgas: float or array Molar mass of the background atmosphere """ epsilon = Mvap / Mgas psat = Psat_T(T, T_ref, Psat_ref, c1, c2, delta_cp_R) qsat = Qsat(psat, P, epsilon) dlpsat_dlT, Lvap = dlnPsat_dlnT(T, Latent_heat_vaporization, T_ref, delta_cp, Rvap) dqsat_dt = qsat**2 * P * dlpsat_dlT / (epsilon*psat*T) return Lvap, qsat, dqsat_dt