# -*- coding: utf-8 -*-
"""
@author: jeremy leconte
"""
import numpy as np
import numba
from .convection import dry_convective_adjustment_numba, turbulent_diffusion_numba
from exo_k.util.molar_mass import Molar_mass
from exo_k.util.cst import PI, N_A
from exo_k.aerosol.util_aerosol import mmr_to_number_density_ratio
[docs]
class Tracers(object):
def __init__(self, settings, tracers=None, tracer_values=None, Dmol=0.,
bg_vmr=None, M_bg=None, Nlay=None, **kwargs):
"""Class that deals with tracers.
Fills out the tracers.qarray and creates tracer_names, a table of correspondence
between tracer names and indices in tracers.qarray.
"""
self.settings=settings
if tracers is None:
self.Ntrac=1
tracers = {'inactive_tracer':{}}
else:
self.Ntrac = len(tracers)
if tracer_values is None:
tracer_values = {}
if Nlay is not None:
self.Nlay = Nlay
else:
raise RuntimeError("We need to know Nlay to initialize Tracers")
self.bg_vmr = bg_vmr.copy()
self.gas_vmr = bg_vmr.copy()
self.var_gas_idx = list()
self.aerosols_idx = list()
self.var_gas_names = list()
self.aerosols_names = list()
self.gas_molar_masses = list()
self.aerosols_bulk_density = list()
self.aerosols_r_eff = list()
self.some_var_gases = False
self.some_active_aerosol = False
self.dico=tracers
self.namelist = list(tracers.keys())
self.idx = dict()
self.qarray = np.empty((self.Ntrac, self.Nlay))
self.qsurf = np.zeros(self.Ntrac)
self.qdeep = - np.ones(self.Ntrac)
for ii, name in enumerate(self.namelist):
self.idx[name]=ii
if name in tracer_values.keys():
self.qarray[ii]=np.copy(tracer_values[name])
elif 'init_value' in self.dico[name].keys():
self.qarray[ii]=np.ones(self.Nlay)*self.dico[name]['init_value']
else:
self.qarray[ii]=np.zeros(self.Nlay)
if 'type' in self.dico[name]:
if self.dico[name]['type'] in ('gas', 'vapor'):
self.some_var_gases = True
self.var_gas_idx.append(ii)
self.var_gas_names.append(name)
self.gas_molar_masses.append(Molar_mass().fetch(name))
if self.dico[name]['type'] in ('aerosol'):
self.some_active_aerosol = True
self.aerosols_idx.append(ii)
self.aerosols_names.append(name)
self.aerosols_bulk_density.append(self.dico[name]['bulk_density'])
self.aerosols_r_eff.append(self.dico[name]['r_eff'])
if 'q_deep' in self.dico[name]:
self.qdeep[ii] = np.copy(self.dico[name]['q_deep'])
if 'surface_reservoir' in self.dico[name]: # should be reserved to condensates
self.qsurf[ii] = np.copy(self.dico[name]['surface_reservoir'])
self.settings.set_parameters(surface_reservoir = True)
self.var_gas_idx = np.array(self.var_gas_idx)
self.var_gas_names = np.array(self.var_gas_names)
self.gas_molar_masses = np.array(self.gas_molar_masses)
self.aerosols_idx = np.array(self.aerosols_idx)
self.aerosols_names = np.array(self.aerosols_names)
self.aerosols_bulk_density = np.array(self.aerosols_bulk_density)
self.aerosols_r_eff = np.array(self.aerosols_r_eff)
self.M_bg = M_bg
self.update_gas_composition()
self.Dmol = Dmol
[docs]
def update_gas_composition(self, update_vmr=True):
"""Performs mass to volume mixing ratio conversion,
computes the new molar mass of the total gas
and transmits new composition to the radiative model.
!!!Small discrepancies with case without tracers!!!
Maybe a problem with the background gas. JL22 => is this still true ?
Parameters
----------
atm: exo_k.Atm object
If atm is provided, its composition will be
updated with the new composition.
"""
qvar = np.zeros(self.Nlay)
ovMvar = np.zeros(self.Nlay)
for ii, idx in enumerate(self.var_gas_idx):
qvar += self.qarray[idx]
ovMvar += self.qarray[idx]/self.gas_molar_masses[ii]
ovMbg=1./self.M_bg
self.Mgas = 1./((1.-qvar)*ovMbg + ovMvar)
if update_vmr:
var_vmr_tot = 0.
for ii, idx in enumerate(self.var_gas_idx):
vmr_tmp = np.core.umath.maximum(
self.Mgas * self.qarray[idx] / self.gas_molar_masses[ii], 0.)
#conversion to vmr
self.gas_vmr[self.var_gas_names[ii]] = vmr_tmp
var_vmr_tot += vmr_tmp
#print(vmr_tmp, var_vmr_tot, isinstance(vmr_tmp, (np.ndarray)))
var_vmr_tot = 1. - var_vmr_tot
#print(var_vmr_tot)
for mol, vmr in self.bg_vmr.items():
#print(mol, vmr, var_vmr_tot, isinstance(var_vmr_tot, (np.ndarray)))
self.gas_vmr[mol] = vmr * var_vmr_tot
#JL23 careful, if we have a variable gas that is also in the background mix, they are not added.
# vmr for bg gas replaces the other.
[docs]
def update_aerosol_properties(self, atm):
if self.some_active_aerosol:
aer_reffs_densities = {}
#atm.compute_mass_density()
for ii, idx in enumerate(self.aerosols_idx):
mmr_rad = np.sqrt(self.qarray[idx,1:] * self.qarray[idx,:-1])
aer_reffs_densities[self.aerosols_names[ii]] = [self.aerosols_r_eff[ii],
mmr_to_number_density_ratio(mmr_rad, atm.Mgas_rad, self.aerosols_r_eff[ii], self.aerosols_bulk_density[ii])]
else:
aer_reffs_densities = None
return aer_reffs_densities
[docs]
def dry_convective_adjustment(self, timestep, Htot, atm, verbose=False):
"""Computes convective adjustement.
Parameters
----------
timestep: float
physical timestep of the current step (in s/cp).
(needs to be converted before it is sent to `turbulent diffusion)
Htot: array, np.ndarray
Total heating rate (in W/kg) of all physical processes
already computed
atm: :class:`Atm` object
The Atm object used in the radiative transfer which
contains many state variables.
"""
new_t = atm.tlay + timestep * Htot
if self.settings['moist_inhibition']:
Mgas_tmp = self.Mgas
else:
Mgas_tmp = np.mean(self.Mgas)
Mgas_tmp = np.full_like(self.Mgas, Mgas_tmp)
H_conv, q_array, self.index_dry_convective_top = dry_convective_adjustment_numba(timestep, self.Nlay, new_t,
atm.exner, atm.dmass, self.qarray, Mgas_tmp, verbose=verbose)
if self.settings['convective_transport']:
self.qarray=q_array
return H_conv
[docs]
def update_surface_reservoir(self, condensing_pairs_idx = False,
surf_layer_mass = 0.):
"""Update surface reservoirs of tracers.
E.g. deals with removing excess condensates from the first layer or
putting more when everything as been evaporated.
"""
if condensing_pairs_idx: # an empty list is False
qcond_surf_layer = self.settings['qcond_surf_layer']
for idx_vap, idx_cond in condensing_pairs_idx:
dm_surf_layer = (qcond_surf_layer - self.qarray[idx_cond, -1]) * surf_layer_mass
if dm_surf_layer > self.qsurf[idx_cond]: #empties surface reservoir
self.qarray[idx_cond, -1] += self.qsurf[idx_cond] / surf_layer_mass
self.qsurf[idx_cond] = 0.
else: # general case
self.qarray[idx_cond, -1] = qcond_surf_layer
self.qsurf[idx_cond] -= dm_surf_layer
[docs]
def interpolate_tracers_profile(self, new_logplay, old_logplay):
"""Re interpolates tracers on a new pressure grid
Parameers
---------
new_logplay: array
new pressure grid
old_logplay: array
Old pressure grid
"""
self.Nlay = new_logplay.shape[0]
new_qarray = np.empty((self.Ntrac, self.Nlay))
for ii in range(self.Ntrac):
new_qarray[ii] = np.interp(new_logplay, old_logplay, self.qarray[ii])
self.qarray = new_qarray
self.update_gas_composition(update_vmr=True)
def __getitem__(self, key):
return self.qarray[self.idx[key]]