Skip to content

Commit

Permalink
added a function to calculate climatological tropopause
Browse files Browse the repository at this point in the history
  • Loading branch information
@FranziskaWinterstein
FranziskaWinterstein committed Dec 20, 2023
1 parent d7bbb23 commit 5da7f8a
Show file tree
Hide file tree
Showing 2 changed files with 122 additions and 51 deletions.
81 changes: 54 additions & 27 deletions esmvaltool/diag_scripts/lifetime/lifetime.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,6 +273,7 @@
# from iris.analysis.cartography import area_weights
# from iris.analysis.maths import log, exp
from iris.coord_categorisation import add_year
from iris.util import broadcast_to_shape
# from iris.coords import AuxCoord
from matplotlib.gridspec import GridSpec
from matplotlib.ticker import FormatStrFormatter, LogLocator, NullFormatter
Expand Down Expand Up @@ -549,21 +550,6 @@ def _load_data_and_calculate_coefficients(self):
if cube.coords('time', dim_coords=True):
ih.unify_time_coord(cube)

# Set Z-coordinate
if cube.coords('air_pressure', dim_coords=True):
self.z_coord = 'air_pressure'
z_coord = cube.coord('air_pressure', dim_coords=True)
z_coord.attributes['positive'] = 'down'
z_coord.convert_units('hPa')
elif cube.coords('atmosphere_hybrid_sigma_pressure_coordinate', dim_coords=True):
self.z_coord = 'atmosphere_hybrid_sigma_pressure_coordinate'
z_coord = cube.coord('atmosphere_hybrid_sigma_pressure_coordinate', dim_coords=True)
z_coord.attributes['positive'] = 'down'
add_model_level(cube)
elif cube.coords('altitude', dim_coords=True):
self.z_coord = 'alititude'
z_coord = cube.coord('altitude', dim_coords=True)
z_coord.attributes['positive'] = 'up'

# cube for each variable
variables[variable['short_name']] = cube
Expand All @@ -572,13 +558,37 @@ def _load_data_and_calculate_coefficients(self):

oxidant = {ox: variables[ox] for ox in name_oxidant}

input_data_dataset[dataset]['variables'] = variables
input_data_dataset[dataset]['reaction'] = (
reaction = (
self._calculate_reaction(oxidant,
rho,
variables['ta'],
name_reactant))
input_data_dataset[dataset]['weight'] = self._define_weight(variables)

weight = self._define_weight(variables)

# Set Z-coordinate
if reaction.coords('air_pressure', dim_coords=True):
self.z_coord = 'air_pressure'
z_coord = reaction.coord('air_pressure', dim_coords=True)
z_coord.attributes['positive'] = 'down'
z_coord.convert_units('hPa')
elif reaction.coords('atmosphere_hybrid_sigma_pressure_coordinate',
dim_coords=True):
self.z_coord = 'atmosphere_hybrid_sigma_pressure_coordinate'
z_coord = reaction.coord('atmosphere_hybrid_sigma_pressure_coordinate',
dim_coords=True)
z_coord.attributes['positive'] = 'down'
add_model_level(reaction)
add_model_level(weight)
else:
raise NotImplementedError(f"Lifetime calculation is not implemented"
" for the present type of vertical coordinate.")

input_data_dataset[dataset]['z_coord'] = z_coord

input_data_dataset[dataset]['variables'] = variables
input_data_dataset[dataset]['reaction'] = reaction
input_data_dataset[dataset]['weight'] = weight

return input_data_dataset

Expand Down Expand Up @@ -622,16 +632,17 @@ def _calculate_rho(self, variables):
Calculate number density with respect to the given input
variables.
"""
# model levels
if 'grmassdry' in variables and 'grvol' in variables:
#if self.z_coord == 'lev' and 'grmassdry' in variables and 'grvol' in variables:
#if self.z_coord == 'atmosphere_hybrid_sigma_pressure_coordinate'
# and 'grmassdry' in variables and 'grvol' in variables:
rho = self._number_density_dryair_by_grid(
variables['grmassdry'],
variables['grvol'])
elif ('press' in variables and
'ta' in variables and
# pressure levels
elif ('ta' in variables and
'hus' in variables):
rho = self._number_density_dryair_by_press(
variables['press'],
variables['ta'],
variables['hus'])
else:
Expand All @@ -640,7 +651,7 @@ def _calculate_rho(self, variables):
"Provide either:\n"
" - grmassdry and grvol\n"
" or\n"
" - press, ta and hus")
" - ta and hus")
# dummy
rho = variables['press']

Expand All @@ -664,20 +675,31 @@ def _calculate_reaction(self, oxidant, rho, ta, name_reactant):

return reaction

def _number_density_dryair_by_press(self, press, ta, hus):
def _number_density_dryair_by_press(self, ta, hus, press=None):
"""
Calculate number density of dry air.
Used to convert from mol / mol_dry into molec / cm3
by using present pressure, temperature and
humidity.
by using present temperature and humidity.
##qqq
Should there be an option to provide the simulated pressure field
rather than the derived pressure from interpolation?
###
Used constants:
- N_A Avogrado constant from scipy
- R gas constant from scipy
- m_air Molarmass of Air
- m_h2o Molarmass of watervapor
"""

if not press:
press = broadcast_to_shape(
ta.coord('air_pressure').points,
ta.shape,
ta.coord_dims('air_pressure')
)
rho = ((press * N_A * 10**(-6)) /
(R * ta * (1. + (self.m_air / self.m_h2o - 1.) * hus)))
# correct metadata
Expand Down Expand Up @@ -729,7 +751,7 @@ def _calculate_reaction_rate(self, ta, ox, re):
reaction_rate.data = coeff_a * np.exp(coeff_b * np.log(ta.data)
- (coeff_er / ta.data))

# standard reaction rate
# standard reaction rate (arrhenius)
else:
reaction_rate.data = coeff_a * np.exp(-(coeff_er / ta.data))

Expand Down Expand Up @@ -1065,6 +1087,11 @@ def create_timeseries_plot(self, region, input_data, base_datasets):
# Plot original time series
plot_kwargs = self._get_plot_kwargs(plot_type, base_datasets[label])
plot_kwargs['axes'] = axes

if self.plots[plot_type]['display_mean'] is not False:
mean = cube.collapsed('time', iris.analysis.MEAN)
plot_kwargs['label'] = f"{plot_kwargs['label']} ({mean.data:.2f})"

iris.plot.plot(cube, **plot_kwargs)

# Plot annual means if desired
Expand Down
92 changes: 68 additions & 24 deletions esmvaltool/diag_scripts/lifetime/lifetime_base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
import os
import re

from pprint import pprint

import sys
import cartopy
import iris
Expand Down Expand Up @@ -39,35 +41,64 @@ def calculate_lifetime(dataset, plot_type, region):
def extract_region(dataset, region, case='reaction'):
"""Return cube with everything outside region set to zero.
If z_coord is defined as:
- air_pressure, use:
- ptp and air_pressure
- tp_clim and air_pressure
- atmosphere_hybrid_sigma_pressure_coordinate, use:
- tp_i and model_level_number
- ptp and (derived) air_pressure
- tp_clim and (derived) air_pressure
Current aware regions:
- TROP: troposphere (excl. tropopause), requires tropopause pressure
- STRA: stratosphere (incl. tropopause), requires tropopause pressure
"""
var = dataset[case]
z_coord = dataset['z_coord']

pprint(dataset)

if not 'ptp' in dataset['variables'] and not 'tp_i' in dataset['variables']:
tp_clim = climatological_tropopause(broadcast_to_shape(
var.coord('latitude').points,
var.shape,
var.coord_dims(z_coord)
))

if region in ['TROP', 'STRA']:
## - can I choose the troposphere by the preprocessor?
if 'tp_i' in dataset['variables']:
tp = dataset['variables']['tp_i']
z_coord = 'model_level_number'
elif 'tp' in dataset['variables']:
tp = dataset['variables']['tp']
z_coord = 'air_pressure'
else:
raise NotImplementedError(
"Guessing the tropopause pressure is currently not implemented."
" You need to provide a variabel containing the tropopause pressure.")

use_z_coord = 'air_pressure'
if z_coord.name() == 'air_pressure':
if 'ptp' in dataset['variables']:
tp = dataset['variables']['ptp']
else:
tp = tp_clim
elif z_coord.name() == 'atmosphere_hybrid_sigma_pressure_coordinate':
if 'tp_i' in dataset['variables']:
tp = dataset['variables']['tp_i']
use_z_coord = 'model_level_number'
elif 'ptp' in dataset['variables']:
tp = dataset['variables']['ptp']
else:
tp = tp_clim

z_4d = broadcast_to_shape(
var.coord(z_coord).points,
var.shape,
var.coord_dims(z_coord)
)
tp_4d = broadcast_to_shape(
tp.data,
var.coord(use_z_coord).points,
var.shape,
var.coord_dims('time') + var.coord_dims('latitude') + var.coord_dims('longitude'),
var.coord_dims(use_z_coord)
)
# Problem: ValueError: dim_map must have an entry for every dimension of the input array
# Maybe there is a more elegant way to do this
if not var.shape == tp.shape:
tp_4d = broadcast_to_shape(
tp.data,
var.shape,
var.coord_dims('time') + var.coord_dims('latitude') + var.coord_dims('longitude'),
)
else:
tp_4d = tp

if region == 'TROP':
var.data = np.ma.array(
Expand All @@ -84,16 +115,25 @@ def extract_region(dataset, region, case='reaction'):

return var

def climatological_tropopause(latitude):
"""Return cube with climatological tropopause pressure.
"""
tp_clim = 300. - 215. * (np.cos(latitude) ** 2)

return tp_clim

def sum_up_to_plot_dimensions(var, plot_type):
"""
Returns the cube summed over the appropriate dimensions
"""
if var.coords('air_pressure', dim_coords=True):
z_coord = var.coords('air_pressure', dim_coords=True)
elif var.coords('lev', dim_coords=True):
z_coord = var.coords('lev', dim_coords=True)
elif var.coords('atmosphere_hybrid_sigma_pressure_coordinate', dim_coords=True):
z_coord = var.coords('atmosphere_hybrid_sigma_pressure_coordinate', dim_coords=True)[0]
if plot_type in ['timeseries', 'annual_cycle']:
if var.coords('air_pressure', dim_coords=True):
z_coord = var.coords('air_pressure', dim_coords=True)
elif var.coords('lev', dim_coords=True):
z_coord = var.coords('lev', dim_coords=True)
elif var.coords('atmosphere_hybrid_sigma_pressure_coordinate', dim_coords=True):
z_coord = var.coords('atmosphere_hybrid_sigma_pressure_coordinate', dim_coords=True)[0]

if plot_type == 'timeseries':
cube = var.collapsed(['longitude', 'latitude', z_coord], iris.analysis.SUM)
Expand All @@ -103,7 +143,11 @@ def sum_up_to_plot_dimensions(var, plot_type):
cube = var.collapsed(['longitude', 'latitude'], iris.analysis.SUM)
elif plot_type == 'annual_cycle':
# TODO!
cube = var.collapsed(['longitude', 'latitude', z_coord], iris.analysis.SUM)
# not iris.analysis.SUM but some kind of mean
#cube = var.collapsed(['longitude', 'latitude', z_coord], iris.analysis.SUM)
raise NotImplementedError(f"The sum to plot dimensions for plot_type {plot_type}"
" is currently not implemented")


return cube

Expand Down

0 comments on commit 5da7f8a

Please sign in to comment.