Make derivation of total column ozone (toz) more flexible and add derivation of stratospheric and tropospheric column ozone

esmvalcore/cmor/tables/custom/CMOR_soz.dat

@@ -0,0 +1,22 @@
+!============
+variable_entry:    soz
+!============
+modeling_realm:    atmos
+!----------------------------------
+! Variable attributes:
+!----------------------------------
+standard_name:     equivalent_thickness_at_stp_of_atmosphere_ozone_content
+units:             m
+cell_methods:      time: mean
+cell_measures:     area: areacella
+long_name:         Stratospheric Ozone Column (O3 mole fraction >= 125 ppb)
+comment:           stratospheric ozone column calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU. Here, the stratosphere is defined as the region where O3 mole fraction >= 125 ppb.
+!----------------------------------
+! Additional variable information:
+!----------------------------------
+dimensions:        longitude latitude time
+type:              real
+valid_min:         0.0
+valid_max:         5000.0
+!----------------------------------
+!

esmvalcore/cmor/tables/custom/CMOR_toz.dat

@@ -1,17 +1,17 @@
-SOURCE: CCMI1
+SOURCE: CMIP6
 !============
 variable_entry:    toz
 !============
 modeling_realm:    atmos
 !----------------------------------
 ! Variable attributes:
 !----------------------------------
-standard_name:
-units:             DU
+standard_name:     equivalent_thickness_at_stp_of_atmosphere_ozone_content
+units:             m
 cell_methods:      time: mean
 cell_measures:     area: areacella
 long_name:         Total Ozone Column
-comment:           total ozone column in DU
+comment:           Total ozone column calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU.
 !----------------------------------
 ! Additional variable information:
 !----------------------------------

esmvalcore/cmor/tables/custom/CMOR_tozStderr.dat

@@ -1,17 +1,17 @@
-SOURCE: CCMI1
+SOURCE: CMIP6
 !============
 variable_entry:    tozStderr
 !============
 modeling_realm:    atmos
 !----------------------------------
 ! Variable attributes:
 !----------------------------------
-standard_name:
-units:             DU
+standard_name:     equivalent_thickness_at_stp_of_atmosphere_ozone_content
+units:             m
 cell_methods:      time: mean
 cell_measures:     area: areacella
 long_name:         Total Ozone Column Error
-comment:           total ozone column in DU
+comment:           Total ozone column error calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU.
 !----------------------------------
 ! Additional variable information:
 !----------------------------------

esmvalcore/cmor/tables/custom/CMOR_troz.dat

@@ -0,0 +1,22 @@
+!============
+variable_entry:    troz
+!============
+modeling_realm:    atmos
+!----------------------------------
+! Variable attributes:
+!----------------------------------
+standard_name:     equivalent_thickness_at_stp_of_atmosphere_ozone_content
+units:             m
+cell_methods:      time: mean
+cell_measures:     area: areacella
+long_name:         Tropospheric Ozone Column (O3 mole fraction < 125 ppb)
+comment:           tropospheric ozone column calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU. Here, the troposphere is defined as the region where O3 mole fraction < 125 ppb.
+!----------------------------------
+! Additional variable information:
+!----------------------------------
+dimensions:        longitude latitude time
+type:              real
+valid_min:         0.0
+valid_max:         5000.0
+!----------------------------------
+!

esmvalcore/preprocessor/_derive/_shared.py

@@ -170,15 +170,15 @@ def _create_pressure_array(cube, ps_cube, top_limit):
 def _get_pressure_level_widths(array, air_pressure_axis=1):
     """Compute pressure level widths.
 
-    For a 1D array with pressure level columns, return a 1D array with
-    pressure level widths.
+    For array with pressure level columns, return array with pressure
+    level widths.
 
     """
     array = np.copy(array)
     if np.any(np.diff(array, axis=air_pressure_axis) > 0.0):
         raise ValueError("Pressure level value increased with height")
 
-    # Calculate centers
+    # Calculate array of centers between two neighboring pressure levels
     indices = [slice(None)] * array.ndim
     array_shifted = np.roll(array, -1, axis=air_pressure_axis)
     index_0 = deepcopy(indices)

esmvalcore/preprocessor/_derive/soz.py

@@ -0,0 +1,107 @@
+"""Derivation of variable ``soz``."""
+
+import dask.array as da
+import iris
+
+from ._baseclass import DerivedVariableBase
+from .toz import DerivedVariable as Toz
+from .toz import add_longitude_coord, interpolate_hybrid_plevs
+
+# O3 mole fraction threshold (in ppb) that is used for the definition of the
+# stratosphere (stratosphere = region where O3 mole fraction is at least as
+# high as the threshold value)
+STRATOSPHERIC_O3_THRESHOLD = 125.0
+
+
+class DerivedVariable(DerivedVariableBase):
+    """Derivation of variable ``soz``."""
+
+    @staticmethod
+    def required(project):
+        """Declare the variables needed for derivation."""
+        if project == 'CMIP6':
+            required = [{'short_name': 'o3'}]
+        else:
+            required = [{'short_name': 'tro3'}]
+        return required
+
+    @staticmethod
+    def calculate(cubes):
+        """Compute stratospheric column ozone.
+
+        Note
+        ----
+        Here, the stratosphere is defined as the region in which the O3 mole
+        fraction is at least as high as the given threshold
+        (``STRATOSPHERIC_O3_THRESHOLD``).
+
+        In the calculation of ``toz``, the surface air pressure (``ps``) is
+        used to determine the pressure level width of the lowest layer. For
+        ``soz``, this lowest layer can be ignored since it is not located in
+        the stratosphere (it will be masked out due to the O3 mole fraction
+        threshold). Thus, the surface air pressure (``ps``) is not necessary
+        for the derivation of ``soz`` and is simply replaced with the lowest
+        pressure level in the data to be able to use the ``toz`` derivation
+        function.
+
+        The calculation of ``soz`` consists of three steps:
+        (1) Mask out O3 mole fractions smaller than given threshold.
+        (2) Cut out the lowest pressure level from the data and use it as
+            surface air pressure (``toz``).
+        (3) Use derivation function of ``toz`` to calculate ``soz`` (using the
+            masked data).
+
+        """
+        o3_cube = cubes.extract_cube(
+            iris.Constraint(name='mole_fraction_of_ozone_in_air')
+        )
+
+        # If o3 is given on hybrid pressure levels (e.g., from Table AERmon),
+        # interpolate it to regular pressure levels
+        if len(o3_cube.coord_dims('air_pressure')) > 1:
+            o3_cube = interpolate_hybrid_plevs(o3_cube)
+
+        # To support zonal mean o3 (e.g., from Table AERmonZ), add longitude
+        # coordinate if necessary
+        if not o3_cube.coords('longitude'):
+            o3_cube = add_longitude_coord(o3_cube)
+
+        # (1) Mask O3 mole fraction using the given threshold
+        o3_cube.convert_units('1e-9')
+        mask = o3_cube.lazy_data() < STRATOSPHERIC_O3_THRESHOLD
+        mask |= da.ma.getmaskarray(o3_cube.lazy_data())
+        o3_cube.data = da.ma.masked_array(o3_cube.lazy_data(), mask=mask)
+
+        # (2) Add surrogate for the surface air pressure (ps) cube using the
+        # lowest pressure level available in the data (this is fine since the
+        # the lowest pressure level is far away from the stratosphere).
+
+        # Get dummy ps cube with correct dimensions
+        ps_dims = (o3_cube.coord_dims('time') +
+                   o3_cube.coord_dims('latitude') +
+                   o3_cube.coord_dims('longitude'))
+        idx_to_extract_ps = [0] * o3_cube.ndim
+        for ps_dim in ps_dims:
+            idx_to_extract_ps[ps_dim] = slice(None)
+        ps_cube = o3_cube[tuple(idx_to_extract_ps)].copy()
+
+        # Set ps data using lowest pressure level available and add correct
+        # metadata
+        lowest_plev = o3_cube.coord('air_pressure').points.max()
+        ps_data = da.broadcast_to(lowest_plev, ps_cube.shape)
+        ps_cube.data = ps_data
+        ps_cube.var_name = 'ps'
+        ps_cube.standard_name = 'surface_air_pressure'
+        ps_cube.long_name = 'Surface Air Pressure'
+        ps_cube.units = o3_cube.coord('air_pressure').units
+
+        # Cut lowest pressure level from o3_cube
+        z_dim = o3_cube.coord_dims('air_pressure')[0]
+        idx_to_cut_lowest_plev = [slice(None)] * o3_cube.ndim
+        idx_to_cut_lowest_plev[z_dim] = slice(1, None)
+        o3_cube = o3_cube[tuple(idx_to_cut_lowest_plev)]
+
+        # (3) Use derivation function of toz to calculate soz using the masked
+        # o3 cube and the surrogate ps cube
+        cubes = iris.cube.CubeList([o3_cube, ps_cube])
+        return Toz.calculate(cubes)

esmvalcore/preprocessor/_derive/toz.py

@@ -1,9 +1,14 @@
-"""Derivation of variable `toz`."""
+"""Derivation of variable ``toz``."""
+
+import warnings
 
 import cf_units
 import iris
 from scipy import constants
 
+from esmvalcore.cmor.table import CMOR_TABLES
+
+from .._regrid import extract_levels, regrid
 from ._baseclass import DerivedVariableBase
 from ._shared import pressure_level_widths
 
@@ -19,16 +24,53 @@
 DOBSON_UNIT = cf_units.Unit('2.69e20 m^-2')
 
 
+def add_longitude_coord(cube, ps_cube=None):
+    """Add dimensional ``longitude`` coordinate of length 1 to cube."""
+    lon_coord = iris.coords.DimCoord([180.0], bounds=[[0.0, 360.0]],
+                                     var_name='lon',
+                                     standard_name='longitude',
+                                     long_name='longitude',
+                                     units='degrees_east')
+    new_dim_coords = [(c, cube.coord_dims(c)) for c in cube.dim_coords]
+    new_dim_coords.append((lon_coord, cube.ndim))
+    new_aux_coords = [(c, cube.coord_dims(c)) for c in cube.aux_coords]
+    new_cube = iris.cube.Cube(
+        cube.core_data()[..., None],
+        dim_coords_and_dims=new_dim_coords,
+        aux_coords_and_dims=new_aux_coords,
+    )
+    new_cube.metadata = cube.metadata
+    return new_cube
+
+
+def interpolate_hybrid_plevs(cube):
+    """Interpolate hybrid pressure levels."""
+    # Use CMIP6's plev19 target levels (in Pa)
+    target_levels = CMOR_TABLES['CMIP6'].coords['plev19'].requested
+    cube.coord('air_pressure').convert_units('Pa')
+    cube = extract_levels(cube, target_levels, 'linear',
+                          coordinate='air_pressure')
+    return cube
+
+
 class DerivedVariable(DerivedVariableBase):
-    """Derivation of variable `toz`."""
+    """Derivation of variable ``toz``."""
 
     @staticmethod
     def required(project):
         """Declare the variables needed for derivation."""
+        # TODO: make get_required _derive/__init__.py use variables as argument
+        # and make this dependent on mip
         if project == 'CMIP6':
-            required = [{'short_name': 'o3'}, {'short_name': 'ps'}]
+            required = [
+                {'short_name': 'o3'},
+                {'short_name': 'ps', 'mip': 'Amon'},
+            ]
         else:
-            required = [{'short_name': 'tro3'}, {'short_name': 'ps'}]
+            required = [
+                {'short_name': 'tro3'},
+                {'short_name': 'ps'},
+            ]
         return required
 
     @staticmethod
@@ -41,24 +83,53 @@ def calculate(cubes):
         upper integration bound of 0 Pa is used.
 
         """
-        tro3_cube = cubes.extract_cube(
-            iris.Constraint(name='mole_fraction_of_ozone_in_air'))
+        o3_cube = cubes.extract_cube(
+            iris.Constraint(name='mole_fraction_of_ozone_in_air')
+        )
         ps_cube = cubes.extract_cube(
-            iris.Constraint(name='surface_air_pressure'))
+            iris.Constraint(name='surface_air_pressure')
+        )
+
+        # If o3 is given on hybrid pressure levels (e.g., from Table AERmon),
+        # interpolate it to regular pressure levels
+        if len(o3_cube.coord_dims('air_pressure')) > 1:
+            o3_cube = interpolate_hybrid_plevs(o3_cube)
+
+        # To support zonal mean o3 (e.g., from Table AERmonZ), add longitude
+        # coordinate and collapsed ps cube if necessary to ensure that they
+        # have correct shapes
+        if not o3_cube.coords('longitude'):
+            o3_cube = add_longitude_coord(o3_cube)
+            ps_cube = ps_cube.collapsed('longitude', iris.analysis.MEAN)
+            ps_cube.remove_coord('longitude')
+            ps_cube = add_longitude_coord(ps_cube)
+
+        # If the horizontal dimensions of ps and o3 differ, regrid ps
+        # Note: regrid() checks if the regridding is really necessary before
+        # running the actual interpolation
+        ps_cube = regrid(ps_cube, o3_cube, 'linear')
 
-        p_layer_widths = pressure_level_widths(tro3_cube,
+        # Actual derivation of toz using o3 mole fraction and pressure level
+        # widths
+        p_layer_widths = pressure_level_widths(o3_cube,
                                                ps_cube,
                                                top_limit=0.0)
-        toz_cube = (tro3_cube * p_layer_widths / STANDARD_GRAVITY * MW_O3 /
+        toz_cube = (o3_cube * p_layer_widths / STANDARD_GRAVITY * MW_O3 /
                     MW_AIR)
-        toz_cube = toz_cube.collapsed('air_pressure', iris.analysis.SUM)
-        toz_cube.units = (tro3_cube.units * p_layer_widths.units /
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', category=UserWarning,
+                message='Collapsing a non-contiguous coordinate')
+            toz_cube = toz_cube.collapsed('air_pressure', iris.analysis.SUM)
+        toz_cube.units = (o3_cube.units * p_layer_widths.units /
                           STANDARD_GRAVITY_UNIT * MW_O3_UNIT / MW_AIR_UNIT)
 
-        # Convert from kg m^-2 to Dobson unit (2.69e20 m^-2 )
+        # Convert from kg m^-2 to Dobson units DU (2.69e20 m^-2 ) and from
+        # DU to m (1 mm = 100 DU)
         toz_cube = toz_cube / MW_O3 * AVOGADRO_CONST
         toz_cube.units = toz_cube.units / MW_O3_UNIT * AVOGADRO_CONST_UNIT
         toz_cube.convert_units(DOBSON_UNIT)
-        toz_cube.units = 'DU'
+        toz_cube.data = toz_cube.core_data() * 1e-5
+        toz_cube.units = 'm'
 
         return toz_cube