Benchmarking recipes (Lauer et al.)

doc/sphinx/source/recipes/index.rst

@@ -21,6 +21,7 @@ large variety of input data.
 .. toctree::
    :maxdepth: 1
 
+   recipe_benchmarking
    recipe_model_evaluation
    recipe_monitor
    recipe_psyplot

doc/sphinx/source/recipes/recipe_benchmarking.rst

@@ -0,0 +1,140 @@
+.. _recipe_benchmarking:
+
+Model Benchmarking
+==================
+
+Overview
+--------
+
+These recipes and diagnostics are based on :ref:`recipe_monitor <recipe_monitor>`: that allow plotting arbitrary preprocessor output, i.e., arbitrary variables from arbitrary datasets. An extension of these diagnostics is used to benchmark a model simulation with other datasets (e.g. CMIP6). The benchmarking features are described in `Lauer et al.`_:.
+
+.. _`Lauer et al.`: A. Lauer, Bock, L., Hassler, B., Jöckel, P., Ruhe, L., and Schlund, M.: Monitoring and benchmarking Earth System Model simulations with ESMValTool v2.12.0, Geosci. Model Dev. (submitted).
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in `recipes/model_evaluation`
+
+* recipe_model_benchmarking_annual_cycle.yml
+* recipe_model_benchmarking_boxplots.yml
+* recipe_model_benchmarking_diurnal_cycle.yml
+* recipe_model_benchmarking_maps.yml
+* recipe_model_benchmarking_timeseries.yml
+* recipe_model_benchmarking_zonal.yml
+
+Diagnostics are stored in `diag_scripts/monitor/`
+
+* :ref:`multi_datasets.py
+  <api.esmvaltool.diag_scripts.monitor.multi_datasets>`:
+  Monitoring diagnostic to show multiple datasets in one plot (incl. biases).
+
+
+Recipe settings
+~~~~~~~~~~~~~~~
+
+See :ref:`multi_datasets.py<api.esmvaltool.diag_scripts.monitor.multi_datasets>`: for a list of all possible configuration options that can be specified in the recipe.
+
+.. note::
+   Please note that exactly one dataset (the dataset to be benchmarked) needs to specify the facet ``benchmark_dataset: True`` in the dataset entry of the recipe. For line plots (i.e. annual cycle, seasonal cycle, diurnal cycle, time series), it is recommended, to specify a particular line color and line style in the ``scripts`` section of the recipe for the dataset to be benchmarked (``benchmark_dataset: True``) so that this dataset is easy to identify in the plot. In the example below, MIROC6 is the dataset to be benchmarked and ERA5 is used as a reference dataset.
+
+.. code-block:: yaml
+
+   scripts:
+         allplots:
+           script: monitor/multi_datasets.py
+           plot_folder: '{plot_dir}'
+           plot_filename: '{plot_type}_{real_name}_{mip}'
+           group_variables_by: variable_group
+           facet_used_for_labels: alias
+           plots:
+             diurnal_cycle:
+               annual_mean_kwargs: False
+               legend_kwargs:
+                 loc: upper right
+               plot_kwargs:
+                 'MIROC6':
+                   color: red
+                   label: '{alias}'
+                   linestyle: '-'
+                   linewidth: 2
+                   zorder: 4
+                 ERA5:
+                   color: black
+                   label: '{dataset}'
+                   linestyle: '-'
+                   linewidth: 2
+                   zorder: 3
+                 MultiModelPercentile10:
+                   color: gray
+                   label: '{dataset}'
+                   linestyle: '--'
+                   linewidth: 1
+                   zorder: 2
+                 MultiModelPercentile90:
+                   color: gray
+                   label: '{dataset}'
+                   linestyle: '--'
+                   linewidth: 1
+                   zorder: 2
+                 default:
+                   color: lightgray
+                   label: null
+                   linestyle: '-'
+                   linewidth: 1
+                   zorder: 1
+
+Variables
+---------
+
+Any, but the variables' number of dimensions should match the ones expected by each plot.
+
+References
+----------
+
+* Lauer, A., L. Bock, B. Hassler, P. Jöckel, L. Ruhe, and M. Schlund: Monitoring and benchmarking Earth System Model simulations with ESMValTool v2.12.0, Geosci. Model Dev., xx, xxxx-xxxx,
+  doi: xxx, 202x.
+
+Example plots
+-------------
+
+.. _fig_benchmarking_annual_cycle:
+.. figure::  /recipes/figures/benchmarking/annual_cycle.png
+   :align:   center
+   :width:   16cm
+
+(Left) Multi-year global mean (2000-2004) of the seasonal cycle of near-surface temperature in K from a simulation of MIROC6 and the reference dataset HadCRUT5 (black). The thin gray lines show individual CMIP6 models used for comparison, the dashed gray lines show the 10% and 90% percentiles of these CMIP6 models. (Right) same as (left) but for area-weighted RMSE of near-surface temperature. The light blue shading shows the range of the 10% to 90% percentiles of RMSE values from the ensemble of CMIP6 models used for comparison. Created with recipe_model_benchmarking_annual_cycle.yml.
+
+.. _fig_benchmarking_boxplots:
+.. figure::  /recipes/figures/benchmarking/boxplots.png
+   :align:   center
+   :width:   16cm
+
+(Left) Global area-weighted RMSE (smaller=better), (middle) weighted Pearson’s correlation coefficient (higher=better) and (right) weighted Earth mover’s distance (smaller=better) of the geographical pattern of 5-year means of different variables from a simulation of MIROC6 (red cross) in comparison to the CMIP6 ensemble (boxplot). Reference datasets for calculating the three metrics are: near-surface temperature (tas): HadCRUT5, surface temperature (ts): HadISST, precipitation (pr): GPCP-SG, air pressure at sea level (psl): ERA5, shortwave (rsut) longwave (rlut) radiative fluxes at TOA and shortwave (swcre) and longwave (lwcre) cloud radiative effects: CERES-EBAF. Each box indicates the range from the first quartile to the third quartile, the vertical lines show the median, and the whiskers the minimum and maximum values, excluding the outliers. Outliers are defined as being outside 1.5 times the interquartile range. Created with recipe_model_benchmarking_boxplots.yml.
+
+.. _fig_benchmarking_diurn_cycle:
+.. figure::  /recipes/figures/benchmarking/diurnal_cycle.png
+   :align:   center
+   :width:   10cm
+
+Area-weighted RMSE of the annual mean diurnal cycle (year 2000) of precipitation averaged over the tropical ocean (ocean grid cells in the latitude belt 30°S to 30°N) from a simulation of MIROC6 averaged compared with ERA5 data (black). The light blue shading shows the range of the 10% to 90% percentiles of RMSE values from the ensemble of CMIP6 models used for comparison. Created with recipe_benchmarking_diurnal_cycle.yml.
+
+.. _fig_benchmarking_map:
+.. figure::  /recipes/figures/benchmarking/map.png
+   :align:   center
+   :width:   10cm
+
+5-year annual mean (2000-2004) area-weighted RMSE of the precipitation rate in mm day-1 from a simulation of MIROC6 compared with GPCP-SG data. The stippled areas mask grid cells where the RMSE is smaller than the 90% percentile of RMSE values from an ensemble of CMIP6 models. Created with recipe_model_benchmarking_maps.yml
+
+.. _fig_benchmarking_timeseries:
+.. figure::  /recipes/figures/benchmarking/timeseries.png
+   :align:   center
+   :width:   16cm
+
+(Left) Time series from 2000 through 2014 of global average monthly mean temperature anomalies (reference period 2000-2009) of the near-surface temperature in K from a simulation of MIROC6 (red) and the reference dataset HadCRUT5 (black). The thin gray lines show individual CMIP6 models used for comparison, the dashed gray lines show the 10% and 90% percentiles of these CMIP6 models. (Right) same as (left) but for area-weighted RMSE of the near-surface air temperature. The light blue shading shows the range of the 10% to 90% percentiles of RMSE values from the ensemble of CMIP6 models used for comparison. Created with recipe_model_benchmarking_timeseries.yml.
+
+.. _fig_benchmarking_zonal:
+.. figure::  /recipes/figures/benchmarking/zonal.png
+   :align:   center
+   :width:   10cm
+
+5-year annual mean bias (2000-2004) of the zonally averaged temperature in K from a historical simulation of MIROC6 compared with ERA5 reanalysis data. The stippled areas mask grid cells where the absolute BIAS (${\abs{BIAS}}$) is smaller than the maximum of the absolute 10% (${\abs{p10}}$) and the absolute 90% (${\abs{p90}}$) percentiles from an ensemble of CMIP6 models, i.e. ${\abs{BIAS} \geq max( \abs{p10}, \abs{p90})}$. Created with recipe_model_benchmarking_zonal.yml.

doc/sphinx/source/recipes/recipe_model_evaluation.rst

@@ -96,3 +96,10 @@ Zonal mean precipitation.
    :width:   14cm
 
 Annual cycle of Southern Ocean total cloud cover.
+
+.. _fig_6:
+.. figure::  /recipes/figures/model_evaluation/diurnal_cycle_clt_sepacific_3hr.png
+   :align:   center
+   :width:   14cm
+
+Diurnal cycle of Southeast Pacific total cloud cover.

doc/sphinx/source/recipes/recipe_monitor.rst

@@ -189,6 +189,18 @@ Timeseries of tas including a reference dataset.
 
 Annual cycle of tas including a reference dataset.
 
+.. _fig_diurnal_cycle:
+.. figure::  /recipes/figures/monitor/diurnalcycle_pr_tropics_EC-Earth3_3hr_historical_r1i1p1f1.png
+   :align:   center
+   :width:   14cm
+
+.. _fig_diurnal_cycle_with_ref:
+.. figure::  /recipes/figures/monitor/diurnal_cycle_clt_tropics_3hr.png
+   :align:   center
+   :width:   14cm
+
+Diurnal cycle of clt including a reference dataset.
+
 .. _fig_map_with_ref:
 .. figure::  /recipes/figures/monitor/map_with_ref.png
    :align:   center

esmvaltool/config-references.yml

@@ -774,6 +774,7 @@ projects:
   crescendo: EU H2020 project CRESCENDO
   dlrveu2: DLR project VEU2
   dlrveu: DLR project VEU
+  dlrmabak: DLR project MABAK
   embrace: EU FP7 project EMBRACE
   esm2025: EU H2020 project ESM2025 - Earth system models for the future
   esmval: DLR project ESMVal

esmvaltool/diag_scripts/clouds/clouds.ncl

@@ -114,7 +114,8 @@ begin
 
   variables = metadata_att_as_array(variable_info, "short_name")
   if (.not. any(variables .eq. var0)) then
-    errstr = "diagnostic " + diag + " requires the following variable: " + var0
+    errstr = "diagnostic " + DIAG_SCRIPT \
+             + " requires the following variable: " + var0
     error_msg("f", DIAG_SCRIPT, "", errstr)
   end if
 
@@ -539,6 +540,10 @@ begin
       res@cnLevels            = ispan(0, 60, 5)
     end if
 
+    if (var0.eq."ts") then
+      res@cnLevels            = ispan(274, 304, 2)
+    end if
+
 ;    res@lbLabelBarOn         = False
     res@gsnRightString       = ""
 

esmvaltool/diag_scripts/monitor/monitor.py

@@ -27,6 +27,11 @@
       produce multi panel plots for data with `shape_id` or `region`
       coordinates of length > 1. Supported coordinates: `time`, `shape_id`
       (optional) and `region` (optional).
+    - Diurnal cycle (plot type ``diurnal_cycle``): Generate a diurnal cycle
+      plot (timeseries like climatological from 0 to 24 hours). It will
+      produce multi panel plots for data with `shape_id` or `region`
+      coordinates of length > 1. Supported coordinates: `time`, `shape_id`
+      (optional) and `region` (optional).
 
 Configuration options in recipe
 -------------------------------
@@ -39,10 +44,10 @@
     monitor configuration file can be found :ref:`here <monitor_config_file>`.
 plots: dict, optional
     Plot types plotted by this diagnostic (see list above). Dictionary keys
-    must be ``clim``, ``seasonclim``, ``monclim``, ``timeseries`` or
-    ``annual_cycle``. Dictionary values are dictionaries used as options for
-    the corresponding plot. The allowed options for the different plot types
-    are given below.
+    must be ``clim``, ``seasonclim``, ``monclim``, ``timeseries``,
+    ``annual_cycle`` or ``diurnal_cycle``. Dictionary values are dictionaries
+    used as options for the corresponding plot. The allowed options for the
+    different plot types are given below.
 plot_filename: str, optional
     Filename pattern for the plots.
     Defaults to ``{plot_type}_{real_name}_{dataset}_{mip}_{exp}_{ensemble}``.
@@ -98,6 +103,10 @@
 ----------------------------------------------------
 None
 
+Configuration options for plot type ``diurnal_cycle``
+-----------------------------------------------------
+None
+
 .. hint::
 
    Extra arguments given to the recipe are ignored, so it is safe to use yaml
@@ -166,6 +175,7 @@ def compute(self):
 
                 self.timeseries(cube, var_info)
                 self.plot_annual_cycle(cube, var_info)
+                self.plot_diurnal_cycle(cube, var_info)
                 self.plot_monthly_climatology(cube, var_info)
                 self.plot_seasonal_climatology(cube, var_info)
                 self.plot_climatology(cube, var_info)
@@ -280,6 +290,57 @@ def plot_annual_cycle(self, cube, var_info):
             caption=caption,
         )
 
+    def plot_diurnal_cycle(self, cube, var_info):
+        """Plot the diurnal cycle according to configuration.
+
+        The key 'diurnal_cycle' must be passed to the 'plots' option in the
+        configuration.
+
+        Parameters
+        ----------
+        cube: iris.cube.Cube
+            Data to plot. Must be 1D with time or 2D with an extra 'shape_id'
+            or 'region' coordinate. In that case, the plot will be a multiple
+            one with one figure for each region
+        var_info: dict
+            Variable's metadata from ESMValTool
+
+        Warning
+        -------
+        The hourly climatology is done inside the function so the users can
+        plot both the timeseries and the diurnal cycle in one go
+        """
+        if 'diurnal_cycle' not in self.plots:
+            return
+        cube = climate_statistics(cube, period='hour')
+
+        plotter = PlotSeries()
+        plotter.outdir = self.get_plot_folder(var_info)
+        plotter.img_template = self.get_plot_path('diurnalcycle', var_info,
+                                                  add_ext=False)
+        plotter.filefmt = self.cfg['output_file_type']
+        region_coords = ('shape_id', 'region')
+        options = {
+            'xlabel': '',
+            'xlimits': None,
+            'suptitle': 'Diurnal cycle',
+        }
+        for region_coord in region_coords:
+            if cube.coords(region_coord):
+                plotter.multiplot_cube(cube, 'month', region_coord, **options)
+                return
+        plotter.plot_cube(cube, 'hour', **options)
+        caption = (f"Diurnal cycle of {var_info[n.LONG_NAME]} of "
+                   f"dataset {var_info[n.DATASET]} (project "
+                   f"{var_info[n.PROJECT]}) from {var_info[n.START_YEAR]} to "
+                   f"{var_info[n.END_YEAR]}.")
+        self.record_plot_provenance(
+            self.get_plot_path('diurnalcycle', var_info),
+            var_info,
+            'Diurnal cycle',
+            caption=caption,
+        )
+
     def plot_monthly_climatology(self, cube, var_info):
         """Plot the monthly climatology as a multipanel plot.