diff --git a/doc/sphinx/source/recipes/figures/ocean_south/barotropic_streamfunction_ACCESS-OM2.png b/doc/sphinx/source/recipes/figures/ocean_south/barotropic_streamfunction_ACCESS-OM2.png
new file mode 100644
index 0000000000..0aa085a9e3
Binary files /dev/null and b/doc/sphinx/source/recipes/figures/ocean_south/barotropic_streamfunction_ACCESS-OM2.png differ
diff --git a/doc/sphinx/source/recipes/figures/ocean_south/hovmoller_ACCESS-OM2.png b/doc/sphinx/source/recipes/figures/ocean_south/hovmoller_ACCESS-OM2.png
new file mode 100644
index 0000000000..d89a00e11a
Binary files /dev/null and b/doc/sphinx/source/recipes/figures/ocean_south/hovmoller_ACCESS-OM2.png differ
diff --git a/doc/sphinx/source/recipes/figures/ocean_south/kinetic_energy_ACCESS-OM2-025.png b/doc/sphinx/source/recipes/figures/ocean_south/kinetic_energy_ACCESS-OM2-025.png
new file mode 100644
index 0000000000..02fec495d7
Binary files /dev/null and b/doc/sphinx/source/recipes/figures/ocean_south/kinetic_energy_ACCESS-OM2-025.png differ
diff --git a/doc/sphinx/source/recipes/recipe_ocean_south.rst b/doc/sphinx/source/recipes/recipe_ocean_south.rst
new file mode 100644
index 0000000000..d9bd876df7
--- /dev/null
+++ b/doc/sphinx/source/recipes/recipe_ocean_south.rst
@@ -0,0 +1,70 @@
+.. _recipes_ocean_south:
+
+Simple southern ocean diagnostics
+=====
+
+Overview
+--------
+
+A few simple recipes recreated in ESMValTool for running on bulk datasets e.g. CMIP6, from the `COSIMA cookbook <https://cosima-recipes.readthedocs.io/en/latest/index.html>`_.
+
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in esmvaltool/recipes/
+
+* recipe_ocean_south.yml
+
+Diagnostics are stored in esmvaltool/diag_scripts/ocean_south/
+
+* plot_barotropic.py: for plotting barotropic streamfunction
+* plot_hovmoller.py: for plotting Hovmoller diagrams of ocean temperature and salinity
+* plot_kinetic_energy.py: for plotting ocean kinetic energy
+
+
+User settings in recipe
+-----------------------
+
+*None*
+
+
+Variables
+---------
+
+* umo (ocean, monthly mean, longitude latitude time)
+* so (ocean, monthly mean, longitude latitude time depth)
+* thetao (ocean, monthly mean, longitude latitude time depth)
+* uo (ocean, monthly mean, longitude latitude time depth)
+* vo (ocean, monthly mean, longitude latitude time depth)
+
+
+
+References
+----------
+
+* https://cosima-recipes.readthedocs.io/en/latest/02-Appetisers/Barotropic_Streamfunction.html
+* https://cosima-recipes.readthedocs.io/en/latest/02-Appetisers/Hovmoller_Temperature_Depth.html
+* https://cosima-recipes.readthedocs.io/en/latest/03-Mains/Eddy-Mean_Kinetic_Energy_Decomposition.html
+
+
+Example plots
+-------------
+
+.. _fig_1:
+.. figure::  /recipes/figures/ocean_south/hovmoller_ACCESS-OM2.png
+   :align:   center
+
+   Hovmoller diagram of ocean temperature and salinity for the ACCESS-OM2 model.
+
+.. _fig_2:
+.. figure::  /recipes/figures/ocean_south/kinetic_energy_ACCESS-OM2-025.png
+   :align:   center
+
+   Kinetic energy around Antarctica for the ACCESS-OM2-025 model.
+
+.. _fig_3:
+.. figure::  /recipes/figures/ocean_south/barotropic_streamfunction_ACCESS-OM2.png
+   :align:   center
+
+   Barotropic streamfunction around Antarctica for the ACCESS-OM2 model.
diff --git a/esmvaltool/diag_scripts/ocean_south/plot_barotropic.py b/esmvaltool/diag_scripts/ocean_south/plot_barotropic.py
new file mode 100755
index 0000000000..6b00a695c5
--- /dev/null
+++ b/esmvaltool/diag_scripts/ocean_south/plot_barotropic.py
@@ -0,0 +1,110 @@
+"""diagnostic script to plot"""
+
+import iris
+import os
+import logging
+import numpy as np
+import iris.plot as iplt
+import matplotlib.path as mpath
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import cmocean
+from esmvaltool.diag_scripts.shared import (
+    run_diagnostic,
+    save_figure,
+    group_metadata,
+    select_metadata,
+)
+from esmvalcore.preprocessor import (
+    climate_statistics,
+    regrid
+)
+
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def circumpolar_map():
+    fig = plt.figure(figsize = (12, 8))
+    ax = plt.axes(projection = ccrs.SouthPolarStereo())
+    ax.set_extent([-180, 180, -80, -40], crs = ccrs.PlateCarree())
+    ax.set_facecolor('lightgrey')
+    # Map the plot boundaries to a circle
+    theta = np.linspace(0, 2 * np.pi, 100)
+    center, radius = [0.5, 0.5], 0.5
+    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
+    circle = mpath.Path(verts * radius + center)
+    ax.set_boundary(circle, transform = ax.transAxes)
+
+    return fig, ax
+
+
+def get_provenance_record(caption, ancestor_files):
+    """Create a provenance record describing the diagnostic data and plot."""
+
+    record = {
+        "caption": caption,
+        "statistics": ["other"],
+        "domains": ["polar"],
+        "plot_types": ["map"],
+        "authors": [
+            "chun_felicity",
+        ],
+        "references": [
+            "access-nri",
+        ],
+        "ancestors": ancestor_files,
+    }
+    return record
+
+
+def main(cfg):
+    """# in script cumsum, units handling, climate mean"""
+
+    input_data = cfg["input_data"].values()
+
+    for dataset in input_data:
+        # Load the data
+        input_file = dataset['filename']
+        name = dataset['dataset']
+        cube = iris.load_cube(input_file)
+
+        # units cumulative sum and climate mean
+        p0 = 1035.0
+        cube = cube.collapsed('depth', iris.analysis.SUM)
+        cube.data = cube.data / p0  # divide by density for volume
+        cube.data = cube.data / 1e6  # then divide by 10^6 for Sv
+        cube = regrid(cube, target_grid='0.5x0.5', scheme="linear")
+        cube.data = cube.data.cumsum(1)  #latitude index 1
+        logger.info(cube.summary(shorten=True))
+        cube.units = 'Sv'
+        cube = climate_statistics(cube, period='full', operator='mean')
+
+        # plot
+        fig, ax = circumpolar_map()
+
+        p1 = iplt.contourf(cube,
+                    levels = np.arange(-50, 150, 10),
+                    extend = 'both',
+                    cmap = cmocean.cm.speed,
+                    )
+        plt.colorbar(p1, orientation='vertical', label='$\psi$ (Sv)')
+
+        plt.title(f'Barotropic streamfunction from {name}')
+        # Save output
+        prov_record = get_provenance_record(
+            f'Barotropic streamfunction from {name}.',
+            [input_file],
+        )
+        save_figure(
+            f"barotropic_streamfunction_{name}",
+            prov_record,
+            cfg,
+            figure=fig,
+            dpi=300,
+            )
+
+if __name__ == "__main__":
+    with run_diagnostic() as config:
+        main(config)
diff --git a/esmvaltool/diag_scripts/ocean_south/plot_hovmoller.py b/esmvaltool/diag_scripts/ocean_south/plot_hovmoller.py
new file mode 100755
index 0000000000..c665505f2e
--- /dev/null
+++ b/esmvaltool/diag_scripts/ocean_south/plot_hovmoller.py
@@ -0,0 +1,177 @@
+"""diagnostic script to plot"""
+
+import iris
+import os
+import logging
+import numpy as np
+import iris.plot as iplt
+from matplotlib.gridspec import GridSpec
+from matplotlib import ticker
+import matplotlib.pyplot as plt
+
+import cmocean
+from esmvaltool.diag_scripts.shared import (
+    run_diagnostic,
+    save_figure,
+    group_metadata,
+    select_metadata,
+)
+from esmvalcore.preprocessor import (
+    area_statistics,
+    anomalies
+)
+
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def plot_hovmoller(fsize = 14):
+    """Setup of figures properties."""
+    plt.rcParams['font.size'] = fsize
+    plt.rcParams['xtick.labelsize'] = fsize-2
+    plt.rcParams['ytick.labelsize'] = fsize-2
+    
+    fig = plt.figure(figsize = (8, 6))
+    grid = GridSpec(100, 100)
+    
+    ax = [fig.add_subplot(grid[:30, :45]), fig.add_subplot(grid[:30, 48:93]),
+          fig.add_subplot(grid[32:, :45]), fig.add_subplot(grid[32:, 48:93])]
+    
+    for i in range(len(ax)):
+        # ax[i].xaxis.set_major_formatter(date_format)
+        ax[i].tick_params(axis='x', labelrotation=45)
+        if i != 0 and i != 2:
+            ax[i].set_yticklabels([]) 
+
+    return fig, ax
+
+
+def hovmoller(ds_cube, startyr):
+    """Calculate anomalies summed over area for depth vs time plot."""
+    cube = anomalies(ds_cube, period='full', reference={'start_year':startyr, 'start_month':1,'start_day':1,
+                                                            'end_year':startyr, 'end_month':12, 'end_day':31})
+    # anomalies doesn't retain cell measure required for area statistics
+    cube.add_cell_measure(ds_cube.cell_measure('cell_area'), data_dims=[2,3])
+    # sum anomalies
+    cube = area_statistics(cube, operator='sum')
+
+    #divide by total area cube (anomalies are summed.. weighted by cell_area(multiplied)
+    tarea = totalarea_depth(ds_cube)
+    return cube / tarea
+
+
+def totalarea_depth(cube):
+    """Calculate total area at each depth level."""
+    tmp1 = cube[0] #get first time step
+    # divide by itself for 1s mask 3-dimensional
+    tmp1.data = tmp1.data / tmp1.data
+    # check cell_measure exists?
+    if cube.cell_measure('cell_area'):
+        return area_statistics(tmp1, operator='sum')
+
+
+def plot_details(ax, cf_temp, cf_salt):
+    for i in range(len(ax)):
+        if i < 2:
+            ax[i].set_ylim(500, 0)
+            ax[i].set_xticklabels([])
+        else:
+            ax[i].set_xlabel("")
+            ax[i].set_ylim(5000, 500)
+    ax[0].set_ylabel("Depth [m]")
+    ax[1].set_ylabel("")
+    ax[2].set_ylabel("Depth [m]")
+    ax[3].set_ylabel("")
+
+    # Colorbars
+    bar = plt.axes([0.11, 0.97, 0.35, 0.02])
+    cbar_1 = plt.colorbar(cf_temp, cax = bar, orientation = 'horizontal', extend='both', format= '%.2f')  
+    cbar_1.set_label("Temperature [$\degree$C]")
+
+    bar = plt.axes([0.50, 0.97, 0.35, 0.02])
+    cbar_2 = plt.colorbar(cf_salt, cax = bar, orientation = 'horizontal', extend='both', format= '%.2f')  
+    cbar_2.set_label("Salinity [psu]")
+
+    for cbar in [cbar_1, cbar_2]:
+        tick_locator = ticker.MaxNLocator(nbins=3, prune='both') ## The ticker needs to called within the loop
+        cbar.locator = tick_locator
+        cbar.update_ticks()
+
+def get_provenance_record(caption, ancestor_files):
+    """Create a provenance record describing the diagnostic data and plot."""
+
+    record = {
+        "caption": caption,
+        "statistics": ["other"],
+        "domains": ["global"],
+        "plot_types": ["vert"],
+        "authors": [
+            "chun_felicity",
+        ],
+        "references": [
+            "access-nri",
+        ],
+        "ancestors": ancestor_files,
+    }
+    return record
+
+
+def main(cfg):
+    """Create Hovmoller diagrams."""
+
+    input_data = cfg["input_data"].values()
+
+    # group by dataset
+    ds_groups = group_metadata(
+        input_data,
+        "dataset",
+    )
+
+    levels_tempsal = {'thetao':np.arange(-0.3, 0.31, 0.01), 'so':np.arange(-0.03, 0.031, 0.001)}
+    cmapcm = {'thetao':cmocean.cm.balance, 'so':cmocean.cm.curl}
+    axi = {'thetao':0, 'so':1}
+    # for each ds plt ts, sal 
+    for grp, var_attr in ds_groups.items():
+        fig, ax = plot_hovmoller(fsize = 14)
+        files=[]
+        for metadata in var_attr:  #[thetao, so]
+            name = metadata["dataset"]
+            shortname = metadata["short_name"]
+            input_file = metadata['filename']
+            cube = iris.load_cube(metadata['filename'])
+            hov_cube = hovmoller(cube, metadata['start_year'])
+
+            # plot twice
+            for i in [axi[shortname], axi[shortname]+2]:
+                p1 = iplt.contourf(hov_cube,
+                            levels = levels_tempsal[shortname],
+                            extend = 'both',
+                            cmap = cmapcm[shortname],
+                            axes = ax[i], #0,2 temp
+                            )
+            files.append(input_file)
+
+            if shortname == 'so':
+                cf_salt = p1
+            else:
+                cf_temp = p1
+        
+        plot_details(ax, cf_temp, cf_salt)
+        # Save output
+        prov_record = get_provenance_record(
+            f'Depth-Time Temperature and Salinity from {name}.',
+            files,
+        )
+        save_figure(
+            f"hovmoller_{name}",
+            prov_record,
+            cfg,
+            figure=fig,
+            dpi=300,
+            )
+
+
+if __name__ == "__main__":
+    with run_diagnostic() as config:
+        main(config)
diff --git a/esmvaltool/diag_scripts/ocean_south/plot_kinetic_energy.py b/esmvaltool/diag_scripts/ocean_south/plot_kinetic_energy.py
new file mode 100755
index 0000000000..96af49befe
--- /dev/null
+++ b/esmvaltool/diag_scripts/ocean_south/plot_kinetic_energy.py
@@ -0,0 +1,152 @@
+"""diagnostic script to plot"""
+
+import iris
+import os
+import logging
+import numpy as np
+import iris.plot as iplt
+import matplotlib.pyplot as plt
+
+import cmocean
+import cartopy.crs as ccrs
+import cartopy.feature as feature
+import matplotlib.path as mpath
+from esmvaltool.diag_scripts.shared import (
+    run_diagnostic,
+    save_figure,
+    group_metadata,
+    save_data,
+)
+from esmvalcore.preprocessor import (
+    add_supplementary_variables,
+    axis_statistics,
+    regrid,
+)
+
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def ke_calc(u_vel, v_vel, thkcel):
+    KE = 0.5*(u_vel**2 + v_vel**2)
+    KE = add_supplementary_variables(KE, [thkcel])
+    return KE
+
+def kinetic_energys(u_vel, v_vel, thkcel):
+    """Calculate kinetic energy from u and v components."""
+    
+    u_mean = axis_statistics(u_vel, operator='mean', axis='t')
+    v_mean = axis_statistics(v_vel, operator='mean', axis='t')
+    mke = ke_calc(u_mean, v_mean, thkcel)
+    mke = mke.collapsed('depth', iris.analysis.SUM, weights='cell_thickness')
+    mke.long_name = 'Mean Kinetic Energy'
+    # eke
+    u_transient = u_vel - u_mean
+    v_transient = v_vel - v_mean
+    eke = ke_calc(u_transient, v_transient, thkcel) 
+    eke = eke.collapsed('depth', iris.analysis.SUM, weights='cell_thickness')
+    eke = axis_statistics(eke, operator='mean', axis='t')
+    eke.long_name = 'Eddy Kinetic Energy'
+    # tke
+    tke = ke_calc(u_vel, v_vel, thkcel)
+    tke = axis_statistics(tke, operator='mean', axis='t')
+    tke = tke.collapsed('depth', iris.analysis.SUM, weights='cell_thickness')
+    tke.long_name = 'Total Kinetic Energy'
+    
+    KEall = [regrid(ke, '0.5x0.5', scheme='nearest') for ke in [tke, mke, eke]]  # regrid to map
+    return KEall # tke,mke,eke
+
+def circumpolar_map(fig, i): # make 3 subplots
+    # fig = plt.figure(figsize = (12, 8)) #move
+    # ax = plt.axes(projection = ccrs.SouthPolarStereo())
+    ax = fig.add_subplot(i, projection = ccrs.SouthPolarStereo()) # i 221
+    ax.set_extent([-180, 180, -80, -50], crs = ccrs.PlateCarree())
+    ax.set_facecolor('lightgrey')
+    # Map the plot boundaries to a circle
+    theta = np.linspace(0, 2 * np.pi, 100)
+    center, radius = [0.5, 0.5], 0.5
+    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
+    circle = mpath.Path(verts * radius + center)
+    ax.set_boundary(circle, transform = ax.transAxes)
+
+    return ax
+
+
+def get_provenance_record(caption, ancestor_files):
+    """Create a provenance record describing the diagnostic data and plot."""
+
+    record = {
+        "caption": caption,
+        "statistics": ["other"],
+        "domains": ["shpolar"],
+        "plot_types": ["map"],
+        "authors": [
+            "chun_felicity",
+        ],
+        "references": [
+            "access-nri",
+        ],
+        "ancestors": ancestor_files,
+    }
+    return record
+
+
+def main(cfg):
+    """Create kinetic energy diagrams."""
+
+    input_data = cfg["input_data"].values()
+
+    # group by dataset
+    ds_groups = group_metadata(
+        input_data,
+        "dataset",
+    )
+
+    #
+    for grp, var_attr in ds_groups.items():
+        fig = plt.figure(figsize = (12, 10))
+        
+        files, data_cubes = [], {}
+        for metadata in var_attr:  #{vo,uo, thkcello}
+            
+            shortname = metadata["short_name"]
+            input_file = metadata['filename']
+            cube = iris.load_cube(metadata['filename'])
+            data_cubes[shortname] = cube
+            files.append(input_file)
+            thkness = cube.ancillary_variable('cell_thickness')
+
+        tke, mke, eke = kinetic_energys(data_cubes['uo'], data_cubes['vo'], thkness)
+
+        plt_titles = ['Total Kinetic Energy (TKE)', 'Mean Kinetic Energy (MKE)', 'Eddy Kinetic Energy (EKE)']
+        # plot TKE, MKE, EKE (3subplts)
+        for i, ke_cube in enumerate([tke, mke, eke]):
+            axs = circumpolar_map(fig, 221+i)
+            cf1 = iplt.contourf(ke_cube, levels=np.arange(0,26,0.5), extend='max', axes=axs, cmap = cmocean.cm.ice)
+            axs.set_title(plt_titles[i])
+
+            data_prov = get_provenance_record(
+                f'{plt_titles[i]} for {grp}.',
+                files,
+            ) # save data out
+            save_data(f"{plt_titles[i][-4:-1]}_{grp}", data_prov, cfg, ke_cube)
+        
+        fig.colorbar(cf1, cax=fig.add_axes([0.85, 0.08, 0.02, 0.4]), label='m$^3$ s$^{-2}$',ticks=np.arange(0,26,5), shrink=0.6) # share colorbar
+        # Save figure
+        prov_record = get_provenance_record(
+            f'Kinetic Energy {grp}.',
+            files,
+        )
+        save_figure(
+            f"kinetic_energy_{grp}",
+            prov_record,
+            cfg,
+            figure=fig,
+            dpi=300,
+            )
+
+
+if __name__ == "__main__":
+    with run_diagnostic() as config:
+        main(config)
diff --git a/esmvaltool/recipes/recipe_ocean_south.yml b/esmvaltool/recipes/recipe_ocean_south.yml
new file mode 100644
index 0000000000..e057ea9932
--- /dev/null
+++ b/esmvaltool/recipes/recipe_ocean_south.yml
@@ -0,0 +1,66 @@
+# ESMValTool
+# recipe_ocean_south.yml
+---
+documentation:
+  description: recreating simple cosima recipes for bulk datasets
+  title: Recreating simple cosima recipes
+  authors:
+    - chun_felicity
+  maintainer:
+    - chun_felicity
+  
+  
+diagnostics:
+
+  barotropic_streamfunction:
+    description:
+     Calculate barotropic streamfunction in the Southern Ocean
+    variables:
+      sst_east:
+        short_name: umo
+        mip: Omon
+        additional_datasets: #additional datasets can be added for each diagnostic
+          - {dataset: ACCESS-OM2, project: CMIP6, exp: omip2, activity: OMIP,
+              ensemble: r1i1p1f1, grid: gn, start_year: 0356, end_year: 0366}
+          - {dataset: ACCESS-OM2-025, project: CMIP6, exp: omip2, activity: OMIP,
+              ensemble: r1i1p1f1, grid: gn, start_year: 0356, end_year: 0366}
+    scripts:
+      plot_script:
+        script: ocean_south/plot_barotropic.py
+
+  hovmoller_diagnostic:
+    description: Hovmoller diagram of temperature and salininty anomalies
+    variables:
+      so:
+        mip: Omon
+      thetao:
+        mip: Omon
+    additional_datasets:
+      - {dataset: ACCESS-ESM1-5, project: CMIP6, exp: historical,
+          ensemble: r1i1p1f1, grid: gn, start_year: 1950, end_year: 2014,
+          supplementary_variables: [{short_name: areacello, mip: Ofx}]}
+      - {dataset: ACCESS-OM2, project: CMIP6, exp: omip2, activity: OMIP,
+        ensemble: r1i1p1f1, grid: gn, timerange: 0336/0356,
+        supplementary_variables: [{short_name: areacello, mip: Ofx}]}
+    scripts:
+      plot_script:
+        script: ocean_south/plot_hovmoller.py
+
+  kinetic_energy:
+    description: |
+      Decomposing the kinetic energy into time-mean and transient components.
+    variables:
+      uo:
+        mip: Omon
+      vo:
+        mip: Omon
+    additional_datasets:
+      - {dataset: ACCESS-OM2-025, project: CMIP6, exp: omip2, activity: OMIP,
+          ensemble: r1i1p1f1, grid: gn, timerange: 0306/0306,
+          supplementary_variables: [{short_name: thkcello, mip: Ofx}]}
+      - {dataset: ACCESS-OM2, project: CMIP6, exp: omip2, activity: OMIP,
+              ensemble: r1i1p1f1, grid: gn, start_year: 0356, end_year: 0366,
+              supplementary_variables: [{short_name: thkcello, mip: Ofx}]}
+    scripts:
+      plot_script:
+        script: ocean_south/plot_kinetic_energy.py
\ No newline at end of file