Add SST bias diagnostic for IPCC AR6 Ch3 Fig. 3.24. Can be used for REF.

[malininae]

• Closes SST bias diagnostic #3943

---

Before you get started

• ☝ Create an issue to discuss what you are going to do

Checklist

It is the responsibility of the author to make sure the pull request is ready to review. The icons indicate whether the item will be subject to the 🛠 Technical or 🧪 Scientific review.

• 🛠 This pull request has a descriptive title
• 🛠 Code is written according to the code quality guidelines
• 🛠 Documentation is available
• 🛠 Tests run successfully
• 🛠 The list of authors is up to date
• 🛠 All checks below this pull request were successful

New or updated recipe/diagnostic

• 🧪 Recipe runs successfully
• 🧪 Recipe is well documented
• 🧪 Figure(s) and data look as expected from literature
• 🛠 Provenance information has been added

---

[malininae]

@bouweandela I made the diagnostic as it was done for IPCC, but I realized that for REF the requirement might be different. It's my first PR for a diagnostic for the main branch of ESMValTool, so I am happy to adjust and can do it quickly.

[malininae]

If someone would be able to run the recipe on DKRZ I would appreciate it. I think I've been removed from the project and I'm not sure it would be able to do it and we don't have a full CMIP5 suit on the local computers.

If you run the tool with the command esmvaltool run ipccwg1ar6ch3/recipe_ipccwg1ar6ch3_fig_3_24.yml --search-esgf=always it will automatically download all the required CMIP5 data.

Thanks. CESM1-CAM5 and CESM1-WACCM are missing on ESGF, but available on DKRZ. I've copied these two from there, but the question is, what to do with them long-term, since as far as I understand CMIP5 is being discontinued soon.

[bouweandela]

I made the diagnostic as it was done for IPCC, but I realized that for REF the requirement might be different.

For the REF we will only be using CMIP7 Fast Track data (and obs4MIPs) initially, we're currently testing with CMIP6 data because CMIP7 data is not yet available. More data, such as high resolution, may be added at a later stage. Indeed it seems unlikely that CMIP5 data will be used. However, we do update the datasets section of recipes when running them with the REF, to make sure new data is processed as it comes in, so it does not matter that the CMIP5 data is there, as long as the recipe runs fine without it. I tried running the recipe with only the cmip6 and eq_mask variables and that seems to work fine, can you look at the plots to check if they look OK if you do that?

I did notice that the captions could be more clear, e.g. what is on the axes? which datasets are used? what does the line and shaded area represent?

[ledm]

This PR is ready for review. The recipe is under IPCC AR6 Ch3, but I made in on time for REF and can be used for SST bias diagnostic. @ledm I completely rewrote you original diagnostic, but I left your original development in the history. I would really appreciate if you could do a code review, because you did the original one.

@bouweandela I update the directory with your new PR with ruff and precommit. I get fails on CircleCI, but I don't understand what it wants from me with the docstrings. I tried all the possible combinations (removing empty strings in the docs, after the docs, moving the """ around etc.) and nothing works. According to documentation on the internet my docstrings answer all the requirements, so I'm a bit confused.

Hi @malininae, I'm pretty far removed from the initial code at this stage - the last time I had sight on it would have been 2021, and there's been a lot of development since I last used ESMValTool - not sure I'm the person you want for the code review today! (The plots are pretty though!)

[malininae]

I made the diagnostic as it was done for IPCC, but I realized that for REF the requirement might be different.

For the REF we will only be using CMIP7 Fast Track data (and obs4MIPs) initially, we're currently testing with CMIP6 data because CMIP7 data is not yet available. More data, such as high resolution, may be added at a later stage. Indeed it seems unlikely that CMIP5 data will be used. However, we do update the datasets section of recipes when running them with the REF, to make sure new data is processed as it comes in, so it does not matter that the CMIP5 data is there, as long as the recipe runs fine without it. I tried running the recipe with only the cmip6 and eq_mask variables and that seems to work fine, can you look at the plots to check if they look OK if you do that?

I did notice that the captions could be more clear, e.g. what is on the axes? which datasets are used? what does the line and shaded area represent?

Thanks! I added the full description to the recipe, and saving the data for the plot into the nc.-files. Let me know if needs to be expanded.

Yes, CMIP5 data is not needed there, I tried to make it as much as possible do one can plug and play with whatever variable groups one wants. The eq_mask or whatever other resolved mask is needed only if the mask type is resolved.

I had to remove the MRI models from CMIP5 because they were super weird and had the craziest values. In the IPCC version they were filtered through

        if new_cube.data.mean() < -200. :
                print ("Warning: this model is broken:", dataset, 'mean:', new_cube.data.mean())
                continue

So they weren't used-used in the plot.

[bouweandela]

I tried running the recipe the way I would try running it for the REF, with just one CMIP6 dataset, but it crashes with the following error from script ocean/diagnostic_biases.py:

INFO:diagnostic_biases.py:There were no statistics calculated for cmip6because only one cube was provided
There were no statistics calculated for cmip6because only one cube was provided
Traceback (most recent call last):
  File "/home/b/b381141/src/esmvalgroup/esmvaltool/esmvaltool/diag_scripts/ocean/diagnostic_biases.py", line 383, in <module>
    main(config)
  File "/home/b/b381141/src/esmvalgroup/esmvaltool/esmvaltool/diag_scripts/ocean/diagnostic_biases.py", line 376, in main
    plot_bias_plot(data_list, cfg)
  File "/home/b/b381141/src/esmvalgroup/esmvaltool/esmvaltool/diag_scripts/ocean/diagnostic_biases.py", line 324, in plot_bias_plot
    iris.plot.fill_between(
  File "/home/b/b381141/src/scitools/iris/lib/iris/plot.py", line 1717, in fill_between
    raise TypeError("y1 must be a cube or a coordinate.")
TypeError: y1 must be a cube or a coordinate.

Recipe:

# ESMValTool
# recipe_ipccwg1ar6ch3_fig_3_24.yml
---
documentation:

  title: Zonal and equatorial SST bias

  description: Recipe to produce figure 3.24 for AR6 IPCC WG1.
    Biases in zonal mean and equatorial sea surface temperature (SST)
    in CMIP5, CMIP6 and HighResMIP models. Panel (a) shows multi-model
    mean zonally averaged SST bias; (b) equatorial SST bias;
    and (c) equatorial SST compared to observed mean SST (black line)
    for 1979–1999. The inter-model 5th and 95th percentiles are
    depicted by the respective shaded range. Model climatologies are
    derived from the 1979–1999 mean of the historical simulations,
    using one simulation per model. The Hadley Centre Sea Ice and
    Sea Surface Temperature version 1 (HadISST) observational
    climatology for 1979–1999 is used as the reference for the error
    calculation in (a) and (b); and for observations in (c).

  authors:
    - malinina_elizaveta
    - demora_lee

  maintainer:
    - malinina_elizaveta

  references:
    - eyring21ipcc

  projects:
    - ipcc_ar6

models_cmip6: &models_cmip6
  - {dataset: ACCESS-ESM1-5 , project: CMIP6, mip: Omon, exp: historical, ensemble: r1i1p1f1, grid: gn}
  # observations/reference dataset
  - {dataset: HadISST, project: OBS, mip: Omon, type: reanaly, version: 1, tier: 2}

preprocessors:
  prep_zonal:
    custom_order: true
    climate_statistics:
      operator: mean
      period: full
    regrid:
      target_grid: 1x1
      scheme: linear
    zonal_statistics:
      operator: mean
    convert_units:
      units: degrees_C

  prep_equatorial:
    custom_order: true
    climate_statistics:
      operator: mean
      period: full
    regrid:
      target_grid: 1x1
      scheme: linear
    extract_region:
      start_longitude: 0.
      end_longitude: 360.
      start_latitude: -10.
      end_latitude: 10.
    meridional_statistics:
      operator: mean
    convert_units:
      units: degrees_C

  prep_eq_mask:
    custom_order: true
    climate_statistics:
      operator: mean
      period: full
    regrid:
      target_grid: 1x1
      scheme: linear
    extract_region:
      start_longitude: 0.
      end_longitude: 360.
      start_latitude: -10.
      end_latitude: 10.
    convert_units:
      units: degrees_C

diagnostics:
  zonal:
    variables:
      cmip6:
        short_name: tos
        preprocessor: prep_zonal
        start_year: 1979
        end_year: 1999
        reference_dataset: HadISST
        additional_datasets: *models_cmip6
    scripts:
      panel_a:
        script: ocean/diagnostic_biases.py
        bias: True # if bias should be calculated (optional), default: false
        mask: # if the data should be masked (optional), default: false
          flag: true
          type: simple # accepted values: 'simple' and 'resolved'
        data_statistics:
          best_guess: mean # best guess statistic
          borders: # borders for the shading
            - operator: percentile
              percent: 5
            - operator: percentile
              percent: 95
        mpl_style: ipcc_ar6_fgd # name of the matplotlib style file (optional)
        caption: (a) Zonal mean SST bias # Figure caption (optional)
        color_style: sst_bias # name of the color style, colors are defined for variable groups

  equatorial:
    variables:
      cmip6:
        short_name: tos
        preprocessor: prep_equatorial
        start_year: 1979
        end_year: 1999
        reference_dataset: HadISST
        additional_datasets: *models_cmip6
      eq_mask:
        short_name: tos
        preprocessor: prep_eq_mask
        start_year: 1979
        end_year: 1999
        additional_datasets:
         - {dataset: HadISST, project: OBS, mip: Omon, type: reanaly, version: 1, tier: 2}

    scripts:
      panel_b:
        script: ocean/diagnostic_biases.py
        bias: True
        mask:
          flag: true
          type: resolved
          group: eq_mask
        data_statistics:
          best_guess: mean
          borders:
            - operator: percentile
              percent: 5
            - operator: percentile
              percent: 95
        mpl_style: ipcc_ar6_fgd
        caption: (b) Equatorial SST bias
        color_style: sst_bias

      panel_c:
        script: ocean/diagnostic_biases.py
        bias: False
        mask:
          flag: true
          type: resolved
          group: eq_mask
        data_statistics:
          best_guess: mean
          borders:
            - operator: percentile
              percent: 5
            - operator: percentile
              percent: 95
        mpl_style: ipcc_ar6_fgd
        caption: (c) Equatorial SST
        color_style: sst_bias

Or does it only make sense to run with all datasets?

[malininae]

I tried running the recipe the way I would try running it for the REF, with just one CMIP6 dataset, but it crashes with the following error from script ocean/diagnostic_biases.py:

INFO:diagnostic_biases.py:There were no statistics calculated for cmip6because only one cube was provided
There were no statistics calculated for cmip6because only one cube was provided
Traceback (most recent call last):
  File "/home/b/b381141/src/esmvalgroup/esmvaltool/esmvaltool/diag_scripts/ocean/diagnostic_biases.py", line 383, in <module>
    main(config)
  File "/home/b/b381141/src/esmvalgroup/esmvaltool/esmvaltool/diag_scripts/ocean/diagnostic_biases.py", line 376, in main
    plot_bias_plot(data_list, cfg)
  File "/home/b/b381141/src/esmvalgroup/esmvaltool/esmvaltool/diag_scripts/ocean/diagnostic_biases.py", line 324, in plot_bias_plot
    iris.plot.fill_between(
  File "/home/b/b381141/src/scitools/iris/lib/iris/plot.py", line 1717, in fill_between
    raise TypeError("y1 must be a cube or a coordinate.")
TypeError: y1 must be a cube or a coordinate.

Recipe:

Or does it only make sense to run with all datasets?

Good catch! I thought Nones would do it. I've added a check, if borders exist.

[bouweandela]

Thanks! Could you also save the data in the figures as NetCDF so we could potentially use it to create a different plot or allow users to select which models/observations they want plotted?

I noticed that there is no legend for the figures, so it is a bit unclear what the different lines and shaded areas represent. Could you add that too?

[bouweandela]

Thanks, looking good!