diff --git a/exercises/00_final/precipitation_climatology.py b/exercises/00_final/precipitation_climatology.py
index 55db37d..dca7ced 100644
--- a/exercises/00_final/precipitation_climatology.py
+++ b/exercises/00_final/precipitation_climatology.py
@@ -12,77 +12,82 @@
 import regionmask
 
 
-def convert_pr_units(darray):
+def convert_precipitation_units(precipitation_in_kg_per_m_squared_s):
     """
     Convert precipitation units from [kg m-2 s-1] to [mm day-1].
 
     Parameters
     ----------
-    darray : xarray.DataArray
-        xarray DataArray containing model precipitation data
+    precipitation_in_kg_per_m_squared_s : xarray.DataArray
+        xarray DataArray containing model precipitation data in kg m-2 s-1
 
     Returns
     -------
-    darray : xarray.DataArray
-        the input DataArray with precipitation units modified
+    precipitation_in_mm_per_day : xarray.DataArray
+        the input DataArray with precipitation units modified to mm day-1
     """
     # density 1000 kg m-3 => 1 kg m-2 == 1 mm
     # There are 60*60*24 = 86400 seconds per day
-    darray.data = darray.data * 86400
-    darray.attrs["units"] = "mm/day"
+    precipitation_in_mm_per_day = precipitation_in_kg_per_m_squared_s * 86400
 
-    assert (
-        darray.data.min() >= 0.0
-    ), "There is at least one negative precipitation value"
-    assert darray.data.max() < 2000, "There is a precipitation value/s > 2000 mm/day"
+    precipitation_in_mm_per_day.attrs["units"] = "mm/day"
 
-    return darray
+    if precipitation_in_mm_per_day.data.min() < 0.0:
+        raise ValueError("There is at least one negative precipitation value")
+    if precipitation_in_mm_per_day.data.max() > 2000:
+        raise ValueError("There is a precipitation value/s > 2000 mm/day")
 
+    return precipitation_in_mm_per_day
 
-def plot_zonal(data):
+
+def plot_zonally_averaged_precipitation(precipitation_data):
     """
     Plot zonally-averaged precipitation data and save to file.
 
     Parameters
     ----------
-    data : xarray.DataArray
-        xarray DataArray containing model data
+    precipitation_data : xarray.DataSet
+        xarray DataSet containing precipitation model data, specifying precipitation in
+        [kg m-2 s-1] at given latitudes, longitudes and time. The Dataset should contain
+        four aligned DataArrays: precipitation, latitude, longitude and time.
 
     Returns
     -------
     None
 
     """
-    zonal_pr = data["pr"].mean("lon", keep_attrs=True)
+    zonal_precipitation = precipitation_data["pr"].mean("lon", keep_attrs=True)
 
-    fig, ax = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
+    figure, axes = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
 
-    zonal_pr.sel(lat=[0]).plot.line(ax=ax[0], hue="lat")
-    zonal_pr.sel(lat=[-20, 20]).plot.line(ax=ax[1], hue="lat")
-    zonal_pr.sel(lat=[-45, 45]).plot.line(ax=ax[2], hue="lat")
-    zonal_pr.sel(lat=[-70, 70]).plot.line(ax=ax[3], hue="lat")
+    zonal_precipitation.sel(lat=[0]).plot.line(ax=axes[0], hue="lat")
+    zonal_precipitation.sel(lat=[-20, 20]).plot.line(ax=axes[1], hue="lat")
+    zonal_precipitation.sel(lat=[-45, 45]).plot.line(ax=axes[2], hue="lat")
+    zonal_precipitation.sel(lat=[-70, 70]).plot.line(ax=axes[3], hue="lat")
 
     plt.tight_layout()
-    for axis in ax:
+    for axis in axes:
         axis.set_ylim(0.0, 1.0e-4)
         axis.grid()
     plt.savefig("zonal.png", dpi=200)  # Save figure to file
 
-    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
+    figure, axes = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
 
-    zonal_pr.T.plot()
+    zonal_precipitation.T.plot()
 
     plt.savefig("zonal_map.png", dpi=200)  # Save figure to file
 
 
-def get_country_ann_avg(data, countries):
+def get_country_annual_average(precipitation_data, countries):
     """
     Calculate annual precipitation averages for countries and save to file.
 
     Parameters
     ----------
-    data : xarray.DataArray
-        xarray DataArray containing model data
+    precipitation_data : xarray.DataSet
+        xarray DataSet containing precipitation model data, specifying precipitation in
+        [kg m-2 s-1] at given latitudes, longitudes and time. The Dataset should contain
+        four aligned DataArrays: precipitation, latitude, longitude and time.
     countries : dict(str: str)
         dictionary mapping country names to regionmask codes. For a list see:
         regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions
@@ -92,51 +97,45 @@ def get_country_ann_avg(data, countries):
     None
 
     """
-    data_avg = data["pr"].groupby("time.year").mean("time", keep_attrs=True)
-    data_avg = convert_pr_units(data_avg)
-
-    land = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(data_avg)
-
-    with open("data.txt", "w", encoding="utf-8") as datafile:
-        for k, v in countries.items():
-            data_avg_mask = data_avg.where(land.cf == v)
-
-            # Debugging - plot countries to make sure mask works correctly
-            # fig, geo_axes = plt.subplots(nrows=1, ncols=1, figsize=(12,5),
-            #                 subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
-            # data_avg_mask.sel(year = 2010).plot.contourf(ax=geo_axes,
-            #                                       extend='max',
-            #                                       transform=ccrs.PlateCarree(),
-            #                                       cbar_kwargs={'label': data_avg_mask.units},
-            #                                       cmap=cmocean.cm.haline_r)
-            # geo_axes.add_feature(cfeature.COASTLINE, lw=0.5)
-            # gl = geo_axes.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
-            # linewidth=2, color='gray', alpha=0.5, linestyle='--')
-            # gl.top_labels = False
-            # gl.left_labels = True
-            # gl.xlocator = mticker.FixedLocator([-180, -90, 0, 90])
-            # gl.ylocator = mticker.FixedLocator([-66, -23, 0, 23, 66])
-            # gl.xformatter = LONGITUDE_FORMATTER
-            # gl.yformatter = LATITUDE_FORMATTER
-            # gl.xlabel_style = {'size': 15, 'color': 'gray'}
-            # gl.ylabel_style = {'size': 15, 'color': 'gray'}
-            # print("show %s" %k)
-            # plt.show()
-
-            for yr in data_avg_mask.year.values:
-                precip = data_avg_mask.sel(year=yr).mean().values
-                datafile.write(f"{k.ljust(25)} {yr} : {precip:2.3f} mm/day\n")
+    annual_average_precipitation = (
+        precipitation_data["pr"].groupby("time.year").mean("time", keep_attrs=True)
+    )
+    annual_average_precipitation = convert_precipitation_units(
+        annual_average_precipitation
+    )
+
+    country_mask = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(
+        annual_average_precipitation
+    )
+
+    with open(
+        "annual_average_precipitation_by_country.txt", "w", encoding="utf-8"
+    ) as datafile:
+        for country_name, country_code in countries.items():
+            country_annual_average_precipitation = annual_average_precipitation.where(
+                country_mask.cf == country_code
+            )
+
+            for year in country_annual_average_precipitation.year.values:
+                precipitation = (
+                    country_annual_average_precipitation.sel(year=year).mean().values
+                )
+                datafile.write(
+                    f"{country_name.ljust(25)} {year} : {precipitation:2.3f} mm/day\n"
+                )
             datafile.write("\n")
 
 
-def plot_enso(data):
+def plot_enso_hovmoller_diagram(precipitation_data):
     """
     Plot Hovmöller diagram of equatorial precipitation to visualise ENSO.
 
     Parameters
     ----------
-    data : xarray.DataArray
-        xarray DataArray containing model data
+    precipitation_data : xarray.DataSet
+       xarray DataSet containing precipitation model data, specifying precipitation in
+        [kg m-2 s-1] at given latitudes, longitudes and time. The Dataset should contain
+        four aligned DataArrays: precipitation, latitude, longitude and time.
 
     Returns
     -------
@@ -144,7 +143,7 @@ def plot_enso(data):
 
     """
     enso = (
-        data["pr"]
+        precipitation_data["pr"]
         .sel(lat=slice(-1, 1))
         .sel(lon=slice(120, 280))
         .mean(dim="lat", keep_attrs=True)
@@ -154,21 +153,29 @@ def plot_enso(data):
     plt.savefig("enso.png", dpi=200)  # Save figure to file
 
 
-def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
+def create_precipitation_climatology_plot(
+    seasonal_average_precipitation,
+    model_name,
+    season,
+    mask=None,
+    plot_gridlines=False,
+    levels=None,
+):
     """
     Plot the precipitation climatology.
 
     Parameters
     ----------
-    clim : xarray.DataArray
-        Precipitation climatology data
-    model : str
+    seasonal_average_precipitation : xarray.DataArray
+        Precipitation climatology data. Seasonally averaged precipitation data.
+    model_name : str
         Name of the climate model
     season : str
         Climatological season (one of DJF, MAM, JJA, SON)
     mask : optional str
         mask to apply to plot (one of "land" or "ocean")
-    gridlines : bool
+    plot_gridlines : bool
+
         Select whether to plot gridlines
     levels : list
         Tick mark values for the colorbar
@@ -188,12 +195,12 @@ def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
         subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)},
     )
 
-    clim.sel(season=season).plot.contourf(
+    seasonal_average_precipitation.sel(season=season).plot.contourf(
         ax=geo_axes,
         levels=levels,
         extend="max",
         transform=ccrs.PlateCarree(),
-        cbar_kwargs={"label": clim.units},
+        cbar_kwargs={"label": seasonal_average_precipitation.units},
         cmap=cmocean.cm.rain,
     )
 
@@ -211,8 +218,8 @@ def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
             alpha=0.75,
         )
 
-    if gridlines:
-        gl = geo_axes.gridlines(
+    if plot_gridlines:
+        gridlines = geo_axes.gridlines(
             crs=ccrs.PlateCarree(),
             draw_labels=True,
             linewidth=2,
@@ -220,27 +227,27 @@ def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
             alpha=0.5,
             linestyle="--",
         )
-        gl.top_labels = False
-        gl.left_labels = True
+        gridlines.top_labels = False
+        gridlines.left_labels = True
         # gl.xlines = False
-        gl.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
-        gl.ylocator = mticker.FixedLocator(
+        gridlines.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
+        gridlines.ylocator = mticker.FixedLocator(
             [-66, -23, 0, 23, 66]
         )  # Tropics & Polar Circles
-        gl.xformatter = LONGITUDE_FORMATTER
-        gl.yformatter = LATITUDE_FORMATTER
-        gl.xlabel_style = {"size": 15, "color": "gray"}
-        gl.ylabel_style = {"size": 15, "color": "gray"}
+        gridlines.xformatter = LONGITUDE_FORMATTER
+        gridlines.yformatter = LATITUDE_FORMATTER
+        gridlines.xlabel_style = {"size": 15, "color": "gray"}
+        gridlines.ylabel_style = {"size": 15, "color": "gray"}
 
-    title = f"{model} precipitation climatology ({season})"
+    title = f"{model_name} precipitation climatology ({season})"
     plt.title(title)
 
 
 def main(
-    pr_file,
+    precipitation_netcdf_file,
     season="DJF",
     output_file="output.png",
-    gridlines=False,
+    plot_gridlines=False,
     mask=None,
     cbar_levels=None,
     countries=None,
@@ -250,13 +257,13 @@ def main(
 
     Parameters
     ----------
-    pr_file : str
+    precipitation_netcdf_file : str
         netCDF filename to read precipitation data from
     season : optional str
         Climatological season (one of DJF, MAM, JJA, SON)
     output_file : optional str
         filename to save main image to
-    gridlines : optional bool
+    plot_gridlines : optional bool
         Select whether to plot gridlines
     mask : optional str
         mask to apply to plot (one of "land" or "ocean")
@@ -274,34 +281,38 @@ def main(
     if countries is None:
         countries = {"United Kingdom": "GB"}
 
-    dset = xr.open_dataset(pr_file)
+    precipitation_data = xr.open_dataset(precipitation_netcdf_file)
 
-    plot_zonal(dset)
-    plot_enso(dset)
-    get_country_ann_avg(dset, countries)
+    plot_zonally_averaged_precipitation(precipitation_data)
+    plot_enso_hovmoller_diagram(precipitation_data)
+    get_country_annual_average(precipitation_data, countries)
 
-    clim = dset["pr"].groupby("time.season").mean("time", keep_attrs=True)
+    seasonal_average_precipitation = (
+        precipitation_data["pr"].groupby("time.season").mean("time", keep_attrs=True)
+    )
 
     try:
-        input_units = clim.attrs["units"]
+        input_units = seasonal_average_precipitation.attrs["units"]
     except KeyError as exc:
         raise KeyError(
             "Precipitation variable in {pr_file} must have a units attribute"
         ) from exc
 
     if input_units == "kg m-2 s-1":
-        clim = convert_pr_units(clim)
+        seasonal_average_precipitation = convert_precipitation_units(
+            seasonal_average_precipitation
+        )
     elif input_units == "mm/day":
         pass
     else:
         raise ValueError("""Input units are not 'kg m-2 s-1' or 'mm/day'""")
 
-    create_plot(
-        clim,
-        dset.attrs["source_id"],
+    create_precipitation_climatology_plot(
+        seasonal_average_precipitation,
+        precipitation_data.attrs["source_id"],
         season,
         mask=mask,
-        gridlines=gridlines,
+        plot_gridlines=plot_gridlines,
         levels=cbar_levels,
     )
 
@@ -316,8 +327,8 @@ def main(
         config_name = "default_config"
     else:
         print(f"Using configuration in '{config_name}.json'.")
-    configfile = config_name + ".json"
-    with open(configfile, encoding="utf-8") as json_file:
+    config_file = config_name + ".json"
+    with open(config_file, encoding="utf-8") as json_file:
         config = json.load(json_file)
 
     output_filename = f"{config_name}_output.png"
@@ -329,6 +340,6 @@ def main(
         season=config["season_to_plot"],
         output_file=output_filename,
         mask=config["mask_id"],
-        gridlines=config["gridlines_on"],
+        plot_gridlines=config["gridlines_on"],
         countries=config["countries_to_record"],
     )
diff --git a/exercises/01_base_code/precipitation_climatology.py b/exercises/01_base_code/precipitation_climatology.py
index cc4d435..4dd1257 100644
--- a/exercises/01_base_code/precipitation_climatology.py
+++ b/exercises/01_base_code/precipitation_climatology.py
@@ -20,8 +20,10 @@ def convert_pr_units(darray):
     darray.data = darray.data * 86400
     darray.attrs['units'] = 'mm/day'
     
-    assert darray.data.min() >= 0.0, 'There is at least one negative precipitation value'
-    assert darray.data.max() < 2000, 'There is a precipitation value/s > 2000 mm/day'
+    if darray.data.min() < 0.0:
+        raise ValueError('There is at least one negative precipitation value')
+    if darray.data.max() > 2000:
+        raise ValueError('There is a precipitation value/s > 2000 mm/day')
     
     return darray
 
@@ -65,7 +67,6 @@ def plot_zonal(data):
 
 
 def get_country_ann_avg(data, countries):
-
     data_avg = data['pr'].groupby('time.year').mean('time', keep_attrs=True)
     data_avg = convert_pr_units(data_avg)
 
@@ -113,7 +114,6 @@ def get_country_ann_avg(data, countries):
 
 
 def plot_enso(data):
-
     enso = data['pr'].sel(lat=slice(-1, 1)).sel(lon=slice(120, 280)).mean(dim="lat", keep_attrs=True)
     # print(enso)
     # .groupby('time.year').mean('time', keep_attrs=True)
@@ -238,7 +238,6 @@ def main(pr_file, season="DJF", output_file="output.png", gridlines=False, mask=
     plt.savefig(output_file, dpi=200)  # Save figure to file
 
 if __name__ == '__main__':
-
     input_file = "../../data/pr_Amon_ACCESS-ESM1-5_historical_r1i1p1f1_gn_201001-201412.nc"
     # season_to_plot = "DJF"
     # season_to_plot = "MAM"
diff --git a/exercises/02_formatting/precipitation_climatology.py b/exercises/02_formatting/precipitation_climatology.py
index cc4d435..4dd1257 100644
--- a/exercises/02_formatting/precipitation_climatology.py
+++ b/exercises/02_formatting/precipitation_climatology.py
@@ -20,8 +20,10 @@ def convert_pr_units(darray):
     darray.data = darray.data * 86400
     darray.attrs['units'] = 'mm/day'
     
-    assert darray.data.min() >= 0.0, 'There is at least one negative precipitation value'
-    assert darray.data.max() < 2000, 'There is a precipitation value/s > 2000 mm/day'
+    if darray.data.min() < 0.0:
+        raise ValueError('There is at least one negative precipitation value')
+    if darray.data.max() > 2000:
+        raise ValueError('There is a precipitation value/s > 2000 mm/day')
     
     return darray
 
@@ -65,7 +67,6 @@ def plot_zonal(data):
 
 
 def get_country_ann_avg(data, countries):
-
     data_avg = data['pr'].groupby('time.year').mean('time', keep_attrs=True)
     data_avg = convert_pr_units(data_avg)
 
@@ -113,7 +114,6 @@ def get_country_ann_avg(data, countries):
 
 
 def plot_enso(data):
-
     enso = data['pr'].sel(lat=slice(-1, 1)).sel(lon=slice(120, 280)).mean(dim="lat", keep_attrs=True)
     # print(enso)
     # .groupby('time.year').mean('time', keep_attrs=True)
@@ -238,7 +238,6 @@ def main(pr_file, season="DJF", output_file="output.png", gridlines=False, mask=
     plt.savefig(output_file, dpi=200)  # Save figure to file
 
 if __name__ == '__main__':
-
     input_file = "../../data/pr_Amon_ACCESS-ESM1-5_historical_r1i1p1f1_gn_201001-201412.nc"
     # season_to_plot = "DJF"
     # season_to_plot = "MAM"
diff --git a/exercises/03_naming/README.md b/exercises/03_naming/README.md
new file mode 100644
index 0000000..abb99ba
--- /dev/null
+++ b/exercises/03_naming/README.md
@@ -0,0 +1,20 @@
+# Exercise 3 - Naming for code clarity
+
+Look through the code for any names of methods or variables that could be improved or
+clarified and update them. Note if you are using an IDE like Intellij or VSCode, 
+you can use automatic renaming.
+Can you find an example from each of the suggestions listed below?
+Does this make the code easier to follow?
+
+Consider the following:
+
+- The name should show the intention, think about how someone else might read it (this could be future you)
+- Use pronounceable names e.g. `mass` not `ms`, `stem` not `stm`
+- avoid abbreviations and single letter variable names where possible
+- One word per concept e.g. choose one of `put`, `insert`, `add` in the same code base
+- Use names that can be searched
+- Describe content rather than storage type
+- Naming booleans, use prefixes like `is`, `has` or `can` and avoid negations like `not_green`
+- Plurals to indicate groups, e.g. a list of dog objects would be `dogs`, not `dog_list`
+- Keep it simple and use technical terms where appropriate
+- Use explaining variables
\ No newline at end of file
diff --git a/exercises/03_linting/precipitation_climatology.py b/exercises/03_naming/precipitation_climatology.py
similarity index 97%
rename from exercises/03_linting/precipitation_climatology.py
rename to exercises/03_naming/precipitation_climatology.py
index fde960c..083b5d0 100644
--- a/exercises/03_linting/precipitation_climatology.py
+++ b/exercises/03_naming/precipitation_climatology.py
@@ -19,10 +19,10 @@ def convert_pr_units(darray):
     darray.data = darray.data * 86400
     darray.attrs["units"] = "mm/day"
 
-    assert (
-        darray.data.min() >= 0.0
-    ), "There is at least one negative precipitation value"
-    assert darray.data.max() < 2000, "There is a precipitation value/s > 2000 mm/day"
+    if darray.data.min() < 0.0:
+        raise ValueError("There is at least one negative precipitation value")
+    if darray.data.max() > 2000:
+        raise ValueError("There is a precipitation value/s > 2000 mm/day")
 
     return darray
 
diff --git a/exercises/03_linting/README.md b/exercises/04_linting/README.md
similarity index 98%
rename from exercises/03_linting/README.md
rename to exercises/04_linting/README.md
index cccd7dc..486d41a 100644
--- a/exercises/03_linting/README.md
+++ b/exercises/04_linting/README.md
@@ -1,4 +1,4 @@
-# Exercise 3 - Linting
+# Exercise 4 - Linting
 
 This exercise contains the original code after applying black.
 
diff --git a/exercises/03_linting/long_func.py b/exercises/04_linting/long_func.py
similarity index 100%
rename from exercises/03_linting/long_func.py
rename to exercises/04_linting/long_func.py
diff --git a/exercises/04_linting/precipitation_climatology.py b/exercises/04_linting/precipitation_climatology.py
new file mode 100644
index 0000000..9bc9fb0
--- /dev/null
+++ b/exercises/04_linting/precipitation_climatology.py
@@ -0,0 +1,318 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import xarray as xr
+import scipy
+import cf_xarray
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import cmocean
+import regionmask
+
+
+import matplotlib.ticker as mticker
+from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
+
+
+def convert_precipitation_units(precipitation_in_kg_per_m_squared_s):
+    """Convert kg m-2 s-1 to mm day-1."""
+
+    precipitation_in_mm_per_day = precipitation_in_kg_per_m_squared_s * 86400
+
+    precipitation_in_mm_per_day.attrs["units"] = "mm/day"
+
+    if precipitation_in_mm_per_day.data.min() < 0.0:
+        raise ValueError("There is at least one negative precipitation value")
+    if precipitation_in_mm_per_day.data.max() > 2000:
+        raise ValueError("There is a precipitation value/s > 2000 mm/day")
+
+    return precipitation_in_mm_per_day
+
+
+def apply_mask(darray, sftlf_file, realm):
+    # Function to mask ocean or land using a sftlf (land surface fraction) file.
+    # Inputs:
+    #  darray: Data to mask
+    #  sftlf_file: Land surface fraction file
+    #  realm: Realm to mask
+
+    # This is now done using cartopy package with a single line.
+    pass
+
+
+def plot_zonally_averaged_precipitation(precipitation_data):
+    # print(data)
+    zonal_precipitation = precipitation_data["pr"].mean("lon", keep_attrs=True)
+
+    figure, axes = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
+
+    zonal_precipitation.sel(lat=[0]).plot.line(ax=axes[0], hue="lat")
+    zonal_precipitation.sel(lat=[-20, 20]).plot.line(ax=axes[1], hue="lat")
+    zonal_precipitation.sel(lat=[-45, 45]).plot.line(ax=axes[2], hue="lat")
+    zonal_precipitation.sel(lat=[-70, 70]).plot.line(ax=axes[3], hue="lat")
+
+    plt.tight_layout()
+    for axis in axes:
+        axis.set_ylim(0.0, 1.0e-4)
+        axis.grid()
+    plt.savefig("zonal.png", dpi=200)  # Save figure to file
+
+    figure, axes = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
+
+    zonal_precipitation.T.plot()
+
+    plt.savefig("zonal_map.png", dpi=200)  # Save figure to file
+
+
+def get_country_annual_average(data, countries):
+    annual_average_precipitation = (
+        data["pr"].groupby("time.year").mean("time", keep_attrs=True)
+    )
+    annual_average_precipitation = convert_precipitation_units(
+        annual_average_precipitation
+    )
+
+    country_mask = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(
+        annual_average_precipitation
+    )
+
+    # List possible locations to plot
+    # [print(k, v) for k, v in regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions.items()]
+
+    with open("annual_average_precipitation_by_country.txt", "w") as datafile:
+        for country_name, country_code in countries.items():
+            # land.plot(ax=geo_axes, add_label=False, fc="white", lw=2, alpha=0.5)
+            # clim = clim.where(ocean == "South Pacific Ocean")
+            country_annual_average_precipitation = annual_average_precipitation.where(
+                country_mask.cf == country_code
+            )
+
+            # Debugging - plot countries to make sure mask works correctly
+            # fig, geo_axes = plt.subplots(nrows=1, ncols=1, figsize=(12,5),
+            #                 subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
+            # data_avg_mask.sel(year = 2010).plot.contourf(ax=geo_axes,
+            #                                       extend='max',
+            #                                       transform=ccrs.PlateCarree(),
+            #                                       cbar_kwargs={'label': data_avg_mask.units},
+            #                                       cmap=cmocean.cm.haline_r)
+            # geo_axes.add_feature(cfeature.COASTLINE, lw=0.5)
+            # gl = geo_axes.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
+            # linewidth=2, color='gray', alpha=0.5, linestyle='--')
+            # gl.top_labels = False
+            # gl.left_labels = True
+            # gl.xlocator = mticker.FixedLocator([-180, -90, 0, 90])
+            # gl.ylocator = mticker.FixedLocator([-66, -23, 0, 23, 66])
+            # gl.xformatter = LONGITUDE_FORMATTER
+            # gl.yformatter = LATITUDE_FORMATTER
+            # gl.xlabel_style = {'size': 15, 'color': 'gray'}
+            # gl.ylabel_style = {'size': 15, 'color': 'gray'}
+            # print("show %s" %k)
+            # plt.show()
+
+            for year in country_annual_average_precipitation.year.values:
+                precipitation = (
+                    country_annual_average_precipitation.sel(year=year).mean().values
+                )
+                datafile.write(
+                    "{} {} : {:2.3f} mm/day\n".format(
+                        country_name.ljust(25), year, precipitation
+                    )
+                )
+            datafile.write("\n")
+
+
+def plot_enso_hovmoller_diagram(precipitation_data):
+    enso = (
+        precipitation_data["pr"]
+        .sel(lat=slice(-1, 1))
+        .sel(lon=slice(120, 280))
+        .mean(dim="lat", keep_attrs=True)
+    )
+    # print(enso)
+    # .groupby('time.year').mean('time', keep_attrs=True)
+
+    # # convert to dataframe:
+    # df = monthly_speed.reset_coords(drop=True).to_dataframe()
+    # # add year and month indices:
+    # df['month']=df.index.month
+    # df['year']=df.index.year
+    # # groupby month and year then mean:
+    # enso = enso.groupby(['time.year','time.month']).mean().unstack().T.droplevel(0)
+    # plot:
+    enso.plot()
+
+    plt.savefig("enso.png", dpi=200)  # Save figure to file
+
+
+def create_precipitation_climatology_plot(
+    seasonal_average_precipitation,
+    model_name,
+    season,
+    mask=None,
+    plot_gridlines=False,
+    levels=None,
+):
+    """Plot the precipitation climatology.
+
+    seasonal_average_precipitation : xarray.DataArray
+        Precipitation climatology data. Seasonally averaged precipitation data.
+    model_name : str
+    season (str): Season
+
+    plot_gridlines (bool): Select whether to plot gridlines
+    levels (list): Tick marks on the colorbar
+
+    """
+
+    # fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12,5), subplot_kw={'projection': "3d"})
+    # clim.sel(season=season).T.plot.surface()
+    # plt.show()
+
+    if not levels:
+        levels = np.arange(0, 13.5, 1.5)
+
+    fig, geo_axes = plt.subplots(
+        nrows=1,
+        ncols=1,
+        figsize=(12, 5),
+        subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)},
+    )
+
+    seasonal_average_precipitation.sel(season=season).plot.contourf(
+        ax=geo_axes,
+        levels=levels,
+        extend="max",
+        transform=ccrs.PlateCarree(),
+        cbar_kwargs={"label": seasonal_average_precipitation.units},
+        cmap=cmocean.cm.rain,
+    )
+
+    geo_axes.add_feature(
+        cfeature.COASTLINE, lw=2
+    )  # Add coastines using cartopy feature
+
+    if mask:
+        # Old approach of adding mask before combining into the below command.
+        # if mask == "ocean":
+        # old mask_feat = cfeature.NaturalEarthFeature("physical", "ocean", "110m")
+        # oldold geo_axes.add_feature(cfeature.NaturalEarthFeature("physical", "ocean", "110m"),
+        #                      ec="red", fc="yellow", lw=2, alpha=1.0)
+        # elif mask == "land":
+        # old mask_feat = cfeature.NaturalEarthFeature("physical", "land", "110m")
+        # oldold # geo_axes.add_feature(cfeature.NaturalEarthFeature("physical", "ocean", "110m"),
+        #                           ec="red", fc="yellow", lw=2, alpha=1.0)
+
+        # oldold else:
+        # oldold pass
+        # oldold raise ValueError("Unknown ")
+
+        # Mask out (fade) using 110m resolution data from cartopy.
+        geo_axes.add_feature(
+            cfeature.NaturalEarthFeature("physical", mask, "110m"),
+            ec=None,
+            fc="white",
+            lw=2,
+            alpha=0.75,
+        )
+
+    if plot_gridlines:
+        # If we want gridlines run the code to do this:
+        gridlines = geo_axes.gridlines(
+            crs=ccrs.PlateCarree(),
+            draw_labels=True,
+            linewidth=2,
+            color="gray",
+            alpha=0.5,
+            linestyle="--",
+        )
+        gridlines.top_labels = False
+        gridlines.left_labels = True
+        # gl.xlines = False
+        gridlines.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
+        gridlines.ylocator = mticker.FixedLocator(
+            [-66, -23, 0, 23, 66]
+        )  # Tropics & Polar Circles
+        gridlines.xformatter = LONGITUDE_FORMATTER
+        gridlines.yformatter = LATITUDE_FORMATTER
+        gridlines.xlabel_style = {"size": 15, "color": "gray"}
+        gridlines.ylabel_style = {"size": 15, "color": "gray"}
+
+    title = "{} precipitation climatology ({})".format(model_name, season)
+    plt.title(title)
+    # print("\n\n{}\n\n".format(clim.mean()))
+
+
+def main(
+    precipitation_netcdf_file,
+    season="DJF",
+    output_file="output.png",
+    plot_gridlines=False,
+    mask=None,
+    cbar_levels=None,
+    countries={"United Kingdom": "GB"},
+):
+    """Run the program."""
+
+    precipitation_data = xr.open_dataset(precipitation_netcdf_file)
+
+    plot_zonally_averaged_precipitation(precipitation_data)
+    plot_enso_hovmoller_diagram(precipitation_data)
+    get_country_annual_average(precipitation_data, countries)
+
+    seasonal_average_precipitation = (
+        precipitation_data["pr"].groupby("time.season").mean("time", keep_attrs=True)
+    )
+
+    try:
+        input_units = seasonal_average_precipitation.attrs["units"]
+    except KeyError:
+        raise KeyError(
+            "Precipitation variable in {pr_file} must have a units attribute"
+        )
+
+    if input_units == "kg m-2 s-1":
+        seasonal_average_precipitation = convert_precipitation_units(
+            seasonal_average_precipitation
+        )
+    elif input_units == "mm/day":
+        pass
+    else:
+        raise ValueError("""Input units are not 'kg m-2 s-1' or 'mm/day'""")
+
+    create_precipitation_climatology_plot(
+        seasonal_average_precipitation,
+        precipitation_data.attrs["source_id"],
+        season,
+        mask=mask,
+        plot_gridlines=plot_gridlines,
+        levels=cbar_levels,
+    )
+
+    plt.savefig(output_file, dpi=200)  # Save figure to file
+
+
+if __name__ == "__main__":
+    input_file = (
+        "../../data/pr_Amon_ACCESS-ESM1-5_historical_r1i1p1f1_gn_201001-201412.nc"
+    )
+    # season_to_plot = "DJF"
+    # season_to_plot = "MAM"
+    season_to_plot = "JJA"
+    # season_to_plot = "SON"
+    output_filename = "output.png"
+    gridlines_on = True
+    mask_id = "ocean"
+    cbar_levels = None
+    countries = {
+        "United Kingdom": "GB",
+        "United States of America": "US",
+        "Antarctica": "AQ",
+        "South Africa": "ZA",
+    }
+
+    main(
+        input_file,
+        season=season_to_plot,
+        mask=mask_id,
+        plot_gridlines=gridlines_on,
+        countries=countries,
+    )
diff --git a/exercises/05_better_code/precipitation_climatology.py b/exercises/05_better_code/precipitation_climatology.py
deleted file mode 100644
index 6c827bd..0000000
--- a/exercises/05_better_code/precipitation_climatology.py
+++ /dev/null
@@ -1,334 +0,0 @@
-"""Routines for analysing precipitation climatology from ESM runs."""
-
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.ticker as mticker
-import xarray as xr
-import cartopy.crs as ccrs
-import cartopy.feature as cfeature
-from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
-import cmocean
-import regionmask
-
-
-def convert_pr_units(darray):
-    """
-    Convert precipitation units from [kg m-2 s-1] to [mm day-1].
-
-    Parameters
-    ----------
-    darray : xarray.DataArray
-        xarray DataArray containing model precipitation data
-
-    Returns
-    -------
-    darray : xarray.DataArray
-        the input DataArray with precipitation units modified
-    """
-    # density 1000 kg m-3 => 1 kg m-2 == 1 mm
-    # There are 60*60*24 = 86400 seconds per day
-    darray.data = darray.data * 86400
-    darray.attrs["units"] = "mm/day"
-
-    assert (
-        darray.data.min() >= 0.0
-    ), "There is at least one negative precipitation value"
-    assert darray.data.max() < 2000, "There is a precipitation value/s > 2000 mm/day"
-
-    return darray
-
-
-def plot_zonal(data):
-    """
-    Plot zonally-averaged precipitation data and save to file.
-
-    Parameters
-    ----------
-    data : xarray.DataArray
-        xarray DataArray containing model data
-
-    Returns
-    -------
-    None
-
-    """
-    zonal_pr = data["pr"].mean("lon", keep_attrs=True)
-
-    fig, ax = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
-
-    zonal_pr.sel(lat=[0]).plot.line(ax=ax[0], hue="lat")
-    zonal_pr.sel(lat=[-20, 20]).plot.line(ax=ax[1], hue="lat")
-    zonal_pr.sel(lat=[-45, 45]).plot.line(ax=ax[2], hue="lat")
-    zonal_pr.sel(lat=[-70, 70]).plot.line(ax=ax[3], hue="lat")
-
-    plt.tight_layout()
-    for axis in ax:
-        axis.set_ylim(0.0, 1.0e-4)
-        axis.grid()
-    plt.savefig("zonal.png", dpi=200)  # Save figure to file
-
-    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
-
-    zonal_pr.T.plot()
-
-    plt.savefig("zonal_map.png", dpi=200)  # Save figure to file
-
-
-def get_country_ann_avg(data, countries):
-    """
-    Calculate annual precipitation averages for countries and save to file.
-
-    Parameters
-    ----------
-    data : xarray.DataArray
-        xarray DataArray containing model data
-    countries : dict(str: str)
-        dictionary mapping country names to regionmask codes. For a list see:
-        regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions
-
-    Returns
-    -------
-    None
-
-    """
-    data_avg = data["pr"].groupby("time.year").mean("time", keep_attrs=True)
-    data_avg = convert_pr_units(data_avg)
-
-    land = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(data_avg)
-
-    with open("data.txt", "w", encoding="utf-8") as datafile:
-        for k, v in countries.items():
-            data_avg_mask = data_avg.where(land.cf == v)
-
-            # Debugging - plot countries to make sure mask works correctly
-            # fig, geo_axes = plt.subplots(nrows=1, ncols=1, figsize=(12,5),
-            #                 subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
-            # data_avg_mask.sel(year = 2010).plot.contourf(ax=geo_axes,
-            #                                       extend='max',
-            #                                       transform=ccrs.PlateCarree(),
-            #                                       cbar_kwargs={'label': data_avg_mask.units},
-            #                                       cmap=cmocean.cm.haline_r)
-            # geo_axes.add_feature(cfeature.COASTLINE, lw=0.5)
-            # gl = geo_axes.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
-            # linewidth=2, color='gray', alpha=0.5, linestyle='--')
-            # gl.top_labels = False
-            # gl.left_labels = True
-            # gl.xlocator = mticker.FixedLocator([-180, -90, 0, 90])
-            # gl.ylocator = mticker.FixedLocator([-66, -23, 0, 23, 66])
-            # gl.xformatter = LONGITUDE_FORMATTER
-            # gl.yformatter = LATITUDE_FORMATTER
-            # gl.xlabel_style = {'size': 15, 'color': 'gray'}
-            # gl.ylabel_style = {'size': 15, 'color': 'gray'}
-            # print("show %s" %k)
-            # plt.show()
-
-            for yr in data_avg_mask.year.values:
-                precip = data_avg_mask.sel(year=yr).mean().values
-                datafile.write(
-                    "{} {} : {:2.3f} mm/day\n".format(k.ljust(25), yr, precip)
-                )
-            datafile.write("\n")
-
-
-def plot_enso(data):
-    """
-    Plot Hovmöller diagram of equatorial precipitation to visualise ENSO.
-
-    Parameters
-    ----------
-    data : xarray.DataArray
-        xarray DataArray containing model data
-
-    Returns
-    -------
-    None
-
-    """
-    enso = (
-        data["pr"]
-        .sel(lat=slice(-1, 1))
-        .sel(lon=slice(120, 280))
-        .mean(dim="lat", keep_attrs=True)
-    )
-
-    enso.plot()
-    plt.savefig("enso.png", dpi=200)  # Save figure to file
-
-
-def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
-    """
-    Plot the precipitation climatology.
-
-    Parameters
-    ----------
-    clim : xarray.DataArray
-        Precipitation climatology data
-    model : str
-        Name of the climate model
-    season : str
-        Climatological season (one of DJF, MAM, JJA, SON)
-    mask : optional str
-        mask to apply to plot (one of "land" or "ocean")
-    gridlines : bool
-        Select whether to plot gridlines
-    levels : list
-        Tick mark values for the colorbar
-
-    Returns
-    -------
-    None
-
-    """
-    if not levels:
-        levels = np.arange(0, 13.5, 1.5)
-
-    fig, geo_axes = plt.subplots(
-        nrows=1,
-        ncols=1,
-        figsize=(12, 5),
-        subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)},
-    )
-
-    clim.sel(season=season).plot.contourf(
-        ax=geo_axes,
-        levels=levels,
-        extend="max",
-        transform=ccrs.PlateCarree(),
-        cbar_kwargs={"label": clim.units},
-        cmap=cmocean.cm.rain,
-    )
-
-    geo_axes.add_feature(
-        cfeature.COASTLINE, lw=2
-    )  # Add coastines using cartopy feature
-
-    if mask:
-        # Mask out (fade) using 110m resolution data from cartopy.
-        geo_axes.add_feature(
-            cfeature.NaturalEarthFeature("physical", mask, "110m"),
-            ec=None,
-            fc="white",
-            lw=2,
-            alpha=0.75,
-        )
-
-    if gridlines:
-        gl = geo_axes.gridlines(
-            crs=ccrs.PlateCarree(),
-            draw_labels=True,
-            linewidth=2,
-            color="gray",
-            alpha=0.5,
-            linestyle="--",
-        )
-        gl.top_labels = False
-        gl.left_labels = True
-        # gl.xlines = False
-        gl.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
-        gl.ylocator = mticker.FixedLocator(
-            [-66, -23, 0, 23, 66]
-        )  # Tropics & Polar Circles
-        gl.xformatter = LONGITUDE_FORMATTER
-        gl.yformatter = LATITUDE_FORMATTER
-        gl.xlabel_style = {"size": 15, "color": "gray"}
-        gl.ylabel_style = {"size": 15, "color": "gray"}
-
-    title = "{} precipitation climatology ({})".format(model, season)
-    plt.title(title)
-
-
-def main(
-    pr_file,
-    season="DJF",
-    output_file="output.png",
-    gridlines=False,
-    mask=None,
-    cbar_levels=None,
-    countries=None,
-):
-    """
-    Run the program for producing precipitation plots.
-
-    Parameters
-    ----------
-    pr_file : str
-        netCDF filename to read precipitation data from
-    season : optional str
-        Climatological season (one of DJF, MAM, JJA, SON)
-    output_file : optional str
-        filename to save main image to
-    gridlines : optional bool
-        Select whether to plot gridlines
-    mask : optional str
-        mask to apply to plot (one of "land" or "ocean")
-    cbar_levels : optional list
-        Tick mark values for the colorbar
-    countries : optional dict(str: str)
-        dictionary mapping country names to regionmask codes. For a list see:
-        regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions
-
-    Returns
-    -------
-    None
-
-    """
-    if countries is None:
-        countries = {"United Kingdom": "GB"}
-
-    dset = xr.open_dataset(pr_file)
-
-    plot_zonal(dset)
-    plot_enso(dset)
-    get_country_ann_avg(dset, countries)
-
-    clim = dset["pr"].groupby("time.season").mean("time", keep_attrs=True)
-
-    try:
-        input_units = clim.attrs["units"]
-    except KeyError as exc:
-        raise KeyError(
-            "Precipitation variable in {pr_file} must have a units attribute"
-        ) from exc
-
-    if input_units == "kg m-2 s-1":
-        clim = convert_pr_units(clim)
-    elif input_units == "mm/day":
-        pass
-    else:
-        raise ValueError("""Input units are not 'kg m-2 s-1' or 'mm/day'""")
-
-    create_plot(
-        clim,
-        dset.attrs["source_id"],
-        season,
-        mask=mask,
-        gridlines=gridlines,
-        levels=cbar_levels,
-    )
-
-    plt.savefig(output_file, dpi=200)
-
-
-if __name__ == "__main__":
-    input_file = (
-        "../../data/pr_Amon_ACCESS-ESM1-5_historical_r1i1p1f1_gn_201001-201412.nc"
-    )
-    season_to_plot = "JJA"
-    output_filename = "output.png"
-    gridlines_on = True
-    mask_id = "ocean"
-    colorbar_levels = None
-    countries_to_record = {
-        "United Kingdom": "GB",
-        "United States of America": "US",
-        "Antarctica": "AQ",
-        "South Africa": "ZA",
-    }
-
-    main(
-        input_file,
-        season=season_to_plot,
-        mask=mask_id,
-        gridlines=gridlines_on,
-        countries=countries_to_record,
-    )
diff --git a/exercises/04_docstrings_and_comments/README.md b/exercises/05_docstrings_and_comments/README.md
similarity index 96%
rename from exercises/04_docstrings_and_comments/README.md
rename to exercises/05_docstrings_and_comments/README.md
index 5ac67b5..c898e4e 100644
--- a/exercises/04_docstrings_and_comments/README.md
+++ b/exercises/05_docstrings_and_comments/README.md
@@ -1,4 +1,4 @@
-# Exercise 4 - Docstrings and Comments
+# Exercise 5 - Docstrings and Comments
 
 This exercise contains the code after addressing some of the issues raised by pylint.
 
diff --git a/exercises/04_docstrings_and_comments/gyroradius.py b/exercises/05_docstrings_and_comments/gyroradius.py
similarity index 100%
rename from exercises/04_docstrings_and_comments/gyroradius.py
rename to exercises/05_docstrings_and_comments/gyroradius.py
diff --git a/exercises/04_docstrings_and_comments/precipitation_climatology.py b/exercises/05_docstrings_and_comments/precipitation_climatology.py
similarity index 62%
rename from exercises/04_docstrings_and_comments/precipitation_climatology.py
rename to exercises/05_docstrings_and_comments/precipitation_climatology.py
index 9d9186f..8150cf8 100644
--- a/exercises/04_docstrings_and_comments/precipitation_climatology.py
+++ b/exercises/05_docstrings_and_comments/precipitation_climatology.py
@@ -9,58 +9,67 @@
 import regionmask
 
 
-def convert_pr_units(darray):
+def convert_precipitation_units(precipitation_in_kg_per_m_squared_s):
     """Convert kg m-2 s-1 to mm day-1."""
 
-    darray.data = darray.data * 86400
-    darray.attrs["units"] = "mm/day"
+    precipitation_in_mm_per_day = precipitation_in_kg_per_m_squared_s * 86400
 
-    assert (
-        darray.data.min() >= 0.0
-    ), "There is at least one negative precipitation value"
-    assert darray.data.max() < 2000, "There is a precipitation value/s > 2000 mm/day"
+    precipitation_in_mm_per_day.attrs["units"] = "mm/day"
 
-    return darray
+    if precipitation_in_mm_per_day.data.min() < 0.0:
+        raise ValueError("There is at least one negative precipitation value")
+    if precipitation_in_mm_per_day.data.max() > 2000:
+        raise ValueError("There is a precipitation value/s > 2000 mm/day")
 
+    return precipitation_in_mm_per_day
 
-def plot_zonal(data):
+
+def plot_zonally_averaged_precipitation(precipitation_data):
     # print(data)
-    zonal_pr = data["pr"].mean("lon", keep_attrs=True)
+    zonal_precipitation = precipitation_data["pr"].mean("lon", keep_attrs=True)
 
-    fig, ax = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
+    figure, axes = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
 
-    zonal_pr.sel(lat=[0]).plot.line(ax=ax[0], hue="lat")
-    zonal_pr.sel(lat=[-20, 20]).plot.line(ax=ax[1], hue="lat")
-    zonal_pr.sel(lat=[-45, 45]).plot.line(ax=ax[2], hue="lat")
-    zonal_pr.sel(lat=[-70, 70]).plot.line(ax=ax[3], hue="lat")
+    zonal_precipitation.sel(lat=[0]).plot.line(ax=axes[0], hue="lat")
+    zonal_precipitation.sel(lat=[-20, 20]).plot.line(ax=axes[1], hue="lat")
+    zonal_precipitation.sel(lat=[-45, 45]).plot.line(ax=axes[2], hue="lat")
+    zonal_precipitation.sel(lat=[-70, 70]).plot.line(ax=axes[3], hue="lat")
 
     plt.tight_layout()
-    for axis in ax:
+    for axis in axes:
         axis.set_ylim(0.0, 1.0e-4)
         axis.grid()
     plt.savefig("zonal.png", dpi=200)  # Save figure to file
 
-    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
+    figure, axes = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
 
-    zonal_pr.T.plot()
+    zonal_precipitation.T.plot()
 
     plt.savefig("zonal_map.png", dpi=200)  # Save figure to file
 
 
-def get_country_ann_avg(data, countries):
-    data_avg = data["pr"].groupby("time.year").mean("time", keep_attrs=True)
-    data_avg = convert_pr_units(data_avg)
+def get_country_annual_average(precipitation_data, countries):
+    annual_average_precipitation = (
+        precipitation_data["pr"].groupby("time.year").mean("time", keep_attrs=True)
+    )
+    annual_average_precipitation = convert_precipitation_units(
+        annual_average_precipitation
+    )
 
-    land = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(data_avg)
+    country_mask = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(
+        annual_average_precipitation
+    )
 
     # List possible locations to plot
     # [print(k, v) for k, v in regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions.items()]
 
-    with open("data.txt", "w", encoding="utf-8") as datafile:
-        for k, v in countries.items():
+    with open("annual_average_precipitation_by_country.txt", "w", encoding="utf-8") as datafile:
+        for country_name, country_code in countries.items():
             # land.plot(ax=geo_axes, add_label=False, fc="white", lw=2, alpha=0.5)
             # clim = clim.where(ocean == "South Pacific Ocean")
-            data_avg_mask = data_avg.where(land.cf == v)
+            country_annual_average_precipitation = annual_average_precipitation.where(
+                country_mask.cf == country_code
+            )
 
             # Debugging - plot countries to make sure mask works correctly
             # fig, geo_axes = plt.subplots(nrows=1, ncols=1, figsize=(12,5),
@@ -84,17 +93,21 @@ def get_country_ann_avg(data, countries):
             # print("show %s" %k)
             # plt.show()
 
-            for yr in data_avg_mask.year.values:
-                precip = data_avg_mask.sel(year=yr).mean().values
+            for year in country_annual_average_precipitation.year.values:
+                precip = (
+                    country_annual_average_precipitation.sel(year=year).mean().values
+                )
                 datafile.write(
-                    "{} {} : {:2.3f} mm/day\n".format(k.ljust(25), yr, precip)
+                    "{} {} : {:2.3f} mm/day\n".format(
+                        country_name.ljust(25), year, precip
+                    )
                 )
             datafile.write("\n")
 
 
-def plot_enso(data):
+def plot_enso_hovmoller_diagram(precipitation_data):
     enso = (
-        data["pr"]
+        precipitation_data["pr"]
         .sel(lat=slice(-1, 1))
         .sel(lon=slice(120, 280))
         .mean(dim="lat", keep_attrs=True)
@@ -115,7 +128,14 @@ def plot_enso(data):
     plt.savefig("enso.png", dpi=200)  # Save figure to file
 
 
-def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
+def create_precipitation_climatology_plot(
+    seasonal_average_precipitation,
+    model_name,
+    season,
+    mask=None,
+    plot_gridlines=False,
+    levels=None,
+):
     """Plot the precipitation climatology.
 
     clim (xarray.DataArray): Precipitation climatology data
@@ -141,12 +161,12 @@ def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
         subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)},
     )
 
-    clim.sel(season=season).plot.contourf(
+    seasonal_average_precipitation.sel(season=season).plot.contourf(
         ax=geo_axes,
         levels=levels,
         extend="max",
         transform=ccrs.PlateCarree(),
-        cbar_kwargs={"label": clim.units},
+        cbar_kwargs={"label": seasonal_average_precipitation.units},
         cmap=cmocean.cm.rain,
     )
 
@@ -178,9 +198,9 @@ def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
             alpha=0.75,
         )
 
-    if gridlines:
+    if plot_gridlines:
         # If we want gridlines run the code to do this:
-        gl = geo_axes.gridlines(
+        gridlines = geo_axes.gridlines(
             crs=ccrs.PlateCarree(),
             draw_labels=True,
             linewidth=2,
@@ -188,28 +208,28 @@ def create_plot(clim, model, season, mask=None, gridlines=False, levels=None):
             alpha=0.5,
             linestyle="--",
         )
-        gl.top_labels = False
-        gl.left_labels = True
+        gridlines.top_labels = False
+        gridlines.left_labels = True
         # gl.xlines = False
-        gl.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
-        gl.ylocator = mticker.FixedLocator(
+        gridlines.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
+        gridlines.ylocator = mticker.FixedLocator(
             [-66, -23, 0, 23, 66]
         )  # Tropics & Polar Circles
-        gl.xformatter = LONGITUDE_FORMATTER
-        gl.yformatter = LATITUDE_FORMATTER
-        gl.xlabel_style = {"size": 15, "color": "gray"}
-        gl.ylabel_style = {"size": 15, "color": "gray"}
+        gridlines.xformatter = LONGITUDE_FORMATTER
+        gridlines.yformatter = LATITUDE_FORMATTER
+        gridlines.xlabel_style = {"size": 15, "color": "gray"}
+        gridlines.ylabel_style = {"size": 15, "color": "gray"}
 
-    title = "{} precipitation climatology ({})".format(model, season)
+    title = "{} precipitation climatology ({})".format(model_name, season)
     plt.title(title)
     # print("\n\n{}\n\n".format(clim.mean()))
 
 
 def main(
-    pr_file,
+    precipitation_netcdf_file,
     season="DJF",
     output_file="output.png",
-    gridlines=False,
+    plot_gridlines=False,
     mask=None,
     cbar_levels=None,
     countries=None,
@@ -219,34 +239,38 @@ def main(
     if countries is None:
         countries = {"United Kingdom": "GB"}
 
-    dset = xr.open_dataset(pr_file)
+    precipitation_data = xr.open_dataset(precipitation_netcdf_file)
 
-    plot_zonal(dset)
-    plot_enso(dset)
-    get_country_ann_avg(dset, countries)
+    plot_zonally_averaged_precipitation(precipitation_data)
+    plot_enso_hovmoller_diagram(precipitation_data)
+    get_country_annual_average(precipitation_data, countries)
 
-    clim = dset["pr"].groupby("time.season").mean("time", keep_attrs=True)
+    seasonal_average_precipitation = (
+        precipitation_data["pr"].groupby("time.season").mean("time", keep_attrs=True)
+    )
 
     try:
-        input_units = clim.attrs["units"]
+        input_units = seasonal_average_precipitation.attrs["units"]
     except KeyError as exc:
         raise KeyError(
             "Precipitation variable in {pr_file} must have a units attribute"
         ) from exc
 
     if input_units == "kg m-2 s-1":
-        clim = convert_pr_units(clim)
+        seasonal_average_precipitation = convert_precipitation_units(
+            seasonal_average_precipitation
+        )
     elif input_units == "mm/day":
         pass
     else:
         raise ValueError("""Input units are not 'kg m-2 s-1' or 'mm/day'""")
 
-    create_plot(
-        clim,
-        dset.attrs["source_id"],
+    create_precipitation_climatology_plot(
+        seasonal_average_precipitation,
+        precipitation_data.attrs["source_id"],
         season,
         mask=mask,
-        gridlines=gridlines,
+        plot_gridlines=plot_gridlines,
         levels=cbar_levels,
     )
 
@@ -276,6 +300,6 @@ def main(
         input_file,
         season=season_to_plot,
         mask=mask_id,
-        gridlines=gridlines_on,
+        plot_gridlines=gridlines_on,
         countries=countries_to_record,
     )
diff --git a/exercises/05_better_code/README.md b/exercises/06_better_code/README.md
similarity index 97%
rename from exercises/05_better_code/README.md
rename to exercises/06_better_code/README.md
index 8225a97..fd3223d 100644
--- a/exercises/05_better_code/README.md
+++ b/exercises/06_better_code/README.md
@@ -1,4 +1,4 @@
-# Exercise 5 - Better Code
+# Exercise 6 - Better Code
 
 This exercise contains the code now complete with docstrings.
 
diff --git a/exercises/05_better_code/pendulum.py b/exercises/06_better_code/pendulum.py
similarity index 79%
rename from exercises/05_better_code/pendulum.py
rename to exercises/06_better_code/pendulum.py
index d314900..6badb17 100644
--- a/exercises/05_better_code/pendulum.py
+++ b/exercises/06_better_code/pendulum.py
@@ -3,13 +3,13 @@
 import numpy as np
 
 
-def get_period(l):
+def get_period(length):
     """
     Calculate the period of a pendulum.
 
     Parameters
     ----------
-    l : float
+    length : float
         length of the pendulum [m]
 
     Returns
@@ -17,16 +17,16 @@ def get_period(l):
     float
         period [s] for a swing of the pendulum
     """
-    return 2.0 * np.pi * np.sqrt(l / 9.81)
+    return 2.0 * np.pi * np.sqrt(length / 9.81)
 
 
-def max_height(l, theta):
+def max_height(length, theta):
     """
     Calculate the maximum height reached by a pendulum.
 
     Parameters
     ----------
-    l : float
+    length : float
         length of the pendulum [m]
     theta : float
         maximum angle of displacment of the pendulum [radians]
@@ -36,16 +36,16 @@ def max_height(l, theta):
     float
         maximum vertical height [m] of the pendulum
     """
-    return l * np.cos(theta)
+    return length * np.cos(theta)
 
 
-def max_speed(l, theta):
+def max_speed(length, theta):
     """
     Calculate the maximum speed of a pendulum.
 
     Parameters
     ----------
-    l : float
+    length : float
         length of the pendulum [m]
     theta : float
         maximum angle of displacment of the pendulum [radians]
@@ -55,18 +55,18 @@ def max_speed(l, theta):
     float
         maximum speed [m/s] of the pendulum
     """
-    return np.sqrt(2.0 * 9.81 * max_height(l, theta))
+    return np.sqrt(2.0 * 9.81 * max_height(length, theta))
 
 
-def energy(m, l, theta):
+def energy(mass, length, theta):
     """
     Calculate the energy of a pendulum.
 
     Parameters
     ----------
-    m : float
+    mass : float
         mass of the pendulum bob [kg]
-    l : float
+    length : float
         length of the pendulum [m]
     theta : float
         maximum angle of displacment of the pendulum [radians]
@@ -76,7 +76,7 @@ def energy(m, l, theta):
     float
         energy [kg . m2 /s2] of the pendulum
     """
-    return m * 9.81 * max_height(l, theta)
+    return mass * 9.81 * max_height(length, theta)
 
 
 def check_small_angle(theta):
@@ -98,13 +98,13 @@ def check_small_angle(theta):
     return False
 
 
-def bpm(l):
+def bpm(length):
     """
     Calculate pendulum frequency in beats per minute.
 
     Parameters
     ----------
-    l : float
+    length : float
         length of the pendulum [m]
 
     Returns
@@ -112,4 +112,4 @@ def bpm(l):
     float
         pendulum frequency in beats per minute [1 / min]
     """
-    return 60.0 / get_period(l)
+    return 60.0 / get_period(length)
diff --git a/exercises/06_better_code/precipitation_climatology.py b/exercises/06_better_code/precipitation_climatology.py
new file mode 100644
index 0000000..0fb0db9
--- /dev/null
+++ b/exercises/06_better_code/precipitation_climatology.py
@@ -0,0 +1,342 @@
+"""Routines for analysing precipitation climatology from ESM runs."""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mticker
+import xarray as xr
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
+import cmocean
+import regionmask
+
+
+def convert_precipitation_units(precipitation_in_kg_per_m_squared_s):
+    """
+    Convert precipitation units from [kg m-2 s-1] to [mm day-1].
+
+    Parameters
+    ----------
+    precipitation_in_kg_per_m_squared_s : xarray.DataArray
+        xarray DataArray containing model precipitation data in kg m-2 s-1
+
+    Returns
+    -------
+    precipitation_in_mm_per_day : xarray.DataArray
+        the input DataArray with precipitation units modified to mm day-1
+    """
+    # density 1000 kg m-3 => 1 kg m-2 == 1 mm
+    # There are 60*60*24 = 86400 seconds per day
+    precipitation_in_mm_per_day = precipitation_in_kg_per_m_squared_s * 86400
+
+    precipitation_in_mm_per_day.attrs["units"] = "mm/day"
+
+    if precipitation_in_mm_per_day.data.min() < 0.0:
+        raise ValueError("There is at least one negative precipitation value")
+    if precipitation_in_mm_per_day.data.max() > 2000:
+        raise ValueError("There is a precipitation value/s > 2000 mm/day")
+
+    return precipitation_in_mm_per_day
+
+
+def plot_zonally_averaged_precipitation(precipitation_data):
+    """
+    Plot zonally-averaged precipitation data and save to file.
+
+    Parameters
+    ----------
+    precipitation_data : xarray.DataArray
+        xarray DataSet containing precipitation model data, specifying precipitation in
+        [kg m-2 s-1] at given latitudes, longitudes and time. The Dataset should contain
+        four aligned DataArrays: precipitation, latitude, longitude and time.
+
+    Returns
+    -------
+    None
+
+    """
+    zonal_precipitation = precipitation_data["pr"].mean("lon", keep_attrs=True)
+
+    figure, axes = plt.subplots(nrows=4, ncols=1, figsize=(12, 8))
+
+    zonal_precipitation.sel(lat=[0]).plot.line(ax=axes[0], hue="lat")
+    zonal_precipitation.sel(lat=[-20, 20]).plot.line(ax=axes[1], hue="lat")
+    zonal_precipitation.sel(lat=[-45, 45]).plot.line(ax=axes[2], hue="lat")
+    zonal_precipitation.sel(lat=[-70, 70]).plot.line(ax=axes[3], hue="lat")
+
+    plt.tight_layout()
+    for axis in axes:
+        axis.set_ylim(0.0, 1.0e-4)
+        axis.grid()
+    plt.savefig("zonal.png", dpi=200)  # Save figure to file
+
+    figure, axes = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))
+
+    zonal_precipitation.T.plot()
+
+    plt.savefig("zonal_map.png", dpi=200)  # Save figure to file
+
+
+def get_country_annual_average(precipitation_data, countries):
+    """
+    Calculate annual precipitation averages for countries and save to file.
+
+    Parameters
+    ----------
+    precipitation_data : xarray.DataArray
+        xarray DataSet containing precipitation model data, specifying precipitation in
+        [kg m-2 s-1] at given latitudes, longitudes and time. The Dataset should contain
+        four aligned DataArrays: precipitation, latitude, longitude and time.
+    countries : dict(str: str)
+        dictionary mapping country names to regionmask codes. For a list see:
+        regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions
+
+    Returns
+    -------
+    None
+
+    """
+    annual_average_precipitation = (
+        precipitation_data["pr"].groupby("time.year").mean("time", keep_attrs=True)
+    )
+    annual_average_precipitation = convert_precipitation_units(
+        annual_average_precipitation
+    )
+
+    country_mask = regionmask.defined_regions.natural_earth_v5_0_0.countries_110.mask(
+        annual_average_precipitation
+    )
+
+    with open(
+        "annual_average_precipitation_by_country.txt", "w", encoding="utf-8"
+    ) as datafile:
+        for country_name, country_code in countries.items():
+            data_avg_mask = annual_average_precipitation.where(
+                country_mask.cf == country_code
+            )
+
+            for year in data_avg_mask.year.values:
+                precip = data_avg_mask.sel(year=year).mean().values
+                datafile.write(
+                    "{} {} : {:2.3f} mm/day\n".format(
+                        country_name.ljust(25), year, precip
+                    )
+                )
+            datafile.write("\n")
+
+
+def plot_enso_hovmoller_diagram(precipitation_data):
+    """
+    Plot Hovmöller diagram of equatorial precipitation to visualise ENSO.
+
+    Parameters
+    ----------
+    precipitation_data : xarray.DataArray
+        xarray DataSet containing precipitation model data, specifying precipitation in
+        [kg m-2 s-1] at given latitudes, longitudes and time. The Dataset should contain
+        four aligned DataArrays: precipitation, latitude, longitude and time.
+
+    Returns
+    -------
+    None
+
+    """
+    enso = (
+        precipitation_data["pr"]
+        .sel(lat=slice(-1, 1))
+        .sel(lon=slice(120, 280))
+        .mean(dim="lat", keep_attrs=True)
+    )
+
+    enso.plot()
+    plt.savefig("enso.png", dpi=200)  # Save figure to file
+
+
+def create_precipitation_climatology_plot(
+    seasonal_average_precipitation,
+    model_name,
+    season,
+    mask=None,
+    plot_gridlines=False,
+    levels=None,
+):
+    """
+    Plot the precipitation climatology.
+
+    Parameters
+    ----------
+    seasonal_average_precipitation : xarray.DataArray
+        Precipitation climatology data
+    model_name : str
+        Name of the climate model
+    season : str
+        Climatological season (one of DJF, MAM, JJA, SON)
+    mask : optional str
+        mask to apply to plot (one of "land" or "ocean")
+    plot_gridlines : bool
+        Select whether to plot gridlines
+    levels : list
+        Tick mark values for the colorbar
+
+    Returns
+    -------
+    None
+
+    """
+    if not levels:
+        levels = np.arange(0, 13.5, 1.5)
+
+    fig, geo_axes = plt.subplots(
+        nrows=1,
+        ncols=1,
+        figsize=(12, 5),
+        subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)},
+    )
+
+    seasonal_average_precipitation.sel(season=season).plot.contourf(
+        ax=geo_axes,
+        levels=levels,
+        extend="max",
+        transform=ccrs.PlateCarree(),
+        cbar_kwargs={"label": seasonal_average_precipitation.units},
+        cmap=cmocean.cm.rain,
+    )
+
+    geo_axes.add_feature(
+        cfeature.COASTLINE, lw=2
+    )  # Add coastines using cartopy feature
+
+    if mask:
+        # Mask out (fade) using 110m resolution data from cartopy.
+        geo_axes.add_feature(
+            cfeature.NaturalEarthFeature("physical", mask, "110m"),
+            ec=None,
+            fc="white",
+            lw=2,
+            alpha=0.75,
+        )
+
+    if plot_gridlines:
+        gridlines = geo_axes.gridlines(
+            crs=ccrs.PlateCarree(),
+            draw_labels=True,
+            linewidth=2,
+            color="gray",
+            alpha=0.5,
+            linestyle="--",
+        )
+        gridlines.top_labels = False
+        gridlines.left_labels = True
+        # gl.xlines = False
+        gridlines.xlocator = mticker.FixedLocator([-180, -90, 0, 90, 180])
+        gridlines.ylocator = mticker.FixedLocator(
+            [-66, -23, 0, 23, 66]
+        )  # Tropics & Polar Circles
+        gridlines.xformatter = LONGITUDE_FORMATTER
+        gridlines.yformatter = LATITUDE_FORMATTER
+        gridlines.xlabel_style = {"size": 15, "color": "gray"}
+        gridlines.ylabel_style = {"size": 15, "color": "gray"}
+
+    title = "{} precipitation climatology ({})".format(model_name, season)
+    plt.title(title)
+
+
+def main(
+    precipitation_netcdf_file,
+    season="DJF",
+    output_file="output.png",
+    plot_gridlines=False,
+    mask=None,
+    cbar_levels=None,
+    countries=None,
+):
+    """
+    Run the program for producing precipitation plots.
+
+    Parameters
+    ----------
+    precipitation_netcdf_file : str
+        netCDF filename to read precipitation data from
+    season : optional str
+        Climatological season (one of DJF, MAM, JJA, SON)
+    output_file : optional str
+        filename to save main image to
+    plot_gridlines : optional bool
+        Select whether to plot gridlines
+    mask : optional str
+        mask to apply to plot (one of "land" or "ocean")
+    cbar_levels : optional list
+        Tick mark values for the colorbar
+    countries : optional dict(str: str)
+        dictionary mapping country names to regionmask codes. For a list see:
+        regionmask.defined_regions.natural_earth_v5_0_0.countries_110.regions
+
+    Returns
+    -------
+    None
+
+    """
+    if countries is None:
+        countries = {"United Kingdom": "GB"}
+
+    precipitation_data = xr.open_dataset(precipitation_netcdf_file)
+
+    plot_zonally_averaged_precipitation(precipitation_data)
+    plot_enso_hovmoller_diagram(precipitation_data)
+    get_country_annual_average(precipitation_data, countries)
+
+    seasonal_average_precipitation = (
+        precipitation_data["pr"].groupby("time.season").mean("time", keep_attrs=True)
+    )
+
+    try:
+        input_units = seasonal_average_precipitation.attrs["units"]
+    except KeyError as exc:
+        raise KeyError(
+            "Precipitation variable in {pr_file} must have a units attribute"
+        ) from exc
+
+    if input_units == "kg m-2 s-1":
+        seasonal_average_precipitation = convert_precipitation_units(
+            seasonal_average_precipitation
+        )
+    elif input_units == "mm/day":
+        pass
+    else:
+        raise ValueError("""Input units are not 'kg m-2 s-1' or 'mm/day'""")
+
+    create_precipitation_climatology_plot(
+        seasonal_average_precipitation,
+        precipitation_data.attrs["source_id"],
+        season,
+        mask=mask,
+        plot_gridlines=plot_gridlines,
+        levels=cbar_levels,
+    )
+
+    plt.savefig(output_file, dpi=200)
+
+
+if __name__ == "__main__":
+    input_file = (
+        "../../data/pr_Amon_ACCESS-ESM1-5_historical_r1i1p1f1_gn_201001-201412.nc"
+    )
+    season_to_plot = "JJA"
+    output_filename = "output.png"
+    gridlines_on = True
+    mask_id = "ocean"
+    colorbar_levels = None
+    countries_to_record = {
+        "United Kingdom": "GB",
+        "United States of America": "US",
+        "Antarctica": "AQ",
+        "South Africa": "ZA",
+    }
+
+    main(
+        input_file,
+        season=season_to_plot,
+        mask=mask_id,
+        plot_gridlines=gridlines_on,
+        countries=countries_to_record,
+    )
diff --git a/slides/_comments_docstrings.qmd b/slides/_comments_docstrings.qmd
index 9b82a2c..2f82aa6 100644
--- a/slides/_comments_docstrings.qmd
+++ b/slides/_comments_docstrings.qmd
@@ -258,9 +258,9 @@ def calculate_gyroradius(mass, v_perp, charge, B, gamma=None):
 
 
 
-## Exercise 4 {.smaller}
+## Exercise 5 {.smaller}
 
-Go to exercise 4 and examine the comments:
+Go to exercise 5 and examine the comments:
 
 - Is there any dead code?
   - How is it best to handle it?
diff --git a/slides/_formatting_with_black.qmd b/slides/_formatting_with_black.qmd
new file mode 100644
index 0000000..c9f0bac
--- /dev/null
+++ b/slides/_formatting_with_black.qmd
@@ -0,0 +1,161 @@
+## Python PEPs {.smaller}
+
+[Python Enhancement Proposals](https://peps.python.org/)
+
+- Technical documentation for the python community
+- Guidelines, standards, and best-practice
+
+Relevant to us today are:
+
+* PEP8 - Python Style Guide [@PEP8]
+* PEP257 - Docstring Conventions [@PEP257]
+* PEP621 - Packaging [@PEP621]
+
+
+
+# PEP8 and Formatting
+
+## PEP8 & Formatting {.smaller}
+
+> _“Readability counts”_\
+> &emsp; &emsp; - Tim Peters in the [Zen of Python](https://peps.python.org/pep-0020/)
+
+By ensuring code aligns with PEP8 we:
+
+- standardise style,
+- conform to best-practices, and
+- improve code readability to
+- make code easier to share, and
+- reduce misinterpretation.
+
+::: {.fragment}
+"But I don't have time to read and memorise all of this..."
+:::
+
+
+## PEP8 & Formatting - Black {.smaller}
+
+![]( https://black.readthedocs.io/en/stable/_static/logo2-readme.png ){.absolute top=25% right=0% height=20%}
+
+Black [@black] - [black.readthedocs.io](https://black.readthedocs.io/en/stable/index.html)
+
+- a PEP 8 compliant formatter
+  - Strict subset of PEP8
+  - _"Opinionated so you don't have to be."_
+- For full details see [Black style](https://black.readthedocs.io/en/stable/the_black_code_style/current_style.html)
+- [Try online](https://black.vercel.app/)
+
+::: {.panel-tabset}
+
+### Linux/macOS
+
+```bash
+(myvenv) $ pip install black
+(myvenv) $ black myfile.py
+(myvenv) $ black mydirectory/
+```
+
+### Windows
+
+```powershell
+(myvenv) PS> pip install black
+(myvenv) PS> black myfile.py
+(myvenv) PS> black mydirectory/
+```
+:::
+
+::: {.notes}
+May look odd at first, but you soon get used to it.\
+Makes life so much easier after a while.
+:::
+
+
+## PEP8 & Formatting - Black - Example {.smaller}
+
+:::: {.columns}
+::: {.column width="50%"}
+```python
+def long_func(x, param_one, param_two=[], param_three=24, param_four=None,
+        param_five="Empty Report", param_six=123456):
+
+
+    val = 12*16 +(24) -10*param_one +  param_six
+
+    if x > 5:
+        
+        print("x is greater than 5")
+
+
+    else:
+        print("x is less than or equal to 5")
+
+
+    if param_four:
+        print(param_five)
+
+
+
+    print('You have called long_func.')
+    print("This function has several params.")
+
+    param_2.append(x*val)
+    return param_2
+
+```
+:::
+::: {.column}
+```python
+def long_func(
+    x,
+    param_one,
+    param_two=[],
+    param_three=24,
+    param_four=None,
+    param_five="Empty Report",
+    param_six=123456,
+):
+    val = 12 * 16 + (24) - 10 * param_one + param_six
+
+    if x > 5:
+        print("x is greater than 5")
+
+    else:
+        print("x is less than or equal to 5")
+
+    if param_four:
+        print(param_five)
+
+    print("You have called long_func.")
+    print("This function has several params.")
+
+    param_2.append(x * val)
+    return param_2
+```
+:::
+::::
+
+
+## PEP8 & Formatting - Black {.smaller}
+
+- I suggest incorporating into your projects now
+  - Well suited to incorporation into continuous integration of git hooks.
+  - Widely-used standard^[[Black - used by](https://github.com/psf/black#used-by)]
+
+- A version for jupyter notebooks exists.
+
+
+<!-- ------------------------------------------------------------------------------ -->
+
+
+## Exercise 2 {.smaller}
+
+Go to exercise 2 and:
+
+- install black
+- run black on `precipitation_climatology.py`
+- examine the output
+  - Is it more readable?
+  - Is there any aspect of the formatting style you find unintuitive?is it better?
+
+
+<!-- ------------------------------------------------------------------------------ -->
\ No newline at end of file
diff --git a/slides/_formatting.qmd b/slides/_formatting_with_pylint.qmd
similarity index 57%
rename from slides/_formatting.qmd
rename to slides/_formatting_with_pylint.qmd
index ac92298..3a4e439 100644
--- a/slides/_formatting.qmd
+++ b/slides/_formatting_with_pylint.qmd
@@ -1,165 +1,3 @@
-## Python PEPs {.smaller}
-
-[Python Enhancement Proposals](https://peps.python.org/)
-
-- Technical documentation for the python community
-- Guidelines, standards, and best-practice
-
-Relevant to us today are:
-
-* PEP8 - Python Style Guide [@PEP8]
-* PEP257 - Docstring Conventions [@PEP257]
-* PEP621 - Packaging [@PEP621]
-
-
-
-# PEP8 and Formatting
-
-## PEP8 & Formatting {.smaller}
-
-> _“Readability counts”_\
-> &emsp; &emsp; - Tim Peters in the [Zen of Python](https://peps.python.org/pep-0020/)
-
-By ensuring code aligns with PEP8 we:
-
-- standardise style,
-- conform to best-practices, and
-- improve code readability to
-- make code easier to share, and
-- reduce misinterpretation.
-
-::: {.fragment}
-"But I don't have time to read and memorise all of this..."
-:::
-
-
-## PEP8 & Formatting - Black {.smaller}
-
-![]( https://black.readthedocs.io/en/stable/_static/logo2-readme.png ){.absolute top=25% right=0% height=20%}
-
-Black [@black] - [black.readthedocs.io](https://black.readthedocs.io/en/stable/index.html)
-
-- a PEP 8 compliant formatter
-  - Strict subset of PEP8
-  - _"Opinionated so you don't have to be."_
-- For full details see [Black style](https://black.readthedocs.io/en/stable/the_black_code_style/current_style.html)
-- [Try online](https://black.vercel.app/)
-
-::: {.panel-tabset}
-
-### Linux/macOS
-
-```bash
-(myvenv) $ pip install black
-(myvenv) $ black myfile.py
-(myvenv) $ black mydirectory/
-```
-
-### Windows
-
-```powershell
-(myvenv) PS> pip install black
-(myvenv) PS> black myfile.py
-(myvenv) PS> black mydirectory/
-```
-:::
-
-::: {.notes}
-May look odd at first, but you soon get used to it.\
-Makes life so much easier after a while.
-:::
-
-
-## PEP8 & Formatting - Black - Example {.smaller}
-
-:::: {.columns}
-::: {.column width="50%"}
-```python
-def long_func(x, param_one, param_two=[], param_three=24, param_four=None,
-        param_five="Empty Report", param_six=123456):
-
-
-    val = 12*16 +(24) -10*param_one +  param_six
-
-    if x > 5:
-        
-        print("x is greater than 5")
-
-
-    else:
-        print("x is less than or equal to 5")
-
-
-    if param_four:
-        print(param_five)
-
-
-
-    print('You have called long_func.')
-    print("This function has several params.")
-
-    param_2.append(x*val)
-    return param_2
-
-```
-:::
-::: {.column}
-```python
-def long_func(
-    x,
-    param_one,
-    param_two=[],
-    param_three=24,
-    param_four=None,
-    param_five="Empty Report",
-    param_six=123456,
-):
-    val = 12 * 16 + (24) - 10 * param_one + param_six
-
-    if x > 5:
-        print("x is greater than 5")
-
-    else:
-        print("x is less than or equal to 5")
-
-    if param_four:
-        print(param_five)
-
-    print("You have called long_func.")
-    print("This function has several params.")
-
-    param_2.append(x * val)
-    return param_2
-```
-:::
-::::
-
-
-## PEP8 & Formatting - Black {.smaller}
-
-- I suggest incorporating into your projects now
-  - Well suited to incorporation into continuous integration of git hooks.
-  - Widely-used standard^[[Black - used by](https://github.com/psf/black#used-by)]
-
-- A version for jupyter notebooks exists.
-
-
-<!-- ------------------------------------------------------------------------------ -->
-
-
-## Exercise 2 {.smaller}
-
-Go to exercise 2 and:
-
-- install black
-- run black on `precipitation_climatology.py`
-- examine the output
-  - Is it more readable?
-  - Is there any aspect of the formatting style you find unintuitive?is it better?
-
-
-<!-- ------------------------------------------------------------------------------ -->
-
 ## PEP8 & Formatting - PyLint {.smaller}
 
 :::: {.columns}
@@ -331,9 +169,9 @@ Get a rating.
   - [emacs](https://www.emacswiki.org/emacs/PythonProgrammingInEmacs)
 
 
-## Exercise 3 {.smaller}
+## Exercise 4 {.smaller}
 
-Go to exercise 3 and:
+Go to exercise 4 and:
 
 - install pylint
 - run pylint on `precipitation_climatology.py`
diff --git a/slides/_fstrings_magic_config.qmd b/slides/_fstrings_magic_config.qmd
index ee6f232..1dbe3b8 100644
--- a/slides/_fstrings_magic_config.qmd
+++ b/slides/_fstrings_magic_config.qmd
@@ -28,6 +28,7 @@ print(f"a={a} and b={b}. Their product is {a * b}, sum is {a + b}, and a/b is {a
 ```
 
 See [Real Python](https://realpython.com/python-f-strings/) for more information.
+Note: pylint W1203 recommends against using f-strings in logging calls.
 
 
 ## Remove Magic Numbers {.smaller}
@@ -187,7 +188,7 @@ print(config)
 <!-- ------------------------------------------------------------------------------ -->
 
 
-## Exercise 5 {.smaller}
+## Exercise 6 {.smaller}
 
 :::: {.columns}
 
diff --git a/slides/_naming_for_clarity.qmd b/slides/_naming_for_clarity.qmd
new file mode 100644
index 0000000..73643af
--- /dev/null
+++ b/slides/_naming_for_clarity.qmd
@@ -0,0 +1,94 @@
+## Naming For Clarity {.smaller}
+
+::: {style="font-size: 90%;"}
+
+It may seem inconsequential, but carefully naming variables and methods can improve the
+readability of code massively and can help to make code self documenting.
+
+A few naming tips and conventions:
+
+::: {.incremental}
+- The name should show the intention, think about how someone else might read it (this could be future you)
+- Use pronounceable names e.g. `mass` not `ms`, `stem` not `stm`
+- avoid abbreviations and single letter variable names where possible
+- One word per concept e.g. choose one of `put`, `insert`, `add` in the same code base
+- Use names that can be searched
+- Describe content rather than storage type
+- Naming booleans, use prefixes like `is`, `has` or `can` and avoid negations like `not_green`
+- Plurals to indicate groups, e.g. a list of dog objects would be `dogs`, not `dog_list`
+- Keep it simple and use technical terms where appropriate
+- Use explaining variables
+:::
+
+:::
+
+## Naming For Clarity {.smaller}
+
+Some examples:
+
+- Pronounceable names without abbreviations: \
+```
+ ms    --> mass
+ chclt --> chocolate
+ stm   --> stem
+ ```
+- Naming for the content, not the type: \
+```
+array       --> dogs
+age_int     --> age
+country_set --> countries
+```
+- Naming booleans: \
+```
+purple    --> is_purple
+not_plant --> is_plant
+sidekick  --> has_sidekick
+```
+
+
+
+## Explaining Variables
+
+Without explaining variables, it is hard to see what this code is doing:
+
+```python
+import re
+
+re.match("^\\+?[1-9][0-9]{7,14}$", "+12223334444")
+```
+
+## With explaining variables:
+
+It is easier to see the intention. The code is more self-documenting.
+
+```python
+import re
+
+validate_phone_number_pattern = "^\\+?[1-9][0-9]{7,14}$"
+re.match(validate_phone_number_pattern, "+12223334444")
+```
+
+## Exercise 3 {.smaller}
+
+::: {style="font-size: 80%;"}
+
+Look through the code for any names of methods or variables that could be improved or
+clarified and update them. Note if you are using an IDE like Intellij or VSCode,
+you can use automatic renaming.
+Can you find an example from each of the suggestions listed below?
+Does this make the code easier to follow?
+
+Consider the following:
+
+- The name should show the intention, think about how someone else might read it (this could be future you)
+- Use pronounceable names e.g. `mass` not `ms`, `stem` not `stm`
+- avoid abbreviations and single letter variable names where possible
+- One word per concept e.g. choose one of `put`, `insert`, `add` in the same code base
+- Use names that can be searched
+- Describe content rather than storage type
+- Naming booleans, use prefixes like `is`, `has` or `can` and avoid negations like `not_green`
+- Plurals to indicate groups, e.g. a list of dog objects would be `dogs`, not `dog_list`
+- Keep it simple and use technical terms where appropriate
+- Use explaining variables
+
+:::
diff --git a/slides/python.qmd b/slides/python.qmd
index 99a71a5..41e9578 100644
--- a/slides/python.qmd
+++ b/slides/python.qmd
@@ -93,10 +93,17 @@ My background in this is domain science and numerical modelling on HPC.
 
 <!-- ------------------------------------------------------------------------------ -->
 
-{{< include _formatting.qmd >}}
+{{< include _formatting_with_black.qmd >}}
 
 <!-- ------------------------------------------------------------------------------ -->
 
+{{< include _naming_for_clarity.qmd >}}
+
+<!-- ------------------------------------------------------------------------------ -->
+
+{{< include _formatting_with_pylint.qmd >}}
+
+<!-- ------------------------------------------------------------------------------ -->
 {{< include _comments_docstrings.qmd >}}
 
 <!-- ------------------------------------------------------------------------------ -->