Add PRIMAVERA Eady Growth Rate diagnostic

esmvaltool/config-references.yml

@@ -245,6 +245,11 @@ authors:
   lledo_llorenc:
     name: Lledo, Llorenc
     institute: BSC, Spain
+    orcid: 
+    email: 
+  loosveldt-tomas_saskia:
+    name: Loosveldt-Tomas, Saskia
+    institute: BSC, Spain
     orcid:
   lorenz_ruth:
     name: Lorenz, Ruth
@@ -338,6 +343,10 @@ authors:
     name: Roehrig, Romain
     institute: MeteoFr, France
     orcid: https://orcid.org/0000-0002-3903-3841
+  sanchez-gomez_emilia:
+    name: Sanchez, Emilia
+    institute: CERFACS, France
+    orcid:
   serva_federico:
     name: Serva, Federico
     institute: CNR, Italy
@@ -690,6 +699,7 @@ references:
   uwisc: "O'Dell et al., J. Clim., doi:10.1175/2007JCLI1958.1, 2008."
   woa: "Locarnini et al., World Ocean Atlas 2013, Volume 1: Temperature, NOAA Atlas NESDIS 73, 2013."
 
+
 projects:
   c3s-magic: Copernicus Climate Change Service 34a Lot 2 (MAGIC) project
   climval: BMBF MiKlip Project ClimVal
@@ -703,6 +713,7 @@ projects:
   eval4cmip: DLR and University of Bremen project funded by the Initiative and Networking Fund of the Helmholtz Society
   ewatercycle: eWaterCycle project
   qa4ecv: QA4ECV
+  primavera: EU H2020 project PRIMAVERA
   trr181: DFG Project TRR-181
   ukesm: UKESM, UK Earth System Model project (NERC)
 

esmvaltool/diag_scripts/primavera/eady_growth_rate.py

@@ -0,0 +1,181 @@
+import os
+import logging
+
+import iris
+import iris.cube
+import iris.analysis
+import iris.util
+
+import numpy as np
+
+from dask import array as da
+
+import esmvaltool.diag_scripts.shared
+import esmvaltool.diag_scripts.shared.names as n
+from esmvaltool.diag_scripts.shared import group_metadata
+from esmvaltool.diag_scripts.shared._base import ProvenanceLogger
+
+
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+class EadyGrowthRate(object):
+    def __init__(self, config):
+        self.cfg = config
+        self.filenames = esmvaltool.diag_scripts.shared.Datasets(self.cfg)
+        self.fill_value = 1e20
+        self.ref_p = 1000.0
+        self.g = 9.80665
+        self.con = 0.3098
+        self.omega = 7.292e-5
+
+    def compute(self):
+        data = group_metadata(self.cfg['input_data'].values(), 'alias')
+        for alias in data:
+            var = group_metadata(data[alias], 'short_name')
+            ta = iris.load_cube(var['ta'][0]['filename'])
+            plev = ta.coord('air_pressure')
+
+            theta = self.potential_temperature(ta, plev)
+
+            del ta
+
+            zg = iris.load_cube(var['zg'][0]['filename'])
+
+            brunt = self.brunt_vaisala_frq(theta, zg)
+
+            lats = zg.coord('latitude')
+            fcor = self.coriolis(lats, zg.shape)
+
+            ua = iris.load_cube(var['ua'][0]['filename'])
+
+            egr = self.eady_growth_rate(fcor, ua, zg, brunt)
+
+            cube_egr = ua.copy(egr * 86400)
+
+            cube_egr = cube_egr.collapsed('time', iris.analysis.MEAN)
+
+            self.save(cube_egr, alias, data)
+
+    def potential_temperature(self, ta, plev):
+        p0 = iris.coords.AuxCoord(self.ref_p,
+                                  long_name='reference_pressure',
+                                  units='hPa')
+        p0.convert_units(plev.units)
+        p = (p0.points/plev.points)**(2/7)
+        theta = ta * iris.util.broadcast_to_shape(
+            p,
+            ta.shape,
+            ta.coord_dims('air_pressure')
+            )
+        theta.long_name = 'potential_air_temperature'
+
+        return theta
+
+    def vertical_integration(self, x, y):
+        # Perform a non-cyclic centered finite-difference to integrate
+        # along pressure levels
+
+        plevs = x.shape[1]
+
+        dxdy_0 = ((x[:, 1, :, :].lazy_data() - x[:, 0, :, :].lazy_data()) /
+                  (y[:, 1, :, :].lazy_data() - y[:, 0, :, :].lazy_data()))
+
+        dxdy_centre = ((x[:, 2:plevs, :, :].lazy_data() -
+                        x[:, 0:plevs-2, :, :].lazy_data()) /
+                       (y[:, 2:plevs, :, :].lazy_data() -
+                        y[:, 0:plevs-2, :, :].lazy_data()))
+
+        dxdy_end = ((x[:, plevs-1, :, :].lazy_data() -
+                     x[:, plevs-2, :, :].lazy_data()) /
+                    (y[:, plevs-1, :, :].lazy_data() -
+                     y[:, plevs-2, :, :].lazy_data()))
+
+        bounds = [dxdy_end, dxdy_0]
+        stacked_bounds = da.stack(bounds, axis=1)
+        total = [dxdy_centre, stacked_bounds]
+
+        # Concatenate arrays where the last slice is dxdy_0
+        dxdy = da.concatenate(total, axis=1)
+
+        # Move dxdy_0 to the beggining of the array
+        dxdy = da.roll(dxdy, 1, axis=1)
+
+        return dxdy
+
+    def brunt_vaisala_frq(self, theta, zg):
+        dthdz = self.vertical_integration(theta, zg)
+        dthdz = da.where(dthdz > 0, dthdz, 0)
+        buoy = (self.g / theta.lazy_data()) * dthdz
+        brunt = da.sqrt(buoy)
+        brunt = da.where(brunt != 0, brunt, self.fill_value)
+
+        return brunt
+
+    def coriolis(self, lats, ndim):
+        fcor = 2.0 * self.omega * np.sin(np.radians(lats.points))
+        fcor = fcor[np.newaxis, np.newaxis, :, np.newaxis]
+        fcor = da.broadcast_to(fcor, ndim)
+
+        return fcor
+
+    def eady_growth_rate(self, fcor, ua, zg, brunt):
+        dudz = self.vertical_integration(ua, zg)
+        egr = self.con * abs(fcor) * abs(dudz) / brunt
+
+        return egr
+
+    def save(self, egr, alias, data):
+        egr.standard_name = None
+        egr.long_name = 'eady_growth_rate'
+        egr.var_name = 'egr'
+        egr.units = 'day-1'
+
+        script = self.cfg[n.SCRIPT]
+        project = data[alias][0]['project']
+        dataset = data[alias][0]['dataset']
+        exp = data[alias][0]['exp']
+        start = data[alias][0]['start_year']
+        end = data[alias][0]['end_year']
+        output_name = '{project}_' \
+                      '{dataset}_' \
+                      '{exp}_' \
+                      '{script}_' \
+                      '{start}_' \
+                      '{end}.nc'.format(project=project,
+                                        dataset=dataset,
+                                        exp=exp,
+                                        script=script,
+                                        start=start,
+                                        end=end)
+        output_file = os.path.join(self.cfg[n.WORK_DIR], output_name)
+        iris.save(egr, output_file)
+
+        caption = ("{script} between {start} and {end}"
+                   "according to {dataset}").format(
+                       script=script.split('_'),
+                       start=start,
+                       end=end,
+                       dataset=dataset
+                   )
+        for i in range(len(data[alias])):
+            ancestors = [data[alias][i]['filename']]
+
+        record = {
+            'caption': caption,
+            'domains': ['global'],
+            'autors': ['sanchez-gomez_emilia'],
+            'references': ['acknow_project'],
+            'ancestors': ancestors
+            }
+        with ProvenanceLogger(self.cfg) as provenance_logger:
+            provenance_logger.log(output_file, record)
+
+
+def main():
+    with esmvaltool.diag_scripts.shared.run_diagnostic() as config:
+        EadyGrowthRate(config).compute()
+
+
+if __name__ == "__main__":
+    main()

esmvaltool/recipes/recipe_eady_growth_rate.yml

@@ -0,0 +1,69 @@
+# ESMValTool
+---
+documentation:
+  description: |
+    Recipe to compute the annual, summer and winter Eady Growth Rates. The
+    output produces netcdf files for each model and season.
+
+
+  authors:
+    - sanchez-gomez_emilia
+
+  maintainer:
+    - loosveldt-tomas_saskia
+
+  projects:
+    - primavera
+
+datasets:
+  - {dataset: CNRM-CM5,  project: CMIP5, exp: historical,  ensemble: r1i1p1,
+     mip: Amon, start_year: 1979,  end_year: 2005}
+
+preprocessors:
+  annual_mean:
+    annual_statistics:
+  summer_mean:
+    extract_season:
+      season: 'JJA'
+    seasonal_statistics:
+  winter_mean:
+    extract_season:
+      season: 'DJF'
+    seasonal_statistics:
+
+diagnostics:
+  annual_egr:
+    variables:
+      ta:
+        preprocessor: annual_mean
+      zg:
+        preprocessor: annual_mean
+      ua:
+        preprocessor: annual_mean
+    scripts:
+      annual_eady_growth_rate:
+        script: primavera/eady_growth_rate.py
+
+  summer_egr:
+    variables:
+      ta:
+        preprocessor: summer_mean
+      zg:
+        preprocessor: summer_mean
+      ua:
+        preprocessor: summer_mean
+    scripts:
+      summer_eady_growth_rate:
+        script: primavera/eady_growth_rate.py
+
+  winter_egr:
+    variables:
+      ta:
+        preprocessor: winter_mean
+      zg:
+        preprocessor: winter_mean
+      ua:
+        preprocessor: winter_mean
+    scripts:
+      winter_eady_growth_rate:
+        script: primavera/eady_growth_rate.py