From 3836dc354c3fafaa42db3e375f6551e670fdcc5f Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Tue, 28 Apr 2026 16:05:59 +1000
Subject: [PATCH 01/50] Add an activity arg to cmip7_nitrogen_dirpath

---
 CMIP7/esm1p6/atmosphere/nitrogen/cmip7_HI_nitrogen_generate.py | 2 +-
 CMIP7/esm1p6/atmosphere/nitrogen/cmip7_PI_nitrogen_generate.py | 2 +-
 CMIP7/esm1p6/atmosphere/nitrogen/cmip7_nitrogen.py             | 3 ++-
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_HI_nitrogen_generate.py b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_HI_nitrogen_generate.py
index d4bcf92..f05e974 100644
--- a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_HI_nitrogen_generate.py
+++ b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_HI_nitrogen_generate.py
@@ -26,7 +26,7 @@ def parse_args():
 
 
 def cmip7_hi_nitrogen_filepath(args, species):
-    dirpath = cmip7_nitrogen_dirpath(args, "mon", species)
+    dirpath = cmip7_nitrogen_dirpath(args, "CMIP", "mon", species)
     filename = (
         f"{species}_input4MIPs_surfaceFluxes_CMIP_"
         f"{args.dataset_version}_gn_"
diff --git a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_PI_nitrogen_generate.py b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_PI_nitrogen_generate.py
index 9c33eef..56abb82 100644
--- a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_PI_nitrogen_generate.py
+++ b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_PI_nitrogen_generate.py
@@ -25,7 +25,7 @@ def parse_args():
 
 
 def cmip7_pi_nitrogen_filepath(args, species):
-    dirpath = cmip7_nitrogen_dirpath(args, "monC", species)
+    dirpath = cmip7_nitrogen_dirpath(args, "CMIP", "monC", species)
     filename = (
         f"{species}_input4MIPs_surfaceFluxes_CMIP_"
         f"{args.dataset_version}_gn_"
diff --git a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_nitrogen.py b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_nitrogen.py
index f3da718..f56b0ee 100644
--- a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_nitrogen.py
+++ b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_nitrogen.py
@@ -17,9 +17,10 @@
 NITROGEN_STASH_ITEM = 884
 
 
-def cmip7_nitrogen_dirpath(args, period, species):
+def cmip7_nitrogen_dirpath(args, activity, period, species):
     return (
         Path(args.cmip7_source_data_dirname)
+        / activity
         / "FZJ"
         / args.dataset_version
         / "atmos"

From 111be596e5ce9a377ee596f6bce9a89d35a989a5 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Tue, 28 Apr 2026 17:20:05 +1000
Subject: [PATCH 02/50] Generate ScenarioMIP nitrogen ancils

---
 .../nitrogen/cmip7_SM_nitrogen_generate.py    | 60 +++++++++++++++++++
 1 file changed, 60 insertions(+)
 create mode 100644 CMIP7/esm1p6/atmosphere/nitrogen/cmip7_SM_nitrogen_generate.py

diff --git a/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_SM_nitrogen_generate.py b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_SM_nitrogen_generate.py
new file mode 100644
index 0000000..120faec
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/nitrogen/cmip7_SM_nitrogen_generate.py
@@ -0,0 +1,60 @@
+from argparse import ArgumentParser
+from pathlib import Path
+
+from cmip7_ancil_argparse import common_parser
+from cmip7_ancil_common import extend_years
+from cmip7_ancil_constants import ANCIL_TODAY
+from nitrogen.cmip7_nitrogen import (
+    cmip7_nitrogen_dirpath,
+    load_cmip7_nitrogen,
+    regrid_cmip7_nitrogen,
+    save_cmip7_nitrogen,
+)
+
+
+def parse_args():
+    parser = ArgumentParser(
+        parents=[common_parser()],
+        prog="cmip7_SM_nitrogen_generate",
+        description=(
+            "Generate input files from CMIP7 ScenarioMIP nitrogen forcings"
+        ),
+    )
+    parser.add_argument("--scenario")
+    parser.add_argument("--dataset-date-range")
+    parser.add_argument("--save-filename")
+    return parser.parse_args()
+
+
+def cmip7_sm_nitrogen_filepath(args, species):
+    dirpath = cmip7_nitrogen_dirpath(args, "ScenarioMIP", "mon", species)
+    filename = (
+        f"{species}_input4MIPs_surfaceFluxes_ScenarioMIP_"
+        f"{args.dataset_version}_gn_"
+        f"{args.dataset_date_range}.nc"
+    )
+    return dirpath / filename
+
+
+def esm_sm_nitrogen_save_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "scenarios"
+        / args.scenario
+        / "atmosphere"
+        / "land"
+        / "biogeochemistry"
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    # Load the CMIP7 datasets
+    nitrogen_cube = load_cmip7_nitrogen(args, cmip7_sm_nitrogen_filepath)
+    # Regrid to match the ESM1.5 mask and extend the time series
+    esm_cube = extend_years(regrid_cmip7_nitrogen(args, nitrogen_cube))
+    # Save the ancillary
+    save_cmip7_nitrogen(args, esm_cube, esm_sm_nitrogen_save_dirpath)

From 41d61177b47b564ccf8b12b64a8b01648c2171d1 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 4 May 2026 12:33:28 +1000
Subject: [PATCH 03/50] Generate Ozone ancillaries for ScenarioMIP

---
 .../ozone/cmip7_SM_ozone_generate.py          | 54 +++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 CMIP7/esm1p6/atmosphere/ozone/cmip7_SM_ozone_generate.py

diff --git a/CMIP7/esm1p6/atmosphere/ozone/cmip7_SM_ozone_generate.py b/CMIP7/esm1p6/atmosphere/ozone/cmip7_SM_ozone_generate.py
new file mode 100644
index 0000000..9afb9eb
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/ozone/cmip7_SM_ozone_generate.py
@@ -0,0 +1,54 @@
+from argparse import ArgumentParser
+from pathlib import Path
+
+from cmip7_ancil_argparse import (
+    grid_parser,
+    path_parser,
+)
+from cmip7_ancil_common import save_ancil
+from cmip7_ancil_constants import ANCIL_TODAY
+from ozone.cmip7_ozone import (
+    fix_cmip7_ozone,
+    load_cmip7_ozone,
+    ozone_parser,
+)
+
+
+def parse_args():
+    parser = ArgumentParser(
+        parents=[path_parser(), grid_parser(), ozone_parser()],
+        prog="cmip7_SM_ozone_generate",
+        description=(
+            "Generate input files from UK CMIP7 ScenarioMIP ozone forcings"
+        ),
+    )
+    parser.add_argument("--scenario")
+    return parser.parse_args()
+
+
+def esm_sm_ozone_save_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "scenarios"
+        / args.scenario
+        / "forcing"
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )
+
+
+def save_cmip7_sm_ozone(args, cube):
+    # Save as an ancillary file
+    save_dirpath = esm_sm_ozone_save_dirpath(args)
+    save_ancil(cube, save_dirpath, args.save_filename, replace_bounds=True)
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    # Load the CMIP7 datasets
+    ozone_cube = load_cmip7_ozone(args)
+    # Match the ESM1.5 mask
+    esm_cube = fix_cmip7_ozone(args, ozone_cube)
+    # Save the ancillary
+    save_cmip7_sm_ozone(args, esm_cube)

From 3f7bcf529035c36227fe2567a4ab556cf4a0cdfc Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Tue, 26 May 2026 15:28:59 +1000
Subject: [PATCH 04/50] Refactor volcanic forcings and add ScenarioMIP

---
 CMIP7/esm1p6/atmosphere/cmip7_SM.py           |  20 +++
 .../volcanic/cmip7_HI_volcanic_generate.py    | 158 ++----------------
 .../volcanic/cmip7_PI_volcanic_generate.py    |   6 +-
 .../volcanic/cmip7_SM_volcanic_generate.py    |  59 +++++++
 .../atmosphere/volcanic/cmip7_volcanic.py     | 152 ++++++++++++++++-
 5 files changed, 244 insertions(+), 151 deletions(-)
 create mode 100644 CMIP7/esm1p6/atmosphere/cmip7_SM.py
 create mode 100644 CMIP7/esm1p6/atmosphere/volcanic/cmip7_SM_volcanic_generate.py

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_SM.py b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
new file mode 100644
index 0000000..821fd36
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
@@ -0,0 +1,20 @@
+from pathlib import Path
+
+from cmip7_ancil_constants import ANCIL_TODAY
+
+CMIP7_SM_BEG_YEAR = 2025
+# Model time interpolation requires an extra year
+CMIP7_SM_END_YEAR = 2100
+CMIP7_SM_NBR_YEARS = CMIP7_SM_END_YEAR + 1 - CMIP7_SM_BEG_YEAR
+
+
+def esm_sm_forcing_save_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "scenarios"
+        / args.scenario
+        / "atmosphere"
+        / "forcing"
+        / "resolution_independent"
+        / ANCIL_TODAY
+    )
diff --git a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_HI_volcanic_generate.py b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_HI_volcanic_generate.py
index 01f6902..21acb48 100644
--- a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_HI_volcanic_generate.py
+++ b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_HI_volcanic_generate.py
@@ -1,25 +1,14 @@
 from argparse import ArgumentParser
 
-import cftime
-import iris
-import numpy as np
 from cmip7_ancil_argparse import dataset_parser, path_parser
-from cmip7_ancil_constants import MONTHS_IN_A_YEAR
 from cmip7_HI import esm_hi_forcing_save_dirpath
 from volcanic.cmip7_volcanic import (
-    SAOD_WAVELENGTH,
     cmip7_volcanic_dirpath,
-    constrain_to_wavelength,
-    mean_over_latitudes,
-    sum_over_height_layers,
+    save_stratospheric_aerosol_optical_depth,
 )
 
 CMIP7_HI_VOLCANIC_BEG_YEAR = 1850
 CMIP7_HI_VOLCANIC_END_YEAR = 2023
-CMIP7_HI_SAOD_TAPER_END_YEAR = 2033
-CMIP7_HI_SAOD_ARRAY_END_YEAR = 2300
-NBR_TAPER_YEARS = CMIP7_HI_SAOD_TAPER_END_YEAR - CMIP7_HI_VOLCANIC_END_YEAR
-NBR_OF_BANDS = 4
 
 
 def parse_args():
@@ -43,152 +32,27 @@ def cmip7_hi_volcanic_filename(args):
     )
 
 
-def constrain_to_year_month(cube, year, month):
-    """
-    Constrain to a given year and month. See #iris.coords.Cell.point in
-    scitools-iris.readthedocs.io/en/stable/generated/api/iris.coords.html
-    """
-    calendar = "proleptic_gregorian"
-    beg_date = cftime.datetime(year, month, 1, calendar=calendar)
-    end_year = year + 1 if month == MONTHS_IN_A_YEAR else year
-    end_month = 1 if month == MONTHS_IN_A_YEAR else month + 1
-    end_date = cftime.datetime(end_year, end_month, 1, calendar=calendar)
-    ym_constraint = iris.Constraint(
-        time=lambda cell: beg_date <= cell.point < end_date
-    )
-    return cube.extract(ym_constraint)
-
-
-def constrain_to_latitude_band(cube, band):
-    """
-    Constrain to one of four equal latitude bands.
-    """
-    lat_bound = [90, 30, 0, -30, -90]
-    lat_constraint = iris.Constraint(
-        latitude=(
-            lambda cell: lat_bound[band] > cell.point >= lat_bound[band + 1]
-        )
-    )
-    return cube.extract(lat_constraint)
-
-
-def taper_saod(saod_for_beg_year, saod_for_end_year):
-    """
-    Interpolate between the saod values in saod_for_beg_year
-    and saod_for_end_year. The SAOD values taper from saod_for_end_year
-    towards saod_for_beg_year until CMIP7_HI_SAOD_TAPER_END_YEAR,
-    and remain at saod_for_beg_year afterwards.
-    """
-    RATIO_ARRAY_LEN = CMIP7_HI_SAOD_ARRAY_END_YEAR - CMIP7_HI_VOLCANIC_END_YEAR
-    saod_array = np.zeros((RATIO_ARRAY_LEN, MONTHS_IN_A_YEAR, NBR_OF_BANDS))
-    ratio_array = np.zeros(RATIO_ARRAY_LEN)
-    for index in range(RATIO_ARRAY_LEN):
-        ratio_array[index] = (index + 1) / float(NBR_TAPER_YEARS)
-    ratio_endpoints = np.array([0.0, 1.0])
-    for month_m1 in range(MONTHS_IN_A_YEAR):
-        # Divide into latitude bands.
-        for lat_band_nbr in range(NBR_OF_BANDS):
-            saod_beg = saod_for_beg_year[month_m1, lat_band_nbr]
-            saod_end = saod_for_end_year[month_m1, lat_band_nbr]
-            saod_endpoints = np.array([saod_end, saod_beg])
-            saod_array[:, month_m1, lat_band_nbr] = np.interp(
-                ratio_array, ratio_endpoints, saod_endpoints
-            )
-    return saod_array
-
-
-def save_hi_year_tapered_saod(year, tapered_saod_array, save_file):
-    """
-    Save one year's worth of interpolated saod values in save_file.
-    """
-    index = year - (CMIP7_HI_VOLCANIC_END_YEAR + 1)
-    for month in range(1, MONTHS_IN_A_YEAR + 1):
-        print(f"{year:4d} {month:4d}", end="", file=save_file)
-
-        # Divide into latitude bands.
-        for lat_band_nbr in range(NBR_OF_BANDS):
-            saod = tapered_saod_array[index, month - 1, lat_band_nbr]
-            print(
-                f"{saod:7.1f}",
-                end="",
-                file=save_file,
-            )
-        print(file=save_file)
-
-
 def save_hi_stratospheric_aerosol_optical_depth(args, dataset_path):
     """
     Calculate the average stratospheric aerosol optical depth (SAOD)
     for each historical month by averaging extinction over latitude,
     and summing over stratospheric layers. Save to the save file.
     """
-    # Load the dataset into an Iris cube.
-    cube = iris.load_cube(dataset_path)
-
-    # Constrain to just the CMIP7 prescribed wavelength.
-    cube = constrain_to_wavelength(cube, SAOD_WAVELENGTH)
-
-    # Replace NaN values with 0.
-    np.nan_to_num(cube.data, copy=False)
-
-    save_dirpath = esm_hi_forcing_save_dirpath(args)
-    # Ensure that the save directory exists.
-    save_dirpath.mkdir(mode=0o755, parents=True, exist_ok=True)
-    save_filepath = save_dirpath / args.save_filename
-    # Keep the BEG_YEAR and END_YEAR SAOD values in arrays.
-    saod_for_beg_year = np.zeros((MONTHS_IN_A_YEAR, NBR_OF_BANDS))
-    saod_for_end_year = np.zeros((MONTHS_IN_A_YEAR, NBR_OF_BANDS))
-    with open(save_filepath, "w") as save_file:
-        # Iterate over years and months.
-        for year in range(
-            CMIP7_HI_VOLCANIC_BEG_YEAR, CMIP7_HI_VOLCANIC_END_YEAR + 1
-        ):
-            for month in range(1, MONTHS_IN_A_YEAR + 1):
-                print(f"{year:4d} {month:4d}", end="", file=save_file)
-                ym_cube = constrain_to_year_month(cube, year, month)
-
-                # Divide into latitude bands.
-                for lat_band_nbr in range(NBR_OF_BANDS):
-                    lat_cube = constrain_to_latitude_band(ym_cube, lat_band_nbr)
-
-                    # Find the mean over all latitudes included in this band,
-                    # weighted by area.
-                    lat_cube = mean_over_latitudes(lat_cube)
-
-                    # Calculate the stratospheric aerosol optical depth
-                    # by summing over stratospheric layers,
-                    # weighted by layer height.
-                    lat_cube = sum_over_height_layers(lat_cube)
-                    saod = lat_cube.data * 10000.0
-                    print(
-                        f"{saod:7.1f}",
-                        end="",
-                        file=save_file,
-                    )
-                    # Save the SAOD values for CMIP7_HI_VOLCANIC_BEG_YEAR.
-                    if year == CMIP7_HI_VOLCANIC_BEG_YEAR:
-                        saod_for_beg_year[month - 1, lat_band_nbr] = saod
-                    # Save the SAOD values for CMIP7_HI_VOLCANIC_END_YEAR.
-                    if year == CMIP7_HI_VOLCANIC_END_YEAR:
-                        saod_for_end_year[month - 1, lat_band_nbr] = saod
-                print(file=save_file)
-        # For years from CMIP7_HI_VOLCANIC_END_YEAR + 1 to
-        # CMIP7_HI_SAOD_ARRAY_END_YEAR interpolate between the saod values
-        # in saod_for_beg_year and saod_for_end_year and save values in
-        # save_file.
-        tapered_saod_array = taper_saod(saod_for_beg_year, saod_for_end_year)
-        for year in range(
-            CMIP7_HI_VOLCANIC_END_YEAR + 1, CMIP7_HI_SAOD_ARRAY_END_YEAR + 1
-        ):
-            save_hi_year_tapered_saod(year, tapered_saod_array, save_file)
+    save_stratospheric_aerosol_optical_depth(
+        args,
+        CMIP7_HI_VOLCANIC_BEG_YEAR,
+        CMIP7_HI_VOLCANIC_END_YEAR,
+        dataset_path,
+        esm_hi_forcing_save_dirpath(args),
+    )
 
 
 if __name__ == "__main__":
     args = parse_args()
 
-    dataset_path = cmip7_volcanic_dirpath(
-        args, period="mon"
-    ) / cmip7_hi_volcanic_filename(args)
+    dirpath = cmip7_volcanic_dirpath(args, "CMIP", "mon")
+    filename = cmip7_hi_volcanic_filename(args)
+    dataset_path = dirpath / filename
 
     # Calculate and save the average stratospheric aerosol optical depth.
     save_hi_stratospheric_aerosol_optical_depth(args, dataset_path)
diff --git a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_PI_volcanic_generate.py b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_PI_volcanic_generate.py
index e600050..54d96d9 100644
--- a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_PI_volcanic_generate.py
+++ b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_PI_volcanic_generate.py
@@ -113,9 +113,9 @@ def cmip7_pi_volcanic_patch(average_saod):
 if __name__ == "__main__":
     args = parse_args()
 
-    dataset_path = cmip7_volcanic_dirpath(
-        args, period="monC"
-    ) / cmip7_pi_volcanic_filename(args)
+    dirpath = cmip7_volcanic_dirpath(args, "CMIP", "monC")
+    filename = cmip7_pi_volcanic_filename(args)
+    dataset_path = dirpath / filename
 
     # Calculate the average stratospheric optical depth.
     average_saod = average_stratospheric_aerosol_optical_depth(dataset_path)
diff --git a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_SM_volcanic_generate.py b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_SM_volcanic_generate.py
new file mode 100644
index 0000000..4ff1587
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_SM_volcanic_generate.py
@@ -0,0 +1,59 @@
+from argparse import ArgumentParser
+
+from cmip7_ancil_argparse import dataset_parser, path_parser
+from cmip7_SM import esm_sm_forcing_save_dirpath
+from volcanic.cmip7_volcanic import (
+    cmip7_volcanic_dirpath,
+    save_stratospheric_aerosol_optical_depth,
+)
+
+CMIP7_SM_VOLCANIC_BEG_YEAR = 2022
+CMIP7_SM_VOLCANIC_END_YEAR = 2100
+
+
+def parse_args():
+    parser = ArgumentParser(
+        prog="cmip7_SM_volcanic_generate",
+        description=(
+            "Generate input files from CMIP7 ScenarioMIP volcanic forcings"
+        ),
+        parents=[path_parser(), dataset_parser()],
+    )
+    parser.add_argument("--scenario")
+    parser.add_argument("--dataset-date-range")
+    parser.add_argument("--save-filename")
+    return parser.parse_args()
+
+
+def cmip7_sm_volcanic_filename(args):
+    return (
+        f"ext_input4MIPs_aerosolProperties_ScenarioMIP_"
+        f"{args.dataset_version}_gnz_"
+        f"{args.dataset_date_range}.nc"
+    )
+
+
+def save_sm_stratospheric_aerosol_optical_depth(args, dataset_path):
+    """
+    Calculate the average stratospheric aerosol optical depth (SAOD)
+    for each historical month by averaging extinction over latitude,
+    and summing over stratospheric layers. Save to the save file.
+    """
+    save_stratospheric_aerosol_optical_depth(
+        args,
+        CMIP7_SM_VOLCANIC_BEG_YEAR,
+        CMIP7_SM_VOLCANIC_END_YEAR,
+        dataset_path,
+        esm_sm_forcing_save_dirpath(args),
+    )
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    dirpath = cmip7_volcanic_dirpath(args, "ScenarioMIP", "mon")
+    filename = cmip7_sm_volcanic_filename(args)
+    dataset_path = dirpath / filename
+
+    # Calculate and save the average stratospheric aerosol optical depth.
+    save_sm_stratospheric_aerosol_optical_depth(args, dataset_path)
diff --git a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py
index caeb046..9125af1 100644
--- a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py
+++ b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py
@@ -1,15 +1,21 @@
 from pathlib import Path
 
+import cftime
 import iris
 import numpy as np
+from cmip7_ancil_constants import MONTHS_IN_A_YEAR
 
+NBR_OF_BANDS = 4
+NBR_TAPER_YEARS = 10
+SAOD_ARRAY_END_YEAR = 2300
 # The prescribed wavelength for stratospheric aerosol optical depth
 SAOD_WAVELENGTH = 550.0 * 1e-9
 
 
-def cmip7_volcanic_dirpath(args, period):
+def cmip7_volcanic_dirpath(args, activity, period):
     return (
         Path(args.cmip7_source_data_dirname)
+        / activity
         / "uoexeter"
         / args.dataset_version
         / "atmos"
@@ -52,3 +58,147 @@ def sum_over_height_layers(cube):
         iris.analysis.SUM,
         weights=height_weights,
     )
+
+
+def constrain_to_year_month(cube, year, month):
+    """
+    Constrain to a given year and month. See #iris.coords.Cell.point in
+    scitools-iris.readthedocs.io/en/stable/generated/api/iris.coords.html
+    """
+    calendar = "proleptic_gregorian"
+    beg_date = cftime.datetime(year, month, 1, calendar=calendar)
+    end_year = year + 1 if month == MONTHS_IN_A_YEAR else year
+    end_month = 1 if month == MONTHS_IN_A_YEAR else month + 1
+    end_date = cftime.datetime(end_year, end_month, 1, calendar=calendar)
+    ym_constraint = iris.Constraint(
+        time=lambda cell: beg_date <= cell.point < end_date
+    )
+    return cube.extract(ym_constraint)
+
+
+def constrain_to_latitude_band(cube, band):
+    """
+    Constrain to one of four equal latitude bands.
+    """
+    lat_bound = [90, 30, 0, -30, -90]
+    lat_constraint = iris.Constraint(
+        latitude=(
+            lambda cell: lat_bound[band] > cell.point >= lat_bound[band + 1]
+        )
+    )
+    return cube.extract(lat_constraint)
+
+
+def taper_saod(
+    volcanic_end_year,
+    saod_for_beg_year,
+    saod_for_end_year,
+):
+    """
+    Interpolate between the saod values in saod_for_beg_year
+    and saod_for_end_year. The SAOD values taper from saod_for_end_year
+    towards saod_for_beg_year for NBR_TAPER_YEARS, and remain at
+    saod_for_beg_year afterwards.
+    """
+    RATIO_ARRAY_LEN = SAOD_ARRAY_END_YEAR - volcanic_end_year
+    saod_array = np.zeros((RATIO_ARRAY_LEN, MONTHS_IN_A_YEAR, NBR_OF_BANDS))
+    ratio_array = np.zeros(RATIO_ARRAY_LEN)
+    for index in range(RATIO_ARRAY_LEN):
+        ratio_array[index] = (index + 1) / float(NBR_TAPER_YEARS)
+    ratio_endpoints = np.array([0.0, 1.0])
+    for month_m1 in range(MONTHS_IN_A_YEAR):
+        # Divide into latitude bands.
+        for lat_band_nbr in range(NBR_OF_BANDS):
+            saod_beg = saod_for_beg_year[month_m1, lat_band_nbr]
+            saod_end = saod_for_end_year[month_m1, lat_band_nbr]
+            saod_endpoints = np.array([saod_end, saod_beg])
+            saod_array[:, month_m1, lat_band_nbr] = np.interp(
+                ratio_array, ratio_endpoints, saod_endpoints
+            )
+    return saod_array
+
+
+def save_year_tapered_saod(
+    year, tapered_saod_array, volcanic_end_year, save_file
+):
+    """
+    Save one year's worth of interpolated saod values in save_file.
+    """
+    index = year - (volcanic_end_year + 1)
+    for month in range(1, MONTHS_IN_A_YEAR + 1):
+        print(f"{year:4d} {month:4d}", end="", file=save_file)
+
+        # Divide into latitude bands.
+        for lat_band_nbr in range(NBR_OF_BANDS):
+            saod = tapered_saod_array[index, month - 1, lat_band_nbr]
+            print(
+                f"{saod:7.1f}",
+                end="",
+                file=save_file,
+            )
+        print(file=save_file)
+
+
+def save_stratospheric_aerosol_optical_depth(
+    args, volcanic_beg_year, volcanic_end_year, dataset_path, save_dirpath
+):
+    """
+    Calculate the average stratospheric aerosol optical depth (SAOD)
+    for each historical month by averaging extinction over latitude,
+    and summing over stratospheric layers. Save to the save file.
+    """
+    # Load the dataset into an Iris cube.
+    cube = iris.load_cube(dataset_path)
+
+    # Constrain to just the CMIP7 prescribed wavelength.
+    cube = constrain_to_wavelength(cube, SAOD_WAVELENGTH)
+
+    # Replace NaN values with 0.
+    np.nan_to_num(cube.data, copy=False)
+
+    # Ensure that the save directory exists.
+    save_dirpath.mkdir(mode=0o755, parents=True, exist_ok=True)
+    save_filepath = save_dirpath / args.save_filename
+    # Keep the BEG_YEAR and END_YEAR SAOD values in arrays.
+    saod_for_beg_year = np.zeros((MONTHS_IN_A_YEAR, NBR_OF_BANDS))
+    saod_for_end_year = np.zeros((MONTHS_IN_A_YEAR, NBR_OF_BANDS))
+    with open(save_filepath, "w") as save_file:
+        # Iterate over years and months.
+        for year in range(volcanic_beg_year, volcanic_end_year + 1):
+            for month in range(1, MONTHS_IN_A_YEAR + 1):
+                print(f"{year:4d} {month:4d}", end="", file=save_file)
+                ym_cube = constrain_to_year_month(cube, year, month)
+
+                # Divide into latitude bands.
+                for lat_band_nbr in range(NBR_OF_BANDS):
+                    lat_cube = constrain_to_latitude_band(ym_cube, lat_band_nbr)
+
+                    # Find the mean over all latitudes included in this band,
+                    # weighted by area.
+                    lat_cube = mean_over_latitudes(lat_cube)
+
+                    # Calculate the stratospheric aerosol optical depth
+                    # by summing over stratospheric layers,
+                    # weighted by layer height.
+                    lat_cube = sum_over_height_layers(lat_cube)
+                    saod = lat_cube.data * 10000.0
+                    print(
+                        f"{saod:7.1f}",
+                        end="",
+                        file=save_file,
+                    )
+                    # Save the SAOD values for volcanic_beg_year.
+                    if year == volcanic_beg_year:
+                        saod_for_beg_year[month - 1, lat_band_nbr] = saod
+                    # Save the SAOD values for volcanic_end_year.
+                    if year == volcanic_end_year:
+                        saod_for_end_year[month - 1, lat_band_nbr] = saod
+                print(file=save_file)
+        # For years from volcanic_end_year + 1 to SAOD_ARRAY_END_YEAR
+        # interpolate between the saod values in saod_for_beg_year and
+        # saod_for_end_year and save values in save_file.
+        tapered_saod_array = taper_saod(saod_for_beg_year, saod_for_end_year)
+        for year in range(volcanic_end_year + 1, SAOD_ARRAY_END_YEAR + 1):
+            save_year_tapered_saod(
+                year, tapered_saod_array, volcanic_end_year, save_file
+            )

From 6b7b5c310efa7b8ad3da652ce56417b8e4d18f71 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Tue, 26 May 2026 17:08:52 +1000
Subject: [PATCH 05/50] Fix the call to taper_saod

---
 CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py
index 9125af1..ae6e8c9 100644
--- a/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py
+++ b/CMIP7/esm1p6/atmosphere/volcanic/cmip7_volcanic.py
@@ -197,7 +197,9 @@ def save_stratospheric_aerosol_optical_depth(
         # For years from volcanic_end_year + 1 to SAOD_ARRAY_END_YEAR
         # interpolate between the saod values in saod_for_beg_year and
         # saod_for_end_year and save values in save_file.
-        tapered_saod_array = taper_saod(saod_for_beg_year, saod_for_end_year)
+        tapered_saod_array = taper_saod(
+            volcanic_end_year, saod_for_beg_year, saod_for_end_year
+        )
         for year in range(volcanic_end_year + 1, SAOD_ARRAY_END_YEAR + 1):
             save_year_tapered_saod(
                 year, tapered_saod_array, volcanic_end_year, save_file

From c0e1ebea21fe1c0704273e1aac6e265995ae3ce3 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 28 May 2026 17:25:12 +1000
Subject: [PATCH 06/50] Refactor solar forcings and add ScenarioMIP

---
 .../solar/cmip7_HI_solar_generate.py          | 64 ++++---------------
 .../solar/cmip7_PI_solar_generate.py          |  9 ++-
 .../solar/cmip7_SM_solar_generate.py          | 58 +++++++++++++++++
 CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py  | 56 +++++++++++++++-
 4 files changed, 129 insertions(+), 58 deletions(-)
 create mode 100644 CMIP7/esm1p6/atmosphere/solar/cmip7_SM_solar_generate.py

diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_HI_solar_generate.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_HI_solar_generate.py
index a926d89..b0f5bb7 100644
--- a/CMIP7/esm1p6/atmosphere/solar/cmip7_HI_solar_generate.py
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_HI_solar_generate.py
@@ -1,19 +1,16 @@
 from argparse import ArgumentParser
 
-import iris
-import numpy as np
 from cmip7_ancil_argparse import common_parser
-from cmip7_ancil_constants import REAL_MISSING_DATA_INDICATOR
 from cmip7_HI import (
     CMIP7_HI_BEG_YEAR,
     CMIP7_HI_END_YEAR,
     esm_hi_forcing_save_dirpath,
 )
-from cmip7_PI import CMIP7_PI_YEAR
-from solar.cmip7_solar import cmip7_solar_dirpath, load_cmip7_solar_cube
-
-SOLAR_ARRAY_BEG_YEAR = 1700
-SOLAR_ARRAY_END_YEAR = 2300
+from solar.cmip7_solar import (
+    cmip7_solar_dirpath,
+    cmip7_solar_save,
+    load_cmip7_solar_cube,
+)
 
 
 def parse_args():
@@ -31,64 +28,27 @@ def parse_args():
     return parser.parse_args()
 
 
-def cmip7_hi_solar_year_mean(cube):
-    """
-    Calculate mean TSI values for each year and save them into an array.
-    """
-    NBR_YEARS = SOLAR_ARRAY_END_YEAR - SOLAR_ARRAY_BEG_YEAR + 1
-    solar_array = np.zeros(NBR_YEARS)
-    # Calculate and save the mean annual TSI for each CMIP7 historical year.
-    historical_year_range = range(CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR + 1)
-    assert CMIP7_PI_YEAR in historical_year_range
-    for year in historical_year_range:
-        year_cons = iris.Constraint(time=lambda cell: cell.point.year == year)
-        year_cube = cube.extract(year_cons)
-        year_mean = year_cube.collapsed("time", iris.analysis.MEAN).data
-        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = year_mean
-        # Save the year mean for the pre-industrial year.
-        if year == CMIP7_PI_YEAR:
-            pi_year_mean = year_mean
-
-    # For the years from SOLAR_ARRAY_BEG_YEAR to CMIP7_HI_BEG_YEAR - 1,
-    # set the saved TSI value to the pre-industrial year mean TSI.
-    for year in range(SOLAR_ARRAY_BEG_YEAR, CMIP7_HI_BEG_YEAR):
-        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = pi_year_mean
-
-    # For the years from CMIP7_HI_END_YEAR + 1 to SOLAR_ARRAY_END_YEAR,
-    # set the saved TSI value to the real missing data indicator
-    for year in range(CMIP7_HI_END_YEAR + 1, SOLAR_ARRAY_END_YEAR + 1):
-        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = REAL_MISSING_DATA_INDICATOR
-    return solar_array
-
-
 def cmip7_hi_solar_save(args, cube):
     """
     Save the TSI values for each year into a text file.
     """
-    solar_array = cmip7_hi_solar_year_mean(cube)
     save_dirpath = esm_hi_forcing_save_dirpath(args)
-    # Ensure that the save directory exists.
-    save_dirpath.mkdir(mode=0o755, parents=True, exist_ok=True)
-    save_filepath = save_dirpath / args.save_filename
-    with open(save_filepath, "w") as save_file:
-        for year in range(SOLAR_ARRAY_BEG_YEAR, SOLAR_ARRAY_END_YEAR + 1):
-            year_mean = solar_array[year - SOLAR_ARRAY_BEG_YEAR]
-            if year_mean == REAL_MISSING_DATA_INDICATOR:
-                print(year, f"{year_mean:.1f}", file=save_file)
-            else:
-                print(year, f"{year_mean:.3f}", file=save_file)
+    cmip7_solar_save(
+        args, cube, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR, save_dirpath
+    )
 
 
 if __name__ == "__main__":
     args = parse_args()
 
-    cmip7_filename = (
+    dirpath = cmip7_solar_dirpath(args, "CMIP", "mon")
+    filename = (
         "multiple_input4MIPs_solar_CMIP_"
         f"{args.dataset_version}_gn_"
         f"{args.dataset_date_range}.nc"
     )
-    cmip7_filepath = cmip7_solar_dirpath(args, "mon") / cmip7_filename
+    dataset_path = dirpath / filename
 
-    solar_irradiance_cube = load_cmip7_solar_cube(cmip7_filepath)
+    solar_irradiance_cube = load_cmip7_solar_cube(dataset_path)
 
     cmip7_hi_solar_save(args, solar_irradiance_cube)
diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py
index 88141f4..2c433fd 100644
--- a/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py
@@ -48,12 +48,11 @@ def cmip7_pi_solar_patch(solar_irradiance):
 if __name__ == "__main__":
     args = parse_args()
 
-    cmip7_filename = (
-        f"multiple_input4MIPs_solar_CMIP_{args.dataset_version}_gn.nc"
-    )
-    cmip7_filepath = cmip7_solar_dirpath(args, "fx") / cmip7_filename
+    dirpath = cmip7_solar_dirpath(args, "ScenarioMIP", "fx")
+    filename = f"multiple_input4MIPs_solar_CMIP_{args.dataset_version}_gn.nc"
+    dataset_path = dirpath / filename
 
-    solar_irradiance_cube = load_cmip7_solar_cube(cmip7_filepath)
+    solar_irradiance_cube = load_cmip7_solar_cube(dataset_path)
     solar_irradiance = solar_irradiance_cube[0].data
 
     # Patch the SC variable in the coupling namelist
diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_SM_solar_generate.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_SM_solar_generate.py
new file mode 100644
index 0000000..c6037e9
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_SM_solar_generate.py
@@ -0,0 +1,58 @@
+from argparse import ArgumentParser
+
+from cmip7_ancil_argparse import common_parser
+from cmip7_SM import esm_sm_forcing_save_dirpath
+from solar.cmip7_solar import (
+    cmip7_solar_dirpath,
+    cmip7_solar_save,
+    load_cmip7_solar_cube,
+)
+
+CMIP7_SM_SOLAR_BEG_YEAR = 2022
+CMIP7_SM_SOLAR_END_YEAR = 2299
+
+
+def parse_args():
+    parser = ArgumentParser(
+        prog="cmip7_SM_solar_generate",
+        description=(
+            "Generate input files from CMIP7 ScenarioMIP solar forcings"
+        ),
+        parents=[
+            common_parser(),
+        ],
+    )
+    parser.add_argument("--scenario")
+    parser.add_argument("--dataset-date-range")
+    parser.add_argument("--save-filename")
+    return parser.parse_args()
+
+
+def cmip7_sm_solar_save(args, cube):
+    """
+    Save the TSI values for each year into a text file.
+    """
+    save_dirpath = esm_sm_forcing_save_dirpath(args)
+    cmip7_solar_save(
+        args,
+        cube,
+        CMIP7_SM_SOLAR_BEG_YEAR,
+        CMIP7_SM_SOLAR_END_YEAR,
+        save_dirpath,
+    )
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    dirpath = cmip7_solar_dirpath(args, "ScenarioMIP", "mon")
+    filename = (
+        "multiple_input4MIPs_solar_ScenarioMIP_"
+        f"{args.dataset_version}_gn_"
+        f"{args.dataset_date_range}.nc"
+    )
+    dataset_path = dirpath / filename
+
+    solar_irradiance_cube = load_cmip7_solar_cube(dataset_path)
+
+    cmip7_sm_solar_save(args, solar_irradiance_cube)
diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
index b1ee902..018c6f1 100644
--- a/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
@@ -1,11 +1,19 @@
 from pathlib import Path
 
 import iris
+import numpy as np
+from cmip7_ancil_constants import REAL_MISSING_DATA_INDICATOR
+from cmip7_PI import CMIP7_PI_YEAR
 
+SOLAR_ARRAY_BEG_YEAR = 1700
+SOLAR_ARRAY_END_YEAR = 2300
+SOLAR_PI_DEFAULT_YEAR_MEAN = 1361.603
 
-def cmip7_solar_dirpath(args, period):
+
+def cmip7_solar_dirpath(args, activity, period):
     return (
         Path(args.cmip7_source_data_dirname)
+        / activity
         / "SOLARIS-HEPPA"
         / args.dataset_version
         / "atmos"
@@ -20,3 +28,49 @@ def load_cmip7_solar_cube(path):
     cubelist = iris.load(path)
     name_constraint = iris.Constraint(name="solar_irradiance")
     return cubelist.extract_cube(name_constraint)
+
+
+def cmip7_solar_year_mean(cube, beg_year, end_year):
+    """
+    Calculate mean TSI values for each year and save them into an array.
+    """
+    NBR_YEARS = SOLAR_ARRAY_END_YEAR - SOLAR_ARRAY_BEG_YEAR + 1
+    solar_array = np.zeros(NBR_YEARS)
+    # Calculate and save the mean annual TSI for each CMIP7 historical year.
+    year_range = range(beg_year, end_year + 1)
+    year_mean = SOLAR_PI_DEFAULT_YEAR_MEAN
+    for year in year_range:
+        year_cons = iris.Constraint(time=lambda cell: cell.point.year == year)
+        year_cube = cube.extract(year_cons)
+        year_mean = year_cube.collapsed("time", iris.analysis.MEAN).data
+        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = year_mean
+        # Save the year mean for the pre-industrial year.
+        if year == CMIP7_PI_YEAR:
+            pi_year_mean = year_mean
+
+    # For the years from SOLAR_ARRAY_BEG_YEAR to beg_year - 1,
+    # set the saved TSI value to the pre-industrial year mean TSI.
+    for year in range(SOLAR_ARRAY_BEG_YEAR, beg_year):
+        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = pi_year_mean
+
+    # For the years from CMIP7_HI_END_YEAR + 1 to SOLAR_ARRAY_END_YEAR,
+    # set the saved TSI value to the real missing data indicator
+    for year in range(end_year + 1, SOLAR_ARRAY_END_YEAR + 1):
+        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = REAL_MISSING_DATA_INDICATOR
+    return solar_array
+
+
+def cmip7_solar_save(args, cube, beg_year, end_year, save_dirpath):
+    """
+    Save the TSI values for each year into a text file.
+    """
+    solar_array = cmip7_solar_year_mean(cube, beg_year, end_year)
+    # Ensure that the save directory exists.
+    save_filepath = save_dirpath / args.save_filename
+    with open(save_filepath, "w") as save_file:
+        for year in range(SOLAR_ARRAY_BEG_YEAR, SOLAR_ARRAY_END_YEAR + 1):
+            year_mean = solar_array[year - SOLAR_ARRAY_BEG_YEAR]
+            if year_mean == REAL_MISSING_DATA_INDICATOR:
+                print(year, f"{year_mean:.1f}", file=save_file)
+            else:
+                print(year, f"{year_mean:.3f}", file=save_file)

From 838b776e69bbe4197ab513658b416609af904569 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 28 May 2026 21:51:30 +1000
Subject: [PATCH 07/50] Ensure that the solar save directory exists

---
 CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
index 018c6f1..3663d7a 100644
--- a/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
@@ -66,6 +66,7 @@ def cmip7_solar_save(args, cube, beg_year, end_year, save_dirpath):
     """
     solar_array = cmip7_solar_year_mean(cube, beg_year, end_year)
     # Ensure that the save directory exists.
+    save_dirpath.mkdir(mode=0o755, parents=True, exist_ok=True)
     save_filepath = save_dirpath / args.save_filename
     with open(save_filepath, "w") as save_file:
         for year in range(SOLAR_ARRAY_BEG_YEAR, SOLAR_ARRAY_END_YEAR + 1):

From c13e8c95f300376de9aa1feaab1539a9b27b90e4 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 28 May 2026 22:02:58 +1000
Subject: [PATCH 08/50] Set pi_year_mean correctly

---
 CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
index 3663d7a..3ba7903 100644
--- a/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_solar.py
@@ -38,7 +38,7 @@ def cmip7_solar_year_mean(cube, beg_year, end_year):
     solar_array = np.zeros(NBR_YEARS)
     # Calculate and save the mean annual TSI for each CMIP7 historical year.
     year_range = range(beg_year, end_year + 1)
-    year_mean = SOLAR_PI_DEFAULT_YEAR_MEAN
+    pi_year_mean = SOLAR_PI_DEFAULT_YEAR_MEAN
     for year in year_range:
         year_cons = iris.Constraint(time=lambda cell: cell.point.year == year)
         year_cube = cube.extract(year_cons)

From 99931c4012f41518c942cbf85973976a05a1c172 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 28 May 2026 22:13:27 +1000
Subject: [PATCH 09/50] Use the correct activity

---
 CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py b/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py
index 2c433fd..d518556 100644
--- a/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py
+++ b/CMIP7/esm1p6/atmosphere/solar/cmip7_PI_solar_generate.py
@@ -48,7 +48,7 @@ def cmip7_pi_solar_patch(solar_irradiance):
 if __name__ == "__main__":
     args = parse_args()
 
-    dirpath = cmip7_solar_dirpath(args, "ScenarioMIP", "fx")
+    dirpath = cmip7_solar_dirpath(args, "CMIP", "fx")
     filename = f"multiple_input4MIPs_solar_CMIP_{args.dataset_version}_gn.nc"
     dataset_path = dirpath / filename
 

From dd1544fcda7d6c800a6175fd86122df53ffadc15 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Fri, 5 Jun 2026 11:30:42 +1000
Subject: [PATCH 10/50] Refactor greenhouse gas forcings and add ScenarioMIP

---
 CMIP7/esm1p6/atmosphere/cmip7_SM.py           |   2 +-
 .../atmosphere/ghg/cmip7_HI_ghg_generate.py   | 187 +-----------------
 .../atmosphere/ghg/cmip7_PI_ghg_generate.py   |   6 +-
 .../atmosphere/ghg/cmip7_SM_ghg_generate.py   |  37 ++++
 CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg.py      |   7 +-
 .../esm1p6/atmosphere/ghg/cmip7_ghg_series.py | 181 +++++++++++++++++
 6 files changed, 234 insertions(+), 186 deletions(-)
 create mode 100644 CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
 create mode 100644 CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_SM.py b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
index 821fd36..44ee485 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_SM.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
@@ -2,7 +2,7 @@
 
 from cmip7_ancil_constants import ANCIL_TODAY
 
-CMIP7_SM_BEG_YEAR = 2025
+CMIP7_SM_BEG_YEAR = 2022
 # Model time interpolation requires an extra year
 CMIP7_SM_END_YEAR = 2100
 CMIP7_SM_NBR_YEARS = CMIP7_SM_END_YEAR + 1 - CMIP7_SM_BEG_YEAR
diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py
index d1fb43f..bbf2ab9 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py
@@ -1,23 +1,14 @@
 from argparse import ArgumentParser
-from collections import OrderedDict
-from pathlib import Path
 
-import cftime
-import f90nml
-import iris
-import numpy as np
 from cmip7_ancil_argparse import dataset_parser, path_parser
 from cmip7_HI import (
     CMIP7_HI_BEG_YEAR,
     CMIP7_HI_END_YEAR,
-    CMIP7_HI_NBR_YEARS,
 )
-from ghg.cmip7_ghg import (
-    GHG_MOLAR_MASS,
-    cmip7_ghg_dirpath,
-    cmip7_ghg_filename,
-    cmip7_ghg_mmr,
-    cmip7_pro_greg_date_constraint_from_years,
+from ghg.cmip7_ghg import GHG_MOLAR_MASS
+from ghg.cmip7_ghg_series import (
+    load_cmip7_ghg_series_mmr,
+    cmip7_ghg_update_namelists_file,
 )
 
 
@@ -33,176 +24,14 @@ def parse_args():
     return parser.parse_args()
 
 
-def load_cmip7_hi_ghg_mmr(args, ghg):
-    cmip7_filepath = cmip7_ghg_dirpath(args, ghg) / cmip7_ghg_filename(
-        args, ghg
-    )
-
-    # Read in the CMIP7 cube.
-    full_cube = iris.load_cube(cmip7_filepath)
-
-    # Check that we have the right greenhouse gas.
-    variable_id = full_cube.metadata.attributes["variable_id"]
-    assert ghg == variable_id
-
-    new_cube = full_cube.copy()
-    # Linearly interpolate to Jan 1 so that UM time interpolation
-    # reproduces annual means
-    new_cube.data[:-1] = 0.5 * (full_cube.data[:-1] + full_cube.data[1:])
-    # Extrapolate the last year
-    new_cube.data[-1] = full_cube.data[-1] + 0.5 * (
-        full_cube.data[-1] - full_cube.data[-2]
-    )
-    cube_time = full_cube.coord("time")
-    units = cube_time.units
-    new_cube_time = new_cube.coord("time")
-    # Can't assign to elements of time.points so create a temporary array
-    tvals = np.array(new_cube_time.points)
-    for i in range(len(tvals)):
-        date = units.num2date(cube_time.points[i])
-        # Interpolated data is Jan 1 of the next year
-        newdate = cftime.DatetimeProlepticGregorian(
-            date.year + 1, 1, 1, 0, 0, 0
-        )
-        tvals[i] = units.date2num(newdate)
-    new_cube_time.points = tvals
-
-    # Extract the historical years.
-    date_constraint = cmip7_pro_greg_date_constraint_from_years(
-        CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR
-    )
-    hi_cube = new_cube.extract(date_constraint)
-
-    # Determine the mass mixing ratio.
-    ghg_mmr_list = []
-    for year in range(CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR + 1):
-        year_constraint = cmip7_pro_greg_date_constraint_from_years(year, year)
-        year_cube = hi_cube.extract(year_constraint)
-        ghg_mmr_list.append(cmip7_ghg_mmr(year_cube, ghg))
-    return ghg_mmr_list
-
-
-def read_namelists_lines_up_to(namelists_filepath, exclude_group):
-    """
-    Read lines from namelists_filepath up to but not including a line that
-    contains the string exclude_group. This function is used to avoid having
-    to use f90nml to reformat an entire namelist file. Versions of f90nml
-    older than v1.5 contain a bug that affects null values in namelists.
-    See https://github.com/marshallward/f90nml/pull/180
-    """
-    if not namelists_filepath.exists():
-        raise FileNotFoundError(
-            f"Namelist file {namelists_filepath} does not exist"
-        )
-    # Read the namelists_filepath file up to but not including
-    # the namelist group given by exclude_group.
-    namelists = []
-    exclude_str = "&" + exclude_group.lower()
-    with open(namelists_filepath) as namelists_file:
-        for line in namelists_file:
-            if exclude_str in line.lower():
-                break
-            namelists.append(line)
-    return "".join(namelists)
-
-
-def format_namelist(namelist, float_format="13.6e"):
-    """
-    Change the namelist formatting to the preferred format.
-    """
-    namelist.float_format = float_format
-    namelist.end_comma = True
-    namelist.false_repr = ".FALSE."
-    namelist.true_repr = ".TRUE."
-    namelist.uppercase = True
-
-
-def cmip7_hi_ghg_namelist_str(ghg_mmr_dict, ghg_namelist_name):
-    """
-    Use the greenhouse gas mass mixing ratios to
-    produce a replacement clmchfcg namelist as a string.
-    """
-    # Map each greenhouse gas to an index in the
-    # historical climate forcing arrays.
-    GHG_HI_NAMELIST_INDEX = {
-        "cfc11": 3,
-        "cfc12": 4,
-        "cfc113": 7,
-        "ch4": 1,
-        "co2": 0,
-        "hcfc22": 8,
-        "hfc125": 9,
-        "hfc134a": 10,
-        "n2o": 2,
-    }
-    GHG_HI_NAMELIST_NBR_SPECIES = 11
-    OLD_REAL_MISSING_DATA_VALUE = -32768.0
-
-    # Create arrays to populate the namelist group.
-    namelist_nyears_shape = (GHG_HI_NAMELIST_NBR_SPECIES,)
-    namelist_nyears = np.full(namelist_nyears_shape, CMIP7_HI_NBR_YEARS)
-    namelist_years_shape = (GHG_HI_NAMELIST_NBR_SPECIES, CMIP7_HI_NBR_YEARS)
-    namelist_years = np.broadcast_to(
-        np.array(range(CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR + 1)),
-        namelist_years_shape,
-    ).T
-    namelist_levls = np.zeros(namelist_years.shape)
-    for ghg in GHG_HI_NAMELIST_INDEX:
-        ghg_index = GHG_HI_NAMELIST_INDEX[ghg]
-        namelist_levls[:, ghg_index] = ghg_mmr_dict[ghg]
-    namelist_rates = np.full(namelist_years.shape, OLD_REAL_MISSING_DATA_VALUE)
-
-    # Create a dictionary to use to patch the namelist group.
-    namelist_dict = {
-        "l_clmchfcg": True,
-        "clim_fcg_nyears": namelist_nyears,
-        "clim_fcg_years": namelist_years,
-        "clim_fcg_levls": namelist_levls,
-        "clim_fcg_rates": namelist_rates,
-    }
-
-    patch = {ghg_namelist_name: OrderedDict(namelist_dict)}
-    patch_namelist = f90nml.namelist.Namelist(patch)
-
-    # Change the namelist arrays to row major.
-    patch_str = str(patch_namelist)
-    parser = f90nml.Parser()
-    parser.row_major = True
-    row_major_patch_namelist = parser.reads(patch_str)
-    # Correctly format the namelist.
-    format_namelist(row_major_patch_namelist)
-    # The format is ignored unless you print the namelist or
-    # convert it to a string.
-    return str(row_major_patch_namelist)
-
-
-def update_namelists_file(ghg_mmr_dict):
-    """
-    Use the greenhouse gas mass mixing ratios in ghg_mmr_dict
-    to replace the greenhouse gas namelist in the relevant namelists file.
-    """
-    namelists_filepath = Path("atmosphere") / "namelists"
-    ghg_namelist_name = "clmchfcg"
-    # Read the original namelists file up to ghg_namelist_name.
-    namelists_str = read_namelists_lines_up_to(
-        namelists_filepath, ghg_namelist_name
-    )
-    # Use ghg_mmr_dict and ghg_namelist_name to create
-    # a replacement namelist as a string.
-    ghg_namelist_str = cmip7_hi_ghg_namelist_str(
-        ghg_mmr_dict, ghg_namelist_name
-    )
-    # Replace the original namelists file.
-    with open(namelists_filepath, "w") as namelists_file:
-        print(namelists_str + ghg_namelist_str, file=namelists_file)
-
-
 if __name__ == "__main__":
     args = parse_args()
 
     ghg_mmr_dict = dict()
     for ghg in GHG_MOLAR_MASS:
-        ghg_mmr_dict[ghg] = load_cmip7_hi_ghg_mmr(args, ghg)
+        ghg_mmr_dict[ghg] = load_cmip7_ghg_series_mmr(
+            args, "CMIP", ghg, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR)
 
     # Patch the greenhouse gas namelist.
-    update_namelists_file(ghg_mmr_dict)
+    cmip7_ghg_update_namelists_file(
+        ghg_mmr_dict, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR)
diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_PI_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_PI_ghg_generate.py
index 358c439..8dd3e6b 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_PI_ghg_generate.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_PI_ghg_generate.py
@@ -28,9 +28,9 @@ def parse_args():
 
 
 def load_cmip7_pi_ghg_mmr(args, ghg):
-    cmip7_filepath = cmip7_ghg_dirpath(args, ghg) / cmip7_ghg_filename(
-        args, ghg
-    )
+    dirpath = cmip7_ghg_dirpath(args, "CMIP", ghg)
+    filename = cmip7_ghg_filename(args, "CMIP", ghg)
+    cmip7_filepath = dirpath / filename
 
     # Read in the CMIP7 cube
     full_cube = iris.load_cube(cmip7_filepath)
diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
new file mode 100644
index 0000000..6186745
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
@@ -0,0 +1,37 @@
+from argparse import ArgumentParser
+
+from cmip7_ancil_argparse import dataset_parser, path_parser
+from cmip7_SM import (
+    CMIP7_SM_BEG_YEAR,
+    CMIP7_SM_END_YEAR,
+)
+from ghg.cmip7_ghg import GHG_MOLAR_MASS
+from ghg.cmip7_ghg_series import (
+    load_cmip7_ghg_series_mmr,
+    cmip7_ghg_update_namelists_file,
+)
+
+
+def parse_args():
+    parser = ArgumentParser(
+        parents=[path_parser(), dataset_parser()],
+        prog="cmip7_SM_ghg_generate",
+        description=(
+            "Generate input files from CMIP7 ScenarioMIP greenhouse gas forcings"
+        ),
+    )
+    parser.add_argument("--dataset-date-range")
+    return parser.parse_args()
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    ghg_mmr_dict = dict()
+    for ghg in GHG_MOLAR_MASS:
+        ghg_mmr_dict[ghg] = load_cmip7_ghg_series_mmr(
+            args, "ScenarioMIP", ghg, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR)
+
+    # Patch the greenhouse gas namelist.
+    cmip7_ghg_update_namelists_file(
+        ghg_mmr_dict, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR)
diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg.py
index 9d9c61d..733901a 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg.py
@@ -35,9 +35,10 @@ def cmip7_ghg_mmr(cube, ghg):
     return conc * ghg_scale * GHG_MOLAR_MASS[ghg] / DRY_AIR_MOLAR_MASS
 
 
-def cmip7_ghg_dirpath(args, ghg):
+def cmip7_ghg_dirpath(args, activity, ghg):
     return (
         Path(args.cmip7_source_data_dirname)
+        / activity
         / "CR"
         / args.dataset_version
         / "atmos"
@@ -48,9 +49,9 @@ def cmip7_ghg_dirpath(args, ghg):
     )
 
 
-def cmip7_ghg_filename(args, ghg):
+def cmip7_ghg_filename(args, activity, ghg):
     return (
-        f"{ghg}_input4MIPs_GHGConcentrations_CMIP_"
+        f"{ghg}_input4MIPs_GHGConcentrations_{activity}_"
         f"{args.dataset_version}_gm_"
         f"{args.dataset_date_range}.nc"
     )
diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py
new file mode 100644
index 0000000..4dbf66d
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py
@@ -0,0 +1,181 @@
+from collections import OrderedDict
+from pathlib import Path
+
+import cftime
+import f90nml
+import iris
+import numpy as np
+
+from ghg.cmip7_ghg import (
+    cmip7_ghg_dirpath,
+    cmip7_ghg_filename,
+    cmip7_ghg_mmr,
+    cmip7_pro_greg_date_constraint_from_years,
+)
+
+
+def load_cmip7_ghg_series_mmr(args, activity, ghg, beg_year, end_year):
+    dirpath = cmip7_ghg_dirpath(args, activity, ghg)
+    filename = cmip7_ghg_filename(args, activity, ghg)
+    cmip7_filepath = dirpath / filename
+
+    # Read in the CMIP7 cube.
+    full_cube = iris.load_cube(cmip7_filepath)
+
+    # Check that we have the right greenhouse gas.
+    variable_id = full_cube.metadata.attributes["variable_id"]
+    assert ghg == variable_id
+
+    new_cube = full_cube.copy()
+    # Linearly interpolate to Jan 1 so that UM time interpolation
+    # reproduces annual means
+    new_cube.data[:-1] = 0.5 * (full_cube.data[:-1] + full_cube.data[1:])
+    # Extrapolate the last year
+    new_cube.data[-1] = full_cube.data[-1] + 0.5 * (
+        full_cube.data[-1] - full_cube.data[-2]
+    )
+    cube_time = full_cube.coord("time")
+    units = cube_time.units
+    new_cube_time = new_cube.coord("time")
+    # Can't assign to elements of time.points so create a temporary array
+    tvals = np.array(new_cube_time.points)
+    for i in range(len(tvals)):
+        date = units.num2date(cube_time.points[i])
+        # Interpolated data is Jan 1 of the next year
+        newdate = cftime.DatetimeProlepticGregorian(
+            date.year + 1, 1, 1, 0, 0, 0
+        )
+        tvals[i] = units.date2num(newdate)
+    new_cube_time.points = tvals
+
+    # Extract the historical years.
+    date_constraint = cmip7_pro_greg_date_constraint_from_years(
+        beg_year, end_year
+    )
+    hi_cube = new_cube.extract(date_constraint)
+
+    # Determine the mass mixing ratio.
+    ghg_mmr_list = []
+    for year in range(beg_year, end_year + 1):
+        year_constraint = cmip7_pro_greg_date_constraint_from_years(year, year)
+        year_cube = hi_cube.extract(year_constraint)
+        ghg_mmr_list.append(cmip7_ghg_mmr(year_cube, ghg))
+    return ghg_mmr_list
+
+
+def read_namelists_lines_up_to(namelists_filepath, exclude_group):
+    """
+    Read lines from namelists_filepath up to but not including a line that
+    contains the string exclude_group. This function is used to avoid having
+    to use f90nml to reformat an entire namelist file. Versions of f90nml
+    older than v1.5 contain a bug that affects null values in namelists.
+    See https://github.com/marshallward/f90nml/pull/180
+    """
+    if not namelists_filepath.exists():
+        raise FileNotFoundError(
+            f"Namelist file {namelists_filepath} does not exist"
+        )
+    # Read the namelists_filepath file up to but not including
+    # the namelist group given by exclude_group.
+    namelists = []
+    exclude_str = "&" + exclude_group.lower()
+    with open(namelists_filepath) as namelists_file:
+        for line in namelists_file:
+            if exclude_str in line.lower():
+                break
+            namelists.append(line)
+    return "".join(namelists)
+
+
+def format_namelist(namelist, float_format="13.6e"):
+    """
+    Change the namelist formatting to the preferred format.
+    """
+    namelist.float_format = float_format
+    namelist.end_comma = True
+    namelist.false_repr = ".FALSE."
+    namelist.true_repr = ".TRUE."
+    namelist.uppercase = True
+
+
+def cmip7_ghg_namelist_str(
+    ghg_mmr_dict, ghg_namelist_name, beg_year, end_year
+):
+    """
+    Use the greenhouse gas mass mixing ratios to
+    produce a replacement clmchfcg namelist as a string.
+    """
+    # Map each greenhouse gas to an index in the
+    # historical climate forcing arrays.
+    GHG_NAMELIST_INDEX = {
+        "cfc11": 3,
+        "cfc12": 4,
+        "cfc113": 7,
+        "ch4": 1,
+        "co2": 0,
+        "hcfc22": 8,
+        "hfc125": 9,
+        "hfc134a": 10,
+        "n2o": 2,
+    }
+    GHG_NAMELIST_NBR_SPECIES = 11
+    OLD_REAL_MISSING_DATA_VALUE = -32768.0
+
+    # Create arrays to populate the namelist group.
+    NBR_YEARS = end_year - beg_year + 1
+    namelist_nyears_shape = (GHG_NAMELIST_NBR_SPECIES,)
+    namelist_nyears = np.full(namelist_nyears_shape, NBR_YEARS)
+    namelist_years_shape = (GHG_NAMELIST_NBR_SPECIES, NBR_YEARS)
+    namelist_years = np.broadcast_to(
+        np.array(range(beg_year, end_year + 1)),
+        namelist_years_shape,
+    ).T
+    namelist_levls = np.zeros(namelist_years.shape)
+    for ghg in GHG_NAMELIST_INDEX:
+        ghg_index = GHG_NAMELIST_INDEX[ghg]
+        namelist_levls[:, ghg_index] = ghg_mmr_dict[ghg]
+    namelist_rates = np.full(namelist_years.shape, OLD_REAL_MISSING_DATA_VALUE)
+
+    # Create a dictionary to use to patch the namelist group.
+    namelist_dict = {
+        "l_clmchfcg": True,
+        "clim_fcg_nyears": namelist_nyears,
+        "clim_fcg_years": namelist_years,
+        "clim_fcg_levls": namelist_levls,
+        "clim_fcg_rates": namelist_rates,
+    }
+
+    patch = {ghg_namelist_name: OrderedDict(namelist_dict)}
+    patch_namelist = f90nml.namelist.Namelist(patch)
+
+    # Change the namelist arrays to row major.
+    patch_str = str(patch_namelist)
+    parser = f90nml.Parser()
+    parser.row_major = True
+    row_major_patch_namelist = parser.reads(patch_str)
+    # Correctly format the namelist.
+    format_namelist(row_major_patch_namelist)
+    # The format is ignored unless you print the namelist or
+    # convert it to a string.
+    return str(row_major_patch_namelist)
+
+
+def cmip7_ghg_update_namelists_file(ghg_mmr_dict, beg_year, end_year):
+    """
+    Use the greenhouse gas mass mixing ratios in ghg_mmr_dict
+    to replace the greenhouse gas namelist in the relevant namelists file.
+    """
+    namelists_filepath = Path("atmosphere") / "namelists"
+    ghg_namelist_name = "clmchfcg"
+    # Read the original namelists file up to ghg_namelist_name.
+    namelists_str = read_namelists_lines_up_to(
+        namelists_filepath, ghg_namelist_name
+    )
+    # Use ghg_mmr_dict and ghg_namelist_name to create
+    # a replacement namelist as a string.
+    ghg_namelist_str = cmip7_ghg_namelist_str(
+        ghg_mmr_dict, ghg_namelist_name, beg_year, end_year
+    )
+    # Replace the original namelists file.
+    with open(namelists_filepath, "w") as namelists_file:
+        print(namelists_str + ghg_namelist_str, file=namelists_file)

From 8585ada70b82f5e12fc6d2ec5bbdba15fb8f28ac Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Fri, 5 Jun 2026 11:41:46 +1000
Subject: [PATCH 11/50] Apply ruff formatting to satisfy pre-commit

---
 CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
index 6186745..f26512d 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
@@ -7,8 +7,8 @@
 )
 from ghg.cmip7_ghg import GHG_MOLAR_MASS
 from ghg.cmip7_ghg_series import (
-    load_cmip7_ghg_series_mmr,
     cmip7_ghg_update_namelists_file,
+    load_cmip7_ghg_series_mmr,
 )
 
 
@@ -17,7 +17,8 @@ def parse_args():
         parents=[path_parser(), dataset_parser()],
         prog="cmip7_SM_ghg_generate",
         description=(
-            "Generate input files from CMIP7 ScenarioMIP greenhouse gas forcings"
+            "Generate input files from CMIP7 ScenarioMIP "
+            "greenhouse gas forcings"
         ),
     )
     parser.add_argument("--dataset-date-range")
@@ -30,8 +31,10 @@ def parse_args():
     ghg_mmr_dict = dict()
     for ghg in GHG_MOLAR_MASS:
         ghg_mmr_dict[ghg] = load_cmip7_ghg_series_mmr(
-            args, "ScenarioMIP", ghg, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR)
+            args, "ScenarioMIP", ghg, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+        )
 
     # Patch the greenhouse gas namelist.
     cmip7_ghg_update_namelists_file(
-        ghg_mmr_dict, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR)
+        ghg_mmr_dict, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+    )

From e5d360fa4a16cb21b641f53a16efed5c58ae5b0e Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Fri, 5 Jun 2026 11:45:49 +1000
Subject: [PATCH 12/50] Apply ruff formatting to satisfy pre-commit

---
 CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py | 8 +++++---
 CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py      | 5 +----
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py
index bbf2ab9..ab8f83b 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py
@@ -7,8 +7,8 @@
 )
 from ghg.cmip7_ghg import GHG_MOLAR_MASS
 from ghg.cmip7_ghg_series import (
-    load_cmip7_ghg_series_mmr,
     cmip7_ghg_update_namelists_file,
+    load_cmip7_ghg_series_mmr,
 )
 
 
@@ -30,8 +30,10 @@ def parse_args():
     ghg_mmr_dict = dict()
     for ghg in GHG_MOLAR_MASS:
         ghg_mmr_dict[ghg] = load_cmip7_ghg_series_mmr(
-            args, "CMIP", ghg, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR)
+            args, "CMIP", ghg, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR
+        )
 
     # Patch the greenhouse gas namelist.
     cmip7_ghg_update_namelists_file(
-        ghg_mmr_dict, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR)
+        ghg_mmr_dict, CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR
+    )
diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py
index 4dbf66d..e3bfe23 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py
@@ -5,7 +5,6 @@
 import f90nml
 import iris
 import numpy as np
-
 from ghg.cmip7_ghg import (
     cmip7_ghg_dirpath,
     cmip7_ghg_filename,
@@ -98,9 +97,7 @@ def format_namelist(namelist, float_format="13.6e"):
     namelist.uppercase = True
 
 
-def cmip7_ghg_namelist_str(
-    ghg_mmr_dict, ghg_namelist_name, beg_year, end_year
-):
+def cmip7_ghg_namelist_str(ghg_mmr_dict, ghg_namelist_name, beg_year, end_year):
     """
     Use the greenhouse gas mass mixing ratios to
     produce a replacement clmchfcg namelist as a string.

From b39bcace52dce4e654ad8cbc2e281ae70d6880aa Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Tue, 9 Jun 2026 17:04:20 +1000
Subject: [PATCH 13/50] Refactor to use dirpath and filename

---
 .../aerosol/cmip7_aerosol_anthro.py           | 22 ++++++++-----------
 1 file changed, 9 insertions(+), 13 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
index d06a438..8d5adef 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
@@ -14,42 +14,38 @@
 )
 
 
-def cmip7_aerosol_anthro_rootpath(args):
+def _anthro_dirpath(args, variable):
     return (
         Path(args.cmip7_source_data_dirname)
+        / "CMIP"
         / "PNNL-JGCRI"
         / args.dataset_version
         / "atmos"
         / "mon"
+        / variable
+        / "gn"
+        / args.dataset_vdate
     )
 
 
 def cmip7_aerosol_air_anthro_filepath(args, species, date_range):
-    rootpath = cmip7_aerosol_anthro_rootpath(args)
+    dirpath = _anthro_dirpath(args, f"{species}_em_AIR_anthro")
     filename = (
         f"{species}-em-AIR-anthro_input4MIPs_emissions_CMIP_"
         f"{args.dataset_version}_gn_"
         f"{date_range}.nc"
     )
-    return (
-        rootpath
-        / f"{species}_em_AIR_anthro"
-        / "gn"
-        / args.dataset_vdate
-        / filename
-    )
+    return dirpath / filename
 
 
 def cmip7_aerosol_anthro_filepath(args, species, date_range):
-    rootpath = cmip7_aerosol_anthro_rootpath(args)
+    dirpath = _anthro_dirpath(args, f"{species}_em_anthro")
     filename = (
         f"{species}-em-anthro_input4MIPs_emissions_CMIP_"
         f"{args.dataset_version}_gn_"
         f"{date_range}.nc"
     )
-    return (
-        rootpath / f"{species}_em_anthro" / "gn" / args.dataset_vdate / filename
-    )
+    return dirpath / filename
 
 
 def load_cmip7_aerosol_anthro(args, species, date_range, constraint):

From 54402033302bb6706afbce5c7033c1d975d11310 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Tue, 9 Jun 2026 17:05:19 +1000
Subject: [PATCH 14/50] Add draft ScenarioMIP scripts in progress

---
 .../aerosol/cmip7_SM_BC_interpolate.py        |  10 ++
 .../aerosol/cmip7_SM_Bio_interpolate.py       |  71 ++++++++++++
 .../aerosol/cmip7_SM_OC_interpolate.py        |  10 ++
 .../aerosol/cmip7_SM_SO2_interpolate.py       |  68 +++++++++++
 .../atmosphere/aerosol/cmip7_SM_aerosol.py    |  25 ++++
 .../aerosol/cmip7_SM_aerosol_anthro.py        | 108 ++++++++++++++++++
 6 files changed, 292 insertions(+)
 create mode 100644 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_BC_interpolate.py
 create mode 100644 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
 create mode 100644 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_OC_interpolate.py
 create mode 100644 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
 create mode 100644 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
 create mode 100644 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_BC_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_BC_interpolate.py
new file mode 100644
index 0000000..6223d9d
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_BC_interpolate.py
@@ -0,0 +1,10 @@
+# Interpolate CMIP7 HI BC emissions to ESM1.6 grid
+from aerosol.cmip7_SM_aerosol_anthro import (
+    cmip7_sm_aerosol_anthro_interpolate,
+    parse_args,
+)
+
+if __name__ == "__main__":
+    args = parse_args(species="BC")
+
+    cmip7_sm_aerosol_anthro_interpolate(args, species="BC", stash_item=129)
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
new file mode 100644
index 0000000..4d88c6e
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -0,0 +1,71 @@
+# Interpolate CMIP7 HI Biomass burning emissions to ESM1.6 grid
+
+from argparse import ArgumentParser
+from ast import literal_eval
+
+from aerosol.cmip7_aerosol_biomass import (
+    load_cmip7_aerosol_biomass,
+    load_cmip7_aerosol_biomass_list,
+    save_cmip7_aerosol_biomass,
+)
+from aerosol.cmip7_SM_aerosol import (
+    CMIP7_SM_AEROSOL_BEG_YEAR,
+    CMIP7_SM_AEROSOL_END_YEAR,
+    esm_sm_aerosol_save_dirpath,
+)
+from cmip7_ancil_argparse import (
+    common_parser,
+    percent_parser,
+)
+from cmip7_ancil_common import cmip7_date_constraint_from_years
+
+
+def parse_args():
+    parser = ArgumentParser(
+        prog="cmip7_SM_Bio_interpolate",
+        description=(
+            "Generate input files from CMIP7 historical biomass forcings"
+        ),
+        parents=[
+            common_parser(),
+            percent_parser(),
+        ],
+    )
+    parser.add_argument("--dataset-date-range-list", type=literal_eval)
+    parser.add_argument("--save-filename")
+    return parser.parse_args()
+
+
+def load_cmip7_sm_aerosol_biomass(args, species):
+    return load_cmip7_aerosol_biomass_list(
+        args,
+        species,
+        args.dataset_date_range_list,
+        cmip7_date_constraint_from_years(
+            CMIP7_SM_AEROSOL_BEG_YEAR,
+            CMIP7_SM_AEROSOL_END_YEAR,
+        ),
+    )
+
+
+def load_cmip7_sm_aerosol_biomass_percentage(args, species):
+    return load_cmip7_aerosol_biomass(
+        args,
+        species,
+        args.percent_date_range,
+        cmip7_date_constraint_from_years(
+            CMIP7_SM_AEROSOL_BEG_YEAR,
+            CMIP7_SM_AEROSOL_END_YEAR,
+        ),
+    )
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    save_cmip7_aerosol_biomass(
+        args,
+        load_cmip7_sm_aerosol_biomass_percentage,
+        load_cmip7_sm_aerosol_biomass,
+        esm_sm_aerosol_save_dirpath(args),
+    )
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_OC_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_OC_interpolate.py
new file mode 100644
index 0000000..7b5c8c5
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_OC_interpolate.py
@@ -0,0 +1,10 @@
+# Interpolate CMIP7 HI OCFF emissions to ESM1.6 grid
+from aerosol.cmip7_SM_aerosol_anthro import (
+    cmip7_sm_aerosol_anthro_interpolate,
+    parse_args,
+)
+
+if __name__ == "__main__":
+    args = parse_args(species="OC")
+
+    cmip7_sm_aerosol_anthro_interpolate(args, species="OC", stash_item=135)
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
new file mode 100644
index 0000000..425bf2c
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -0,0 +1,68 @@
+# Interpolate CMIP7 HI SO2 emissions to ESM1.6 grid
+
+from argparse import ArgumentParser
+from ast import literal_eval
+
+from aerosol.cmip7_SM_aerosol import (
+    CMIP7_SM_AEROSOL_BEG_YEAR,
+    CMIP7_SM_AEROSOL_END_YEAR,
+    esm_sm_aerosol_ancil_dirpath,
+    esm_sm_aerosol_save_dirpath,
+)
+from aerosol.cmip7_SM_aerosol_anthro import load_cmip7_sm_aerosol_anthro
+from aerosol.cmip7_SO2_interpolate import (
+    load_dms,
+    save_cmip7_so2_aerosol_anthro,
+)
+from cmip7_ancil_argparse import (
+    common_parser,
+    dms_filename_parser,
+)
+from cmip7_HI import fix_esm15_sm_ancil_date
+
+
+def parse_args():
+    DMS_ANCIL_FILENAME = "scycl_1849_2015_ESM1_v4.anc"
+
+    parser = ArgumentParser(
+        prog="cmip7_SM_SO2_interpolate",
+        description=("Generate input files from CMIP7 historical SO2 forcings"),
+        parents=[
+            common_parser(),
+            dms_filename_parser(dms_ancil_filename=DMS_ANCIL_FILENAME),
+        ],
+    )
+    parser.add_argument("--dataset-date-range-list", type=literal_eval)
+    parser.add_argument("--save-filename")
+    return parser.parse_args()
+
+
+def load_cmip7_sm_so2_aerosol_anthro(args, species):
+    return load_cmip7_sm_aerosol_anthro(
+        args,
+        species,
+        beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
+        end_year=CMIP7_SM_AEROSOL_END_YEAR,
+    )
+
+
+def load_sm_dms(args):
+    # Use the CMIP6 DMS
+    dms_ancil_dirpath = (
+        esm_sm_aerosol_ancil_dirpath(args.esm15_inputs_dirname)
+        / args.esm_grid_rel_dirname
+        / args.esm15_aerosol_version
+    )
+    return load_dms(args, dms_ancil_dirpath, fix_esm15_sm_ancil_date)
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    save_cmip7_so2_aerosol_anthro(
+        args,
+        load_cmip7_sm_so2_aerosol_anthro,
+        args.dataset_date_range_list,
+        load_sm_dms,
+        esm_sm_aerosol_save_dirpath(args),
+    )
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
new file mode 100644
index 0000000..56b9520
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
@@ -0,0 +1,25 @@
+from pathlib import Path
+
+from cmip7_ancil_constants import ANCIL_TODAY
+from cmip7_HI import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+
+CMIP7_SM_AEROSOL_BEG_YEAR = CMIP7_SM_BEG_YEAR - 1
+CMIP7_SM_AEROSOL_END_YEAR = CMIP7_SM_END_YEAR + 1
+
+
+def esm_sm_aerosol_ancil_dirpath(ancil_root_dirname):
+    return (
+        Path(ancil_root_dirname)
+        / "modern"
+        / "historical"
+        / "atmosphere"
+        / "aerosol"
+    )
+
+
+def esm_sm_aerosol_save_dirpath(args):
+    return (
+        esm_sm_aerosol_ancil_dirpath(args.ancil_target_dirname)
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
new file mode 100644
index 0000000..bc15608
--- /dev/null
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
@@ -0,0 +1,108 @@
+from argparse import ArgumentParser
+from ast import literal_eval
+from pathlib import Path
+
+from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_interpolate
+from aerosol.cmip7_aerosol_common import load_cmip7_aerosol_list
+from aerosol.cmip7_SM_aerosol import (
+    CMIP7_SM_AEROSOL_BEG_YEAR,
+    CMIP7_SM_AEROSOL_END_YEAR,
+    esm_sm_aerosol_save_dirpath,
+)
+from cmip7_ancil_argparse import common_parser
+from cmip7_ancil_common import (
+    cmip7_date_constraint_from_years,
+    fix_coords,
+)
+
+
+def parse_args(species):
+    parser = ArgumentParser(
+        prog=f"cmip7_SM_{species}_interpolate",
+        description=(
+            f"Generate input files from CMIP7 ScenarioMIP {species} forcings"
+        ),
+        parents=[common_parser()],
+    )
+    parser.add_argument("--scenario")
+    parser.add_argument("--dataset-date-range-list", type=literal_eval)
+    parser.add_argument("--save-filename")
+    return parser.parse_args()
+
+
+def _anthro_dirpath(args, variable):
+    return (
+        Path(args.cmip7_source_data_dirname)
+        / "ScenarioMIP"
+        / "IIASA-IAMC"
+        / args.dataset_version
+        / "atmos"
+        / "mon"
+        / variable
+        / "gn"
+        / args.dataset_vdate
+    )
+
+
+def cmip7_sm_aerosol_air_anthro_filepath(args, species, date_range):
+    dirpath = _anthro_dirpath(args, f"{species}_em_AIR_anthro")
+    filename = (
+        f"{species}-em-AIR-anthro_input4MIPs_emissions_ScenarioMIP_"
+        f"{args.dataset_version}_gn_"
+        f"{date_range}.nc"
+    )
+    return dirpath / filename
+
+
+def cmip7_sm_aerosol_anthro_filepath(args, species, date_range):
+    dirpath = _anthro_dirpath(args, f"{species}_em_anthro")
+    filename = (
+        f"{species}-em-anthro_input4MIPs_emissions_ScenarioMIP_"
+        f"{args.dataset_version}_gn_"
+        f"{date_range}.nc"
+    )
+    return dirpath / filename
+
+
+def load_cmip7_sm_aerosol_air_anthro(
+    args,
+    species,
+    beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
+    end_year=CMIP7_SM_AEROSOL_END_YEAR,
+):
+    cube = load_cmip7_aerosol_list(
+        args,
+        cmip7_sm_aerosol_air_anthro_filepath,
+        species,
+        args.dataset_date_range_list,
+        cmip7_date_constraint_from_years(beg_year, end_year),
+    )
+    fix_coords(args, cube)
+    return cube
+
+
+def load_cmip7_sm_aerosol_anthro(
+    args,
+    species,
+    beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
+    end_year=CMIP7_SM_AEROSOL_END_YEAR,
+):
+    cube = load_cmip7_aerosol_list(
+        args,
+        cmip7_sm_aerosol_anthro_filepath,
+        species,
+        args.dataset_date_range_list,
+        cmip7_date_constraint_from_years(beg_year, end_year),
+    )
+    fix_coords(args, cube)
+    return cube
+
+
+def cmip7_sm_aerosol_anthro_interpolate(args, species, stash_item):
+    cmip7_aerosol_anthro_interpolate(
+        args,
+        load_cmip7_sm_aerosol_anthro,
+        species,
+        stash_item,
+        esm_sm_aerosol_save_dirpath(args),
+    )

From 68230ff9e30cb98f30803c9996fb72ef8e6dffa6 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 11 Jun 2026 16:09:44 +1000
Subject: [PATCH 15/50] Use PI DMS for ScenarioMIP SO2 interpolation

---
 .../atmosphere/aerosol/cmip7_HI_SO2_interpolate.py  |  6 +++---
 .../atmosphere/aerosol/cmip7_PI_SO2_interpolate.py  |  6 +++---
 .../atmosphere/aerosol/cmip7_SM_SO2_interpolate.py  | 13 ++++++-------
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py | 12 ++++++------
 4 files changed, 18 insertions(+), 19 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py
index e3aea5d..b4d755b 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py
@@ -20,16 +20,16 @@
 )
 from cmip7_HI import fix_esm15_hi_ancil_date
 
+HI_DMS_ANCIL_FILENAME = "scycl_1849_2015_ESM1_v4.anc"
 
-def parse_args():
-    DMS_ANCIL_FILENAME = "scycl_1849_2015_ESM1_v4.anc"
 
+def parse_args():
     parser = ArgumentParser(
         prog="cmip7_HI_SO2_interpolate",
         description=("Generate input files from CMIP7 historical SO2 forcings"),
         parents=[
             common_parser(),
-            dms_filename_parser(dms_ancil_filename=DMS_ANCIL_FILENAME),
+            dms_filename_parser(dms_ancil_filename=HI_DMS_ANCIL_FILENAME),
         ],
     )
     parser.add_argument("--dataset-date-range-list", type=literal_eval)
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py
index ffe8ee8..8e94943 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py
@@ -14,10 +14,10 @@
 from cmip7_ancil_argparse import common_parser, dms_filename_parser
 from cmip7_PI import fix_esm15_pi_ancil_date
 
+PI_DMS_ANCIL_FILENAME = "scycl_1850_ESM1_v4.anc"
 
-def parse_args():
-    DMS_ANCIL_FILENAME = "scycl_1850_ESM1_v4.anc"
 
+def parse_args():
     parser = ArgumentParser(
         prog="cmip7_PI_SO2_interpolate",
         description=(
@@ -25,7 +25,7 @@ def parse_args():
         ),
         parents=[
             common_parser(),
-            dms_filename_parser(dms_ancil_filename=DMS_ANCIL_FILENAME),
+            dms_filename_parser(dms_ancil_filename=PI_DMS_ANCIL_FILENAME),
         ],
     )
     parser.add_argument("--dataset-date-range")
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
index 425bf2c..b0838d9 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -3,10 +3,11 @@
 from argparse import ArgumentParser
 from ast import literal_eval
 
+from aerosol.cmip7_PI_aerosol import esm_pi_aerosol_ancil_dirpath
+from aerosol.cmip7_PI_SO2_interpolate import PI_DMS_ANCIL_FILENAME
 from aerosol.cmip7_SM_aerosol import (
     CMIP7_SM_AEROSOL_BEG_YEAR,
     CMIP7_SM_AEROSOL_END_YEAR,
-    esm_sm_aerosol_ancil_dirpath,
     esm_sm_aerosol_save_dirpath,
 )
 from aerosol.cmip7_SM_aerosol_anthro import load_cmip7_sm_aerosol_anthro
@@ -18,18 +19,16 @@
     common_parser,
     dms_filename_parser,
 )
-from cmip7_HI import fix_esm15_sm_ancil_date
+from cmip7_PI import fix_esm15_pi_ancil_date
 
 
 def parse_args():
-    DMS_ANCIL_FILENAME = "scycl_1849_2015_ESM1_v4.anc"
-
     parser = ArgumentParser(
         prog="cmip7_SM_SO2_interpolate",
         description=("Generate input files from CMIP7 historical SO2 forcings"),
         parents=[
             common_parser(),
-            dms_filename_parser(dms_ancil_filename=DMS_ANCIL_FILENAME),
+            dms_filename_parser(dms_ancil_filename=PI_DMS_ANCIL_FILENAME),
         ],
     )
     parser.add_argument("--dataset-date-range-list", type=literal_eval)
@@ -49,11 +48,11 @@ def load_cmip7_sm_so2_aerosol_anthro(args, species):
 def load_sm_dms(args):
     # Use the CMIP6 DMS
     dms_ancil_dirpath = (
-        esm_sm_aerosol_ancil_dirpath(args.esm15_inputs_dirname)
+        esm_pi_aerosol_ancil_dirpath(args.esm15_inputs_dirname)
         / args.esm_grid_rel_dirname
         / args.esm15_aerosol_version
     )
-    return load_dms(args, dms_ancil_dirpath, fix_esm15_sm_ancil_date)
+    return load_dms(args, dms_ancil_dirpath, fix_esm15_pi_ancil_date)
 
 
 if __name__ == "__main__":
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
index 56b9520..21620e2 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
@@ -1,17 +1,17 @@
 from pathlib import Path
 
 from cmip7_ancil_constants import ANCIL_TODAY
-from cmip7_HI import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
 
 CMIP7_SM_AEROSOL_BEG_YEAR = CMIP7_SM_BEG_YEAR - 1
 CMIP7_SM_AEROSOL_END_YEAR = CMIP7_SM_END_YEAR + 1
 
 
-def esm_sm_aerosol_ancil_dirpath(ancil_root_dirname):
+def esm_sm_aerosol_ancil_dirpath(args):
     return (
-        Path(ancil_root_dirname)
-        / "modern"
-        / "historical"
+        Path(args.ancil_target_dirname)
+        / "scenarios"
+        / args.scenario
         / "atmosphere"
         / "aerosol"
     )
@@ -19,7 +19,7 @@ def esm_sm_aerosol_ancil_dirpath(ancil_root_dirname):
 
 def esm_sm_aerosol_save_dirpath(args):
     return (
-        esm_sm_aerosol_ancil_dirpath(args.ancil_target_dirname)
+        esm_sm_aerosol_ancil_dirpath(args)
         / args.esm_grid_rel_dirname
         / ANCIL_TODAY
     )

From 41d94aba418ecb693d78c505f6b1d2567ea5a09c Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 11 Jun 2026 17:20:22 +1000
Subject: [PATCH 16/50] Set dirpath correctly

---
 .../esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
index 6077f1a..e30277e 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
@@ -16,24 +16,28 @@
 )
 
 
-def cmip7_aerosol_biomass_rootpath(args):
+def _biomass_dirpath(args, species):
     return (
         Path(args.cmip7_source_data_dirname)
+        / "CMIP"
         / "DRES"
         / args.dataset_version
         / "atmos"
         / "mon"
+        / species
+        / "gn"
+        / args.dataset_vdate
     )
 
 
 def cmip7_aerosol_biomass_filepath(args, species, date_range):
-    rootpath = cmip7_aerosol_biomass_rootpath(args)
+    dirpath = _biomass_dirpath(args, species)
     filename = (
         f"{species}_input4MIPs_emissions_CMIP_"
         f"{args.dataset_version}_gn_"
         f"{date_range}.nc"
     )
-    return rootpath / species / "gn" / args.dataset_vdate / filename
+    return dirpath / filename
 
 
 def load_cmip7_aerosol_biomass(args, species, date_range, constraint):

From b165688fe9709842afbcbc1a674140842bc9866c Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 11 Jun 2026 21:15:42 +1000
Subject: [PATCH 17/50] Use --dataset-date-range

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py | 3 +--
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py | 3 +--
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py  | 3 +--
 3 files changed, 3 insertions(+), 6 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index 4d88c6e..feb0773 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -1,7 +1,6 @@
 # Interpolate CMIP7 HI Biomass burning emissions to ESM1.6 grid
 
 from argparse import ArgumentParser
-from ast import literal_eval
 
 from aerosol.cmip7_aerosol_biomass import (
     load_cmip7_aerosol_biomass,
@@ -31,7 +30,7 @@ def parse_args():
             percent_parser(),
         ],
     )
-    parser.add_argument("--dataset-date-range-list", type=literal_eval)
+    parser.add_argument("--dataset-date-range")
     parser.add_argument("--save-filename")
     return parser.parse_args()
 
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
index b0838d9..8ff0d0a 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -1,7 +1,6 @@
 # Interpolate CMIP7 HI SO2 emissions to ESM1.6 grid
 
 from argparse import ArgumentParser
-from ast import literal_eval
 
 from aerosol.cmip7_PI_aerosol import esm_pi_aerosol_ancil_dirpath
 from aerosol.cmip7_PI_SO2_interpolate import PI_DMS_ANCIL_FILENAME
@@ -31,7 +30,7 @@ def parse_args():
             dms_filename_parser(dms_ancil_filename=PI_DMS_ANCIL_FILENAME),
         ],
     )
-    parser.add_argument("--dataset-date-range-list", type=literal_eval)
+    parser.add_argument("--dataset-date-range")
     parser.add_argument("--save-filename")
     return parser.parse_args()
 
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
index bc15608..1b818e9 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
@@ -1,5 +1,4 @@
 from argparse import ArgumentParser
-from ast import literal_eval
 from pathlib import Path
 
 from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_interpolate
@@ -25,7 +24,7 @@ def parse_args(species):
         parents=[common_parser()],
     )
     parser.add_argument("--scenario")
-    parser.add_argument("--dataset-date-range-list", type=literal_eval)
+    parser.add_argument("--dataset-date-range")
     parser.add_argument("--save-filename")
     return parser.parse_args()
 

From 2386c92965edf802c7bab8b9f5d477597aa9d938 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 11 Jun 2026 21:40:44 +1000
Subject: [PATCH 18/50] Consistently use dataset_date_range

---
 .../atmosphere/aerosol/cmip7_SM_Bio_interpolate.py     |  5 ++---
 .../atmosphere/aerosol/cmip7_SM_SO2_interpolate.py     |  2 +-
 .../atmosphere/aerosol/cmip7_SM_aerosol_anthro.py      | 10 +++++-----
 3 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index feb0773..364f9a3 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -4,7 +4,6 @@
 
 from aerosol.cmip7_aerosol_biomass import (
     load_cmip7_aerosol_biomass,
-    load_cmip7_aerosol_biomass_list,
     save_cmip7_aerosol_biomass,
 )
 from aerosol.cmip7_SM_aerosol import (
@@ -36,10 +35,10 @@ def parse_args():
 
 
 def load_cmip7_sm_aerosol_biomass(args, species):
-    return load_cmip7_aerosol_biomass_list(
+    return load_cmip7_aerosol_biomass(
         args,
         species,
-        args.dataset_date_range_list,
+        args.dataset_date_range,
         cmip7_date_constraint_from_years(
             CMIP7_SM_AEROSOL_BEG_YEAR,
             CMIP7_SM_AEROSOL_END_YEAR,
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
index 8ff0d0a..f5896fe 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -60,7 +60,7 @@ def load_sm_dms(args):
     save_cmip7_so2_aerosol_anthro(
         args,
         load_cmip7_sm_so2_aerosol_anthro,
-        args.dataset_date_range_list,
+        args.dataset_date_range,
         load_sm_dms,
         esm_sm_aerosol_save_dirpath(args),
     )
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
index 1b818e9..e3fd2d9 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
@@ -2,7 +2,7 @@
 from pathlib import Path
 
 from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_interpolate
-from aerosol.cmip7_aerosol_common import load_cmip7_aerosol_list
+from aerosol.cmip7_aerosol_common import load_cmip7_aerosol
 from aerosol.cmip7_SM_aerosol import (
     CMIP7_SM_AEROSOL_BEG_YEAR,
     CMIP7_SM_AEROSOL_END_YEAR,
@@ -69,11 +69,11 @@ def load_cmip7_sm_aerosol_air_anthro(
     beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
     end_year=CMIP7_SM_AEROSOL_END_YEAR,
 ):
-    cube = load_cmip7_aerosol_list(
+    cube = load_cmip7_aerosol(
         args,
         cmip7_sm_aerosol_air_anthro_filepath,
         species,
-        args.dataset_date_range_list,
+        args.dataset_date_range,
         cmip7_date_constraint_from_years(beg_year, end_year),
     )
     fix_coords(args, cube)
@@ -86,11 +86,11 @@ def load_cmip7_sm_aerosol_anthro(
     beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
     end_year=CMIP7_SM_AEROSOL_END_YEAR,
 ):
-    cube = load_cmip7_aerosol_list(
+    cube = load_cmip7_aerosol(
         args,
         cmip7_sm_aerosol_anthro_filepath,
         species,
-        args.dataset_date_range_list,
+        args.dataset_date_range,
         cmip7_date_constraint_from_years(beg_year, end_year),
     )
     fix_coords(args, cube)

From e9607d985656e7a2a236018119691243e2ef078a Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Thu, 11 Jun 2026 21:54:38 +1000
Subject: [PATCH 19/50] Add --scenario

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py | 1 +
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py | 1 +
 2 files changed, 2 insertions(+)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index 364f9a3..159a376 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -29,6 +29,7 @@ def parse_args():
             percent_parser(),
         ],
     )
+    parser.add_argument("--scenario")
     parser.add_argument("--dataset-date-range")
     parser.add_argument("--save-filename")
     return parser.parse_args()
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
index f5896fe..073ddc9 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -30,6 +30,7 @@ def parse_args():
             dms_filename_parser(dms_ancil_filename=PI_DMS_ANCIL_FILENAME),
         ],
     )
+    parser.add_argument("--scenario")
     parser.add_argument("--dataset-date-range")
     parser.add_argument("--save-filename")
     return parser.parse_args()

From e423db3050cf73007d3c71e34c95185031df7997 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Fri, 12 Jun 2026 11:41:56 +1000
Subject: [PATCH 20/50] In load_sector_dict, use the correct filepath

---
 .../atmosphere/aerosol/cmip7_HI_SO2_interpolate.py     |  2 ++
 .../atmosphere/aerosol/cmip7_PI_SO2_interpolate.py     |  2 ++
 .../atmosphere/aerosol/cmip7_SM_SO2_interpolate.py     | 10 ++++++++--
 .../esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py |  9 ++++-----
 4 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py
index b4d755b..84a9f9d 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py
@@ -3,6 +3,7 @@
 from argparse import ArgumentParser
 from ast import literal_eval
 
+from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_filepath
 from aerosol.cmip7_HI_aerosol import (
     CMIP7_HI_AEROSOL_BEG_YEAR,
     CMIP7_HI_AEROSOL_END_YEAR,
@@ -62,6 +63,7 @@ def load_hi_dms(args):
     save_cmip7_so2_aerosol_anthro(
         args,
         load_cmip7_hi_so2_aerosol_anthro,
+        cmip7_aerosol_anthro_filepath,
         args.dataset_date_range_list,
         load_hi_dms,
         esm_hi_aerosol_save_dirpath(args),
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py
index 8e94943..7e1cac8 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py
@@ -2,6 +2,7 @@
 
 from argparse import ArgumentParser
 
+from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_filepath
 from aerosol.cmip7_PI_aerosol import (
     esm_pi_aerosol_ancil_dirpath,
     esm_pi_aerosol_save_dirpath,
@@ -49,6 +50,7 @@ def load_pi_dms(args):
     save_cmip7_so2_aerosol_anthro(
         args,
         load_cmip7_pi_aerosol_anthro,
+        cmip7_aerosol_anthro_filepath,
         args.dataset_date_range,
         load_pi_dms,
         esm_pi_aerosol_save_dirpath(args),
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
index 073ddc9..9421ceb 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -9,7 +9,10 @@
     CMIP7_SM_AEROSOL_END_YEAR,
     esm_sm_aerosol_save_dirpath,
 )
-from aerosol.cmip7_SM_aerosol_anthro import load_cmip7_sm_aerosol_anthro
+from aerosol.cmip7_SM_aerosol_anthro import (
+    cmip7_sm_aerosol_anthro_filepath,
+    load_cmip7_sm_aerosol_anthro,
+)
 from aerosol.cmip7_SO2_interpolate import (
     load_dms,
     save_cmip7_so2_aerosol_anthro,
@@ -24,7 +27,9 @@
 def parse_args():
     parser = ArgumentParser(
         prog="cmip7_SM_SO2_interpolate",
-        description=("Generate input files from CMIP7 historical SO2 forcings"),
+        description=(
+            "Generate input files from CMIP7 ScenarioMIP SO2 forcings"
+        ),
         parents=[
             common_parser(),
             dms_filename_parser(dms_ancil_filename=PI_DMS_ANCIL_FILENAME),
@@ -61,6 +66,7 @@ def load_sm_dms(args):
     save_cmip7_so2_aerosol_anthro(
         args,
         load_cmip7_sm_so2_aerosol_anthro,
+        cmip7_sm_aerosol_anthro_filepath,
         args.dataset_date_range,
         load_sm_dms,
         esm_sm_aerosol_save_dirpath(args),
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
index 993ab20..5629d5f 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
@@ -7,7 +7,6 @@
 import iris
 import netCDF4
 import numpy as np
-from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_filepath
 from aerosol.cmip7_aerosol_common import zero_poles
 from cmip7_ancil_common import (
     INTERPOLATION_SCHEME,
@@ -21,14 +20,14 @@
 )
 
 
-def load_sector_dict(args, date_range_maybe_list):
+def load_sector_dict(args, filepath_fn, date_range_maybe_list):
     date_range = (
         date_range_maybe_list[0]
         if isinstance(date_range_maybe_list, list)
         else date_range_maybe_list
     )
     # Iris doesn't read the sector coordinate so use netCDF4
-    d = netCDF4.Dataset(cmip7_aerosol_anthro_filepath(args, "SO2", date_range))
+    d = netCDF4.Dataset(filepath_fn(args, "SO2", date_range))
     sectord = dict()
     for s in d["sector"].ids.split(";"):
         i, name = s.split(":")
@@ -59,12 +58,12 @@ def tile_yearly_data(from_cube, to_cube):
 
 
 def save_cmip7_so2_aerosol_anthro(
-    args, cmip7_load_fn, date_range, dms_load_fn, save_dirpath
+    args, cmip7_load_fn, filepath_fn, date_range, dms_load_fn, save_dirpath
 ):
     cmip7_so2 = cmip7_load_fn(args, "SO2")
 
     # Iris doesn't read the sector coordinate
-    sectord = load_sector_dict(args, date_range)
+    sectord = load_sector_dict(args, filepath_fn, date_range)
 
     cmip7_so2_high = (
         cmip7_so2[:, sectord["Energy"]]

From 47e33c94b3a058ff4cbefe87938086439675785f Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Fri, 12 Jun 2026 12:25:41 +1000
Subject: [PATCH 21/50] Add debug prints

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
index 5629d5f..817380c 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
@@ -119,5 +119,10 @@ def save_cmip7_so2_aerosol_anthro(
             (dms_lon, 2),
         ],
     )
-
+    print("so2_low")
+    print(so2_low)
+    print("so2_high")
+    print(so2_high)
+    print("dms")
+    print(dms)
     save_ancil([so2_low, so2_high, dms], save_dirpath, args.save_filename)

From 78e425c446dcff9d79eb106635f0b070bb0fe5d4 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Fri, 12 Jun 2026 14:00:41 +1000
Subject: [PATCH 22/50] When calling save_ancil, set replace_bounds=True

---
 .../atmosphere/aerosol/cmip7_SO2_interpolate.py     | 13 ++++++-------
 .../atmosphere/aerosol/cmip7_aerosol_anthro.py      |  2 +-
 .../atmosphere/aerosol/cmip7_aerosol_biomass.py     |  7 ++++++-
 3 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
index 817380c..24f4ff0 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
@@ -119,10 +119,9 @@ def save_cmip7_so2_aerosol_anthro(
             (dms_lon, 2),
         ],
     )
-    print("so2_low")
-    print(so2_low)
-    print("so2_high")
-    print(so2_high)
-    print("dms")
-    print(dms)
-    save_ancil([so2_low, so2_high, dms], save_dirpath, args.save_filename)
+    save_ancil(
+        [so2_low, so2_high, dms],
+        save_dirpath,
+        args.save_filename,
+        replace_bounds=True,
+    )
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
index 8d5adef..5be2a5a 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
@@ -93,4 +93,4 @@ def cmip7_aerosol_anthro_interpolate(
     esm_cube.attributes["STASH"] = iris.fileformats.pp.STASH(
         model=1, section=0, item=stash_item
     )
-    save_ancil(esm_cube, save_dirpath, args.save_filename)
+    save_ancil(esm_cube, save_dirpath, args.save_filename, replace_bounds=True)
diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
index e30277e..3d66643 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
@@ -137,7 +137,12 @@ def save_cmip7_aerosol_biomass(args, load_pc_fn, load_fn, save_dirpath):
         model=1, section=0, item=131
     )
 
-    save_ancil([low_esm, high_esm], save_dirpath, args.save_filename)
+    save_ancil(
+        [low_esm, high_esm],
+        save_dirpath,
+        args.save_filename,
+        replace_bounds=True,
+    )
 
     now = datetime.now()
     print(f"{now}: save_ancil done")

From 7de7fafbdf14e4ebabcadc88cc0946f0cb0edef5 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 12:10:44 +1000
Subject: [PATCH 23/50] Add interpolate_monthly (AI assisted)

---
 .../aerosol/cmip7_SM_aerosol_anthro.py        | 17 ++++-
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 76 +++++++++++++++++++
 2 files changed, 89 insertions(+), 4 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
index e3fd2d9..976231b 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
@@ -8,10 +8,13 @@
     CMIP7_SM_AEROSOL_END_YEAR,
     esm_sm_aerosol_save_dirpath,
 )
+from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
 from cmip7_ancil_argparse import common_parser
 from cmip7_ancil_common import (
     cmip7_date_constraint_from_years,
+    extend_years,
     fix_coords,
+    interpolate_monthly,
 )
 
 
@@ -74,10 +77,13 @@ def load_cmip7_sm_aerosol_air_anthro(
         cmip7_sm_aerosol_air_anthro_filepath,
         species,
         args.dataset_date_range,
-        cmip7_date_constraint_from_years(beg_year, end_year),
+        cmip7_date_constraint_from_years(CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR),
     )
     fix_coords(args, cube)
-    return cube
+    interpolated = interpolate_monthly(
+        cube, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+    )
+    return extend_years(interpolated)
 
 
 def load_cmip7_sm_aerosol_anthro(
@@ -91,10 +97,13 @@ def load_cmip7_sm_aerosol_anthro(
         cmip7_sm_aerosol_anthro_filepath,
         species,
         args.dataset_date_range,
-        cmip7_date_constraint_from_years(beg_year, end_year),
+        cmip7_date_constraint_from_years(CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR),
     )
     fix_coords(args, cube)
-    return cube
+    interpolated = interpolate_monthly(
+        cube, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+    )
+    return extend_years(interpolated)
 
 
 def cmip7_sm_aerosol_anthro_interpolate(args, species, stash_item):
diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index f76bae1..8e1d07f 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -8,6 +8,7 @@
 import cf_units
 import cftime
 import iris
+import iris.analysis
 import mule
 import numpy as np
 from cmip7_ancil_constants import (
@@ -135,6 +136,81 @@ def extend_years(cube):
     return cubelist.concatenate_cube()
 
 
+def interpolate_monthly(cube, beg_year, end_year):
+    # Get original time units and calendar
+    time_coord = cube.coord("time")
+    units = time_coord.units
+    calendar = units.calendar if units.calendar else "standard"
+
+    # Generate target monthly midpoints and bounds
+    tdates = []
+    tbounds = []
+    for year in range(beg_year, end_year + 1):
+        for month in range(1, MONTHS_IN_A_YEAR + 1):
+            beg_date = cftime.datetime(year, month, 1, calendar=calendar)
+            if month == MONTHS_IN_A_YEAR:
+                end_date = cftime.datetime(year + 1, 1, 1, calendar=calendar)
+            else:
+                end_date = cftime.datetime(year, month + 1, 1, calendar=calendar)
+
+            mid_date = beg_date + (end_date - beg_date) / 2
+            tdates.append(mid_date)
+            tbounds.append([beg_date, end_date])
+
+    tpoints = [units.date2num(d) for d in tdates]
+    tbounds_num = [
+        [units.date2num(b[0]), units.date2num(b[1])] for b in tbounds
+    ]
+
+    # Perform linear time-interpolation
+    new_cube = cube.interpolate(
+        [("time", tpoints)], iris.analysis.Linear()
+    )
+
+    # Create new DimCoord with contiguous bounds
+    new_time_coord = iris.coords.DimCoord(
+        tpoints,
+        standard_name="time",
+        bounds=tbounds_num,
+        units=units,
+    )
+    time_dim = cube.coord_dims("time")
+    new_cube.remove_coord("time")
+    new_cube.add_dim_coord(new_time_coord, time_dim)
+
+    # Overwrite interpolated data with original data for months that exist,
+    # matching by the bounds interval (start month and end month).
+    # This works only for cubes that already have time bounds.
+    if time_coord.has_bounds():
+        bounds = time_coord.bounds
+        # Map target bounds intervals:
+        # (start_year, start_month, end_year, end_month) -> index.
+        tbounds_dict = {
+            (b[0].year, b[0].month, b[1].year, b[1].month): index
+            for index, b in enumerate(tbounds)
+        }
+        # Ensure data is a writable array (not read-only memmap).
+        new_cube.data = np.array(new_cube.data)
+        for i in range(len(time_coord.points)):
+            # Convert original bounds of slice i to date objects.
+            beg_date = units.num2date(bounds[i][0])
+            end_date = units.num2date(bounds[i][1])
+            # Match bounds based on start and end year/month components
+            key = (
+                beg_date.year, beg_date.month, end_date.year, end_date.month
+            )
+            print(
+                "In interpolate_monthly(): "
+                f"matched {beg_date.year}/{beg_date.month}, "
+                f"{end_date.year}/{end_date.month}"
+            )
+            if key in tbounds_dict:
+                target_idx = tbounds_dict[key]
+                new_cube.data[target_idx] = cube.data[i]
+
+    return new_cube
+
+
 def set_coord_system(cube):
     coord_system = iris.coord_systems.GeogCS(6371229.0)
     cube.coord("latitude").coord_system = coord_system

From 7e0867ded90e8ffcb85c5b98e959b4affa0e289f Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 13:09:08 +1000
Subject: [PATCH 24/50] In call to DimCoord, copy more parameters from
 time_coord

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 8e1d07f..fa4783c 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -170,9 +170,12 @@ def interpolate_monthly(cube, beg_year, end_year):
     # Create new DimCoord with contiguous bounds
     new_time_coord = iris.coords.DimCoord(
         tpoints,
-        standard_name="time",
-        bounds=tbounds_num,
+        standard_name=time_coord.standard_name,
+        long_name=time_coord.long_name,
+        var_name=time_coord.var_name,
         units=units,
+        bounds=tbounds_num,
+        attributes=time_coord.attributes,
     )
     time_dim = cube.coord_dims("time")
     new_cube.remove_coord("time")

From 5ad101ab1669c47766eb8652c910747ac439d5bc Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 13:22:58 +1000
Subject: [PATCH 25/50] In extend_years, print time coordinates

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index fa4783c..3c82baf 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -131,6 +131,10 @@ def extend_years(cube):
     if time_coord.has_bounds():
         end_year_tc.bounds = end_year_tc.bounds + length_one_year
 
+    # Print time coordinates.
+    print(beg_year_tc)
+    print(time_coord)
+    print(end_year_tc)
     # Return a cube with extended years.
     cubelist = iris.cube.CubeList((beg_year, cube, end_year))
     return cubelist.concatenate_cube()
@@ -202,11 +206,6 @@ def interpolate_monthly(cube, beg_year, end_year):
             key = (
                 beg_date.year, beg_date.month, end_date.year, end_date.month
             )
-            print(
-                "In interpolate_monthly(): "
-                f"matched {beg_date.year}/{beg_date.month}, "
-                f"{end_date.year}/{end_date.month}"
-            )
             if key in tbounds_dict:
                 target_idx = tbounds_dict[key]
                 new_cube.data[target_idx] = cube.data[i]

From 3aa81f364a681a3f3749758113b55deb3af5e762 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 13:45:26 +1000
Subject: [PATCH 26/50] Add a metatdata comparison (AI assisted)

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 3c82baf..98141cc 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -135,6 +135,16 @@ def extend_years(cube):
     print(beg_year_tc)
     print(time_coord)
     print(end_year_tc)
+    # Compare metadata.
+    print("beg_year_tc.metadata == time_coord.metadata:",
+          beg_year_tc.metadata == time_coord.metadata
+    )
+    if beg_year_tc.metadata != time_coord.metadata:
+        for field in beg_year_tc.metadata._fields:
+            v1 = getattr(beg_year_tc.metadata, field)
+            v2 = getattr(time_coord.metadata, field)
+            if v1 != v2:
+                print(f"Diff in {field}: beg={v1!r}, mid={v2!r}")
     # Return a cube with extended years.
     cubelist = iris.cube.CubeList((beg_year, cube, end_year))
     return cubelist.concatenate_cube()

From e43a6fafcc632d195f6d6b826d30065607d188de Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 15:30:50 +1000
Subject: [PATCH 27/50] Add more diagnostics

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 21 +++++++++++--------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 98141cc..1997d4d 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -136,15 +136,18 @@ def extend_years(cube):
     print(time_coord)
     print(end_year_tc)
     # Compare metadata.
-    print("beg_year_tc.metadata == time_coord.metadata:",
-          beg_year_tc.metadata == time_coord.metadata
-    )
-    if beg_year_tc.metadata != time_coord.metadata:
-        for field in beg_year_tc.metadata._fields:
-            v1 = getattr(beg_year_tc.metadata, field)
-            v2 = getattr(time_coord.metadata, field)
-            if v1 != v2:
-                print(f"Diff in {field}: beg={v1!r}, mid={v2!r}")
+    for tc in [beg_year_tc, end_year_tc]:
+        print("tc.metadata == time_coord.metadata:",
+              tc.metadata == time_coord.metadata
+              )
+        if tc.metadata != time_coord.metadata:
+            for field in tc.metadata._fields:
+                v1 = getattr(tc.metadata, field)
+                v2 = getattr(time_coord.metadata, field)
+                if v1 != v2:
+                    print(f"Diff in {field}: beg={v1!r}, mid={v2!r}")
+    iris.util.describe_diff(beg_year, cube)
+    iris.util.describe_diff(end_year, cube)
     # Return a cube with extended years.
     cubelist = iris.cube.CubeList((beg_year, cube, end_year))
     return cubelist.concatenate_cube()

From 55424b47460417568df2fcf3bd6fc84f40740cb8 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 17:02:20 +1000
Subject: [PATCH 28/50] Addm more diagnostics

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 1997d4d..78633a0 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -135,6 +135,15 @@ def extend_years(cube):
     print(beg_year_tc)
     print(time_coord)
     print(end_year_tc)
+    print("--- Dtype Diagnostics ---")
+    print("beg_year_tc points dtype:", beg_year_tc.core_points().dtype)
+    print("time_coord points dtype:", time_coord.core_points().dtype)
+    print("end_year_tc points dtype:", end_year_tc.core_points().dtype)
+    if time_coord.has_bounds():
+        print("beg_year_tc bounds dtype:", beg_year_tc.core_bounds().dtype)
+        print("time_coord bounds dtype:", time_coord.core_bounds().dtype)
+        print("end_year_tc bounds dtype:", end_year_tc.core_bounds().dtype)
+    print("-------------------------")
     # Compare metadata.
     for tc in [beg_year_tc, end_year_tc]:
         print("tc.metadata == time_coord.metadata:",

From 157147bed640b31aa43b1d22aad84ac4331d07fe Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 17:18:35 +1000
Subject: [PATCH 29/50] Use float64 for the interpolated cube in
 interpolate_monthly

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 78633a0..94a9801 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -183,10 +183,11 @@ def interpolate_monthly(cube, beg_year, end_year):
             tdates.append(mid_date)
             tbounds.append([beg_date, end_date])
 
-    tpoints = [units.date2num(d) for d in tdates]
-    tbounds_num = [
-        [units.date2num(b[0]), units.date2num(b[1])] for b in tbounds
-    ]
+    tpoints = np.array([units.date2num(d) for d in tdates], dtype=np.float64)
+    tbounds_num = np.array(
+        [[units.date2num(b[0]), units.date2num(b[1])] for b in tbounds],
+        dtype=np.float64,
+    )
 
     # Perform linear time-interpolation
     new_cube = cube.interpolate(

From 1e8320f5f4f45fd65c9b7fc93ebcf31e9ad3c02f Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 17:39:29 +1000
Subject: [PATCH 30/50] Use np.ma.asarray

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 94a9801..2569ad0 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -220,7 +220,7 @@ def interpolate_monthly(cube, beg_year, end_year):
             for index, b in enumerate(tbounds)
         }
         # Ensure data is a writable array (not read-only memmap).
-        new_cube.data = np.array(new_cube.data)
+        new_cube.data = np.ma.asarray(new_cube.data)
         for i in range(len(time_coord.points)):
             # Convert original bounds of slice i to date objects.
             beg_date = units.num2date(bounds[i][0])
@@ -233,6 +233,7 @@ def interpolate_monthly(cube, beg_year, end_year):
                 target_idx = tbounds_dict[key]
                 new_cube.data[target_idx] = cube.data[i]
 
+    new_cube.data = np.ma.asarray(new_cube.data)
     return new_cube
 
 

From 49d38503817ac24f7e56b8e03ae3e38ee78b542b Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 20:45:35 +1000
Subject: [PATCH 31/50] Try reverting replace_bounds back to False

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
index 3d66643..e30277e 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_biomass.py
@@ -137,12 +137,7 @@ def save_cmip7_aerosol_biomass(args, load_pc_fn, load_fn, save_dirpath):
         model=1, section=0, item=131
     )
 
-    save_ancil(
-        [low_esm, high_esm],
-        save_dirpath,
-        args.save_filename,
-        replace_bounds=True,
-    )
+    save_ancil([low_esm, high_esm], save_dirpath, args.save_filename)
 
     now = datetime.now()
     print(f"{now}: save_ancil done")

From 52ed74a27c78ed038a165ff09d886e92f1a841e1 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 21:03:15 +1000
Subject: [PATCH 32/50] Try reverting replace_bounds back to False (2)

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
index 24f4ff0..044b734 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py
@@ -123,5 +123,4 @@ def save_cmip7_so2_aerosol_anthro(
         [so2_low, so2_high, dms],
         save_dirpath,
         args.save_filename,
-        replace_bounds=True,
     )

From 871dbb0ac0c1b6540776e368a57ffe0652ec9c3a Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 21:09:38 +1000
Subject: [PATCH 33/50] Try reverting replace_bounds back to False (3)

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
index 5be2a5a..8d5adef 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_anthro.py
@@ -93,4 +93,4 @@ def cmip7_aerosol_anthro_interpolate(
     esm_cube.attributes["STASH"] = iris.fileformats.pp.STASH(
         model=1, section=0, item=stash_item
     )
-    save_ancil(esm_cube, save_dirpath, args.save_filename, replace_bounds=True)
+    save_ancil(esm_cube, save_dirpath, args.save_filename)

From 1905719aa6dd4eefedeaefdebdb833705fb179ec Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Sat, 13 Jun 2026 21:19:54 +1000
Subject: [PATCH 34/50] Revert the diagnostics in extend_years

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 26 -------------------
 1 file changed, 26 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 2569ad0..3e1e7aa 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -131,32 +131,6 @@ def extend_years(cube):
     if time_coord.has_bounds():
         end_year_tc.bounds = end_year_tc.bounds + length_one_year
 
-    # Print time coordinates.
-    print(beg_year_tc)
-    print(time_coord)
-    print(end_year_tc)
-    print("--- Dtype Diagnostics ---")
-    print("beg_year_tc points dtype:", beg_year_tc.core_points().dtype)
-    print("time_coord points dtype:", time_coord.core_points().dtype)
-    print("end_year_tc points dtype:", end_year_tc.core_points().dtype)
-    if time_coord.has_bounds():
-        print("beg_year_tc bounds dtype:", beg_year_tc.core_bounds().dtype)
-        print("time_coord bounds dtype:", time_coord.core_bounds().dtype)
-        print("end_year_tc bounds dtype:", end_year_tc.core_bounds().dtype)
-    print("-------------------------")
-    # Compare metadata.
-    for tc in [beg_year_tc, end_year_tc]:
-        print("tc.metadata == time_coord.metadata:",
-              tc.metadata == time_coord.metadata
-              )
-        if tc.metadata != time_coord.metadata:
-            for field in tc.metadata._fields:
-                v1 = getattr(tc.metadata, field)
-                v2 = getattr(time_coord.metadata, field)
-                if v1 != v2:
-                    print(f"Diff in {field}: beg={v1!r}, mid={v2!r}")
-    iris.util.describe_diff(beg_year, cube)
-    iris.util.describe_diff(end_year, cube)
     # Return a cube with extended years.
     cubelist = iris.cube.CubeList((beg_year, cube, end_year))
     return cubelist.concatenate_cube()

From d7abf90f35fddbfa1dc45aa3f4d48fdfcb08e5ed Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 11:06:10 +1000
Subject: [PATCH 35/50] Sum IIASA-IAMC sectors (partly AI generated)

---
 .../aerosol/cmip7_SM_Bio_interpolate.py       | 194 +++++++++++++++---
 1 file changed, 161 insertions(+), 33 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index 159a376..9acef9d 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -1,33 +1,35 @@
-# Interpolate CMIP7 HI Biomass burning emissions to ESM1.6 grid
+# Interpolate CMIP7 SM Biomass burning emissions to ESM1.6 grid
 
 from argparse import ArgumentParser
+from datetime import datetime
+from pathlib import Path
 
-from aerosol.cmip7_aerosol_biomass import (
-    load_cmip7_aerosol_biomass,
-    save_cmip7_aerosol_biomass,
+import iris
+import netCDF4
+from aerosol.cmip7_aerosol_common import (
+    load_cmip7_aerosol,
+    zero_poles,
 )
-from aerosol.cmip7_SM_aerosol import (
-    CMIP7_SM_AEROSOL_BEG_YEAR,
-    CMIP7_SM_AEROSOL_END_YEAR,
-    esm_sm_aerosol_save_dirpath,
+from aerosol.cmip7_SM_aerosol import esm_sm_aerosol_save_dirpath
+from cmip7_ancil_argparse import common_parser
+from cmip7_ancil_common import (
+    INTERPOLATION_SCHEME,
+    cmip7_date_constraint_from_years,
+    esm_grid_mask_cube,
+    extend_years,
+    save_ancil,
+    set_coord_system,
 )
-from cmip7_ancil_argparse import (
-    common_parser,
-    percent_parser,
-)
-from cmip7_ancil_common import cmip7_date_constraint_from_years
+from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
 
 
 def parse_args():
     parser = ArgumentParser(
         prog="cmip7_SM_Bio_interpolate",
         description=(
-            "Generate input files from CMIP7 historical biomass forcings"
+            "Generate input files from CMIP7 ScenarioMIP biomass forcings"
         ),
-        parents=[
-            common_parser(),
-            percent_parser(),
-        ],
+        parents=[common_parser()],
     )
     parser.add_argument("--scenario")
     parser.add_argument("--dataset-date-range")
@@ -35,28 +37,154 @@ def parse_args():
     return parser.parse_args()
 
 
-def load_cmip7_sm_aerosol_biomass(args, species):
-    return load_cmip7_aerosol_biomass(
-        args,
-        species,
-        args.dataset_date_range,
-        cmip7_date_constraint_from_years(
-            CMIP7_SM_AEROSOL_BEG_YEAR,
-            CMIP7_SM_AEROSOL_END_YEAR,
-        ),
+def _biomass_dirpath(args, species):
+    return (
+        Path(args.cmip7_source_data_dirname)
+        / "ScenarioMIP"
+        / "IIASA-IAMC"
+        / args.dataset_version
+        / "atmos"
+        / "mon"
+        / species
+        / "gn"
+        / args.dataset_vdate
+    )
+
+
+def cmip7_sm_aerosol_biomass_filepath(args, species, date_range):
+    dirpath = _biomass_dirpath(args, species)
+    filename = (
+        f"{species}_input4MIPs_emissions_ScenarioMIP_"
+        f"{args.dataset_version}_gn_"
+        f"{date_range}.nc"
+    )
+    return dirpath / filename
+
+
+def load_cmip7_aerosol_biomass(args, species, date_range, constraint):
+    cube = load_cmip7_aerosol(
+        args, cmip7_sm_aerosol_biomass_filepath, species, date_range, constraint
+    )
+    # This data is missing over oceans,
+    # so needs to be filled with zero for the model
+    cube.data = cube.data.filled(0.0)
+    return cube
+
+
+def load_sector_dict(args, filepath_fn):
+    # Iris doesn't read the sector coordinate so use netCDF4
+    date_range = args.dataset_date_range
+    d = netCDF4.Dataset(filepath_fn(args, "BC", date_range))
+    sectord = dict()
+    for s in d["sector"].ids.split(";"):
+        i, name = s.split(":")
+        sectord[name.strip()] = int(i)
+    d.close()
+    return sectord
+
+
+force_load = True
+
+
+def split_sm_low_high(bc_cube, oc_cube, sector_dict):
+    """
+    Based on code generated by Gemini Flash 3.5.
+
+    Slices and sums the 4D ScenarioMIP BC and OC emission cubes to partition
+    low-level and high-level injection height emissions.
+
+    Dimensions of input cubes: (time, sector, lat, lon)
+    Sectors:
+    0: Agricultural Waste Burning;
+    1: Forest Burning;
+    2: Grassland Burning;
+    3: Peat Burning
+    """
+    agri = sector_dict["Agricultural Waste Burning"]
+    fore = sector_dict["Forest Burning"]
+    sava = sector_dict["Grassland Burning"]
+    peat = sector_dict["Peat Burning"]
+
+    # Sum low-level sectors using collapsed:
+    low_sectors = [agri, sava, peat]
+    bc_low = bc_cube[:, low_sectors].collapsed("sector", iris.analysis.SUM)
+    oc_low = oc_cube[:, low_sectors].collapsed("sector", iris.analysis.SUM)
+
+    # Extract high-level sector:
+    bc_high = bc_cube[:, fore]
+    oc_high = oc_cube[:, fore]
+
+    # Remove the scalar sector coordinate from the sliced cubes to prevent
+    # coordinate mismatch warnings when summing or saving
+    for cube in (bc_low, oc_low, bc_high, oc_high):
+        if cube.coords("sector"):
+            cube.remove_coord("sector")
+
+    # Combine the species into bulk low and high emission cubes
+    low = bc_low + oc_low
+    high = bc_high + oc_high
+
+    return low, high
+
+
+def save_cmip7_aerosol_biomass(args, filepath_fn, load_fn, save_dirpath):
+    bc = load_fn(args, "BC")
+    oc = load_fn(args, "OC")
+
+    sector_dict = load_sector_dict(args, filepath_fn)
+
+    low, high = split_sm_low_high(bc, oc, sector_dict)
+
+    if force_load:
+        _ = low.data
+        _ = high.data
+        now = datetime.now()
+        print(f"{now}: LO, HI done")
+
+    # Regrid requires matching coordinate systems
+    set_coord_system(low)
+    set_coord_system(high)
+
+    now = datetime.now()
+    print(f"{now}: set_coord_system done")
+
+    esm_grid_mask = esm_grid_mask_cube(args)
+    low_esm = low.regrid(esm_grid_mask, INTERPOLATION_SCHEME)
+    high_esm = high.regrid(esm_grid_mask, INTERPOLATION_SCHEME)
+
+    now = datetime.now()
+    print(f"{now}: regrid done")
+
+    zero_poles(low_esm)
+    zero_poles(high_esm)
+
+    now = datetime.now()
+    print(f"{now}: zero_poles done")
+
+    low_esm.attributes["STASH"] = iris.fileformats.pp.STASH(
+        model=1, section=0, item=130
     )
+    high_esm.attributes["STASH"] = iris.fileformats.pp.STASH(
+        model=1, section=0, item=131
+    )
+
+    save_ancil([low_esm, high_esm], save_dirpath, args.save_filename)
+
+    now = datetime.now()
+    print(f"{now}: save_ancil done")
 
 
-def load_cmip7_sm_aerosol_biomass_percentage(args, species):
-    return load_cmip7_aerosol_biomass(
+def load_cmip7_sm_aerosol_biomass(args, species):
+    cube = load_cmip7_aerosol_biomass(
         args,
         species,
-        args.percent_date_range,
+        args.dataset_date_range,
         cmip7_date_constraint_from_years(
-            CMIP7_SM_AEROSOL_BEG_YEAR,
-            CMIP7_SM_AEROSOL_END_YEAR,
+            CMIP7_SM_BEG_YEAR,
+            CMIP7_SM_END_YEAR,
         ),
     )
+    return extend_years(cube)
 
 
 if __name__ == "__main__":
@@ -64,7 +192,7 @@ def load_cmip7_sm_aerosol_biomass_percentage(args, species):
 
     save_cmip7_aerosol_biomass(
         args,
-        load_cmip7_sm_aerosol_biomass_percentage,
+        cmip7_sm_aerosol_biomass_filepath,
         load_cmip7_sm_aerosol_biomass,
         esm_sm_aerosol_save_dirpath(args),
     )

From eaf44e3e71862ace87e0575f3eb1683e9354b72d Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 11:56:56 +1000
Subject: [PATCH 36/50] Add and use the _biomass_variable function

---
 .../esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index 9acef9d..291f657 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -37,6 +37,10 @@ def parse_args():
     return parser.parse_args()
 
 
+def _biomass_variable(species):
+    return f"{species}-em-openburning"
+
+
 def _biomass_dirpath(args, species):
     return (
         Path(args.cmip7_source_data_dirname)
@@ -45,7 +49,7 @@ def _biomass_dirpath(args, species):
         / args.dataset_version
         / "atmos"
         / "mon"
-        / species
+        / _biomass_variable(species)
         / "gn"
         / args.dataset_vdate
     )
@@ -54,7 +58,7 @@ def _biomass_dirpath(args, species):
 def cmip7_sm_aerosol_biomass_filepath(args, species, date_range):
     dirpath = _biomass_dirpath(args, species)
     filename = (
-        f"{species}_input4MIPs_emissions_ScenarioMIP_"
+        f"{_biomass_variable(species)}_input4MIPs_emissions_ScenarioMIP_"
         f"{args.dataset_version}_gn_"
         f"{date_range}.nc"
     )

From 500ee240d80ef3f81bf8f3e8555f59b9c8b049ac Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 12:04:17 +1000
Subject: [PATCH 37/50] Correctly use _biomass_variable

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index 291f657..a979c55 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -49,7 +49,7 @@ def _biomass_dirpath(args, species):
         / args.dataset_version
         / "atmos"
         / "mon"
-        / _biomass_variable(species)
+        / _biomass_variable(species).replace("-", "_")
         / "gn"
         / args.dataset_vdate
     )

From b2f048a832babe65b5d273e6bc2f152e39960a31 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 13:18:41 +1000
Subject: [PATCH 38/50] Use interpolate_monthly

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
index a979c55..b3b3abb 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py
@@ -17,6 +17,7 @@
     cmip7_date_constraint_from_years,
     esm_grid_mask_cube,
     extend_years,
+    interpolate_monthly,
     save_ancil,
     set_coord_system,
 )
@@ -188,7 +189,10 @@ def load_cmip7_sm_aerosol_biomass(args, species):
             CMIP7_SM_END_YEAR,
         ),
     )
-    return extend_years(cube)
+    interpolated = interpolate_monthly(
+        cube, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+    )
+    return extend_years(interpolated)
 
 
 if __name__ == "__main__":

From d54eb5212061bf281d5c0bc1276642f5e5c1aa7f Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 13:19:20 +1000
Subject: [PATCH 39/50] Remove unused code and parameters

---
 .../aerosol/cmip7_SM_aerosol_anthro.py        | 43 +------------------
 1 file changed, 2 insertions(+), 41 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
index 976231b..6113b23 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
@@ -3,11 +3,7 @@
 
 from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_interpolate
 from aerosol.cmip7_aerosol_common import load_cmip7_aerosol
-from aerosol.cmip7_SM_aerosol import (
-    CMIP7_SM_AEROSOL_BEG_YEAR,
-    CMIP7_SM_AEROSOL_END_YEAR,
-    esm_sm_aerosol_save_dirpath,
-)
+from aerosol.cmip7_SM_aerosol import esm_sm_aerosol_save_dirpath
 from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
 from cmip7_ancil_argparse import common_parser
 from cmip7_ancil_common import (
@@ -46,16 +42,6 @@ def _anthro_dirpath(args, variable):
     )
 
 
-def cmip7_sm_aerosol_air_anthro_filepath(args, species, date_range):
-    dirpath = _anthro_dirpath(args, f"{species}_em_AIR_anthro")
-    filename = (
-        f"{species}-em-AIR-anthro_input4MIPs_emissions_ScenarioMIP_"
-        f"{args.dataset_version}_gn_"
-        f"{date_range}.nc"
-    )
-    return dirpath / filename
-
-
 def cmip7_sm_aerosol_anthro_filepath(args, species, date_range):
     dirpath = _anthro_dirpath(args, f"{species}_em_anthro")
     filename = (
@@ -66,32 +52,7 @@ def cmip7_sm_aerosol_anthro_filepath(args, species, date_range):
     return dirpath / filename
 
 
-def load_cmip7_sm_aerosol_air_anthro(
-    args,
-    species,
-    beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
-    end_year=CMIP7_SM_AEROSOL_END_YEAR,
-):
-    cube = load_cmip7_aerosol(
-        args,
-        cmip7_sm_aerosol_air_anthro_filepath,
-        species,
-        args.dataset_date_range,
-        cmip7_date_constraint_from_years(CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR),
-    )
-    fix_coords(args, cube)
-    interpolated = interpolate_monthly(
-        cube, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
-    )
-    return extend_years(interpolated)
-
-
-def load_cmip7_sm_aerosol_anthro(
-    args,
-    species,
-    beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
-    end_year=CMIP7_SM_AEROSOL_END_YEAR,
-):
+def load_cmip7_sm_aerosol_anthro(args, species):
     cube = load_cmip7_aerosol(
         args,
         cmip7_sm_aerosol_anthro_filepath,

From 1fd24a7bae2d384e4177f2d392a13a76605d361b Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 14:11:22 +1000
Subject: [PATCH 40/50] Remove unused arguments

---
 .../atmosphere/aerosol/cmip7_SM_SO2_interpolate.py  | 13 ++-----------
 1 file changed, 2 insertions(+), 11 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
index 9421ceb..00a83ad 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py
@@ -4,11 +4,7 @@
 
 from aerosol.cmip7_PI_aerosol import esm_pi_aerosol_ancil_dirpath
 from aerosol.cmip7_PI_SO2_interpolate import PI_DMS_ANCIL_FILENAME
-from aerosol.cmip7_SM_aerosol import (
-    CMIP7_SM_AEROSOL_BEG_YEAR,
-    CMIP7_SM_AEROSOL_END_YEAR,
-    esm_sm_aerosol_save_dirpath,
-)
+from aerosol.cmip7_SM_aerosol import esm_sm_aerosol_save_dirpath
 from aerosol.cmip7_SM_aerosol_anthro import (
     cmip7_sm_aerosol_anthro_filepath,
     load_cmip7_sm_aerosol_anthro,
@@ -42,12 +38,7 @@ def parse_args():
 
 
 def load_cmip7_sm_so2_aerosol_anthro(args, species):
-    return load_cmip7_sm_aerosol_anthro(
-        args,
-        species,
-        beg_year=CMIP7_SM_AEROSOL_BEG_YEAR,
-        end_year=CMIP7_SM_AEROSOL_END_YEAR,
-    )
+    return load_cmip7_sm_aerosol_anthro(args, species)
 
 
 def load_sm_dms(args):

From 7f9d6d3f42290759a2b2fb29def6053adcbb9508 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 14:11:53 +1000
Subject: [PATCH 41/50] In zero_poles, make data writeable

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 3e1e7aa..2a7ea5e 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -232,6 +232,9 @@ def zero_poles(cube):
     # For aerosol emissions they should be zero
     latdim = cube.coord_dims("latitude")
     assert latdim == (1,)
+    # Make data writeable
+    if not cube.data.flags.writeable:
+        cube.data = cube.data.copy()
     cube.data[:, 0] = 0.0
     cube.data[:, -1] = 0.0
 

From 2b89b75da41656d778ab1cd2bab1e4264cc0d49b Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 14:49:01 +1000
Subject: [PATCH 42/50] Move the corrected zero_poles to
 cmip7_aerosol_common.py

---
 .../aerosol/cmip7_aerosol_common.py           |  3 +++
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 24 ++++---------------
 2 files changed, 8 insertions(+), 19 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py
index 5425e1a..1d1bd8f 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py
@@ -32,5 +32,8 @@ def zero_poles(cube):
     # For aerosol emissions they should be zero
     latdim = cube.coord_dims("latitude")
     assert latdim == (1,)
+    # Make data writeable
+    if not cube.data.flags.writeable:
+        cube.data = cube.data.copy()
     cube.data[:, 0] = 0.0
     cube.data[:, -1] = 0.0
diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 2a7ea5e..a7637ad 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -151,7 +151,9 @@ def interpolate_monthly(cube, beg_year, end_year):
             if month == MONTHS_IN_A_YEAR:
                 end_date = cftime.datetime(year + 1, 1, 1, calendar=calendar)
             else:
-                end_date = cftime.datetime(year, month + 1, 1, calendar=calendar)
+                end_date = cftime.datetime(
+                    year, month + 1, 1, calendar=calendar
+                )
 
             mid_date = beg_date + (end_date - beg_date) / 2
             tdates.append(mid_date)
@@ -164,9 +166,7 @@ def interpolate_monthly(cube, beg_year, end_year):
     )
 
     # Perform linear time-interpolation
-    new_cube = cube.interpolate(
-        [("time", tpoints)], iris.analysis.Linear()
-    )
+    new_cube = cube.interpolate([("time", tpoints)], iris.analysis.Linear())
 
     # Create new DimCoord with contiguous bounds
     new_time_coord = iris.coords.DimCoord(
@@ -200,9 +200,7 @@ def interpolate_monthly(cube, beg_year, end_year):
             beg_date = units.num2date(bounds[i][0])
             end_date = units.num2date(bounds[i][1])
             # Match bounds based on start and end year/month components
-            key = (
-                beg_date.year, beg_date.month, end_date.year, end_date.month
-            )
+            key = (beg_date.year, beg_date.month, end_date.year, end_date.month)
             if key in tbounds_dict:
                 target_idx = tbounds_dict[key]
                 new_cube.data[target_idx] = cube.data[i]
@@ -227,18 +225,6 @@ def fix_coords(args, cube):
     ).coord_system
 
 
-def zero_poles(cube):
-    # Polar values should have no longitude dependence
-    # For aerosol emissions they should be zero
-    latdim = cube.coord_dims("latitude")
-    assert latdim == (1,)
-    # Make data writeable
-    if not cube.data.flags.writeable:
-        cube.data = cube.data.copy()
-    cube.data[:, 0] = 0.0
-    cube.data[:, -1] = 0.0
-
-
 def save_ancil(
     cubes, save_dirpath, save_filename, gregorian=True, replace_bounds=False
 ):

From 97c57c93337a220a3f31f5364261ca550eb1ec9c Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 15:33:41 +1000
Subject: [PATCH 43/50] Modify a copy of data then write it back

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py
index 1d1bd8f..09cc010 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py
@@ -32,8 +32,7 @@ def zero_poles(cube):
     # For aerosol emissions they should be zero
     latdim = cube.coord_dims("latitude")
     assert latdim == (1,)
-    # Make data writeable
-    if not cube.data.flags.writeable:
-        cube.data = cube.data.copy()
-    cube.data[:, 0] = 0.0
-    cube.data[:, -1] = 0.0
+    data = cube.data.copy()
+    data[:, 0] = 0.0
+    data[:, -1] = 0.0
+    cube.data = data

From a947a4bd822dc56a6f3684c07920922797d44bd9 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 16:15:58 +1000
Subject: [PATCH 44/50] Try setting CMIP7_SM_BEG_YEAR to 2023

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py | 4 ----
 CMIP7/esm1p6/atmosphere/cmip7_SM.py                 | 4 +---
 2 files changed, 1 insertion(+), 7 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
index 21620e2..03d3cfb 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py
@@ -1,10 +1,6 @@
 from pathlib import Path
 
 from cmip7_ancil_constants import ANCIL_TODAY
-from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
-
-CMIP7_SM_AEROSOL_BEG_YEAR = CMIP7_SM_BEG_YEAR - 1
-CMIP7_SM_AEROSOL_END_YEAR = CMIP7_SM_END_YEAR + 1
 
 
 def esm_sm_aerosol_ancil_dirpath(args):
diff --git a/CMIP7/esm1p6/atmosphere/cmip7_SM.py b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
index 44ee485..094c206 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_SM.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
@@ -2,10 +2,8 @@
 
 from cmip7_ancil_constants import ANCIL_TODAY
 
-CMIP7_SM_BEG_YEAR = 2022
-# Model time interpolation requires an extra year
+CMIP7_SM_BEG_YEAR = 2023
 CMIP7_SM_END_YEAR = 2100
-CMIP7_SM_NBR_YEARS = CMIP7_SM_END_YEAR + 1 - CMIP7_SM_BEG_YEAR
 
 
 def esm_sm_forcing_save_dirpath(args):

From 3933ba7b33f801f96a0b47ba0b55bfa77a560861 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Mon, 15 Jun 2026 18:59:09 +1000
Subject: [PATCH 45/50] Satisfy pre-commit

---
 CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
index 6113b23..6b3aafb 100644
--- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
+++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py
@@ -4,7 +4,6 @@
 from aerosol.cmip7_aerosol_anthro import cmip7_aerosol_anthro_interpolate
 from aerosol.cmip7_aerosol_common import load_cmip7_aerosol
 from aerosol.cmip7_SM_aerosol import esm_sm_aerosol_save_dirpath
-from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
 from cmip7_ancil_argparse import common_parser
 from cmip7_ancil_common import (
     cmip7_date_constraint_from_years,
@@ -12,6 +11,7 @@
     fix_coords,
     interpolate_monthly,
 )
+from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
 
 
 def parse_args(species):

From 8853f0b49d4130fb407ef65071dbf93b2a8dc8f9 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Wed, 17 Jun 2026 09:32:09 +1000
Subject: [PATCH 46/50] Try changing CMIP7_SM_BEG_YEAR back to 2022

---
 CMIP7/esm1p6/atmosphere/cmip7_SM.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_SM.py b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
index 094c206..060215c 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_SM.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_SM.py
@@ -2,7 +2,7 @@
 
 from cmip7_ancil_constants import ANCIL_TODAY
 
-CMIP7_SM_BEG_YEAR = 2023
+CMIP7_SM_BEG_YEAR = 2022
 CMIP7_SM_END_YEAR = 2100
 
 

From 04cf8d48381b3361d99bf6cd7fada201f9c9aa97 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Wed, 17 Jun 2026 10:11:39 +1000
Subject: [PATCH 47/50] Start ScenarioMIP greenhouse gas processing one year
 late

---
 CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
index f26512d..e2c6bac 100644
--- a/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
+++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py
@@ -28,13 +28,14 @@ def parse_args():
 if __name__ == "__main__":
     args = parse_args()
 
+    CMIP7_SM_GHG_BEG_YEAR = CMIP7_SM_BEG_YEAR + 1
     ghg_mmr_dict = dict()
     for ghg in GHG_MOLAR_MASS:
         ghg_mmr_dict[ghg] = load_cmip7_ghg_series_mmr(
-            args, "ScenarioMIP", ghg, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+            args, "ScenarioMIP", ghg, CMIP7_SM_GHG_BEG_YEAR, CMIP7_SM_END_YEAR
         )
 
     # Patch the greenhouse gas namelist.
     cmip7_ghg_update_namelists_file(
-        ghg_mmr_dict, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+        ghg_mmr_dict, CMIP7_SM_GHG_BEG_YEAR, CMIP7_SM_END_YEAR
     )

From 614a9908a3dc0318ed1c1bd3c679a7084b6d3a59 Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Wed, 17 Jun 2026 18:39:16 +1000
Subject: [PATCH 48/50] Interpolate each month separately

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 22 +++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index a7637ad..14fdc3f 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -9,6 +9,7 @@
 import cftime
 import iris
 import iris.analysis
+import iris.coord_categorisation
 import mule
 import numpy as np
 from cmip7_ancil_constants import (
@@ -167,6 +168,27 @@ def interpolate_monthly(cube, beg_year, end_year):
 
     # Perform linear time-interpolation
     new_cube = cube.interpolate([("time", tpoints)], iris.analysis.Linear())
+    new_cube.data = np.ma.asarray(new_cube.data)
+
+    # Add month categorisation to extract each month's series
+    if not cube.coords("month_number"):
+        iris.coord_categorisation.add_month_number(cube, "time", name="month_number")
+
+    # Interpolate each month separately and interleave the data
+    for m in range(1, MONTHS_IN_A_YEAR + 1):
+        month_constraint = iris.Constraint(month_number=m)
+        m_cube = cube.extract(month_constraint)
+        if m_cube is not None:
+            # Select target time points for this month across all years
+            m_tpoints = tpoints[m - 1 :: MONTHS_IN_A_YEAR]
+            # Interpolate across the years for this month
+            m_interpolated = m_cube.interpolate([("time", m_tpoints)], iris.analysis.Linear())
+            # Place the interpolated data back into the interleaved target indices
+            new_cube.data[m - 1 :: MONTHS_IN_A_YEAR] = m_interpolated.data
+
+    # Clean up month_number coordinate if added
+    if cube.coords("month_number"):
+        cube.remove_coord("month_number")
 
     # Create new DimCoord with contiguous bounds
     new_time_coord = iris.coords.DimCoord(

From 1b1bf8adc2c86355160185796ee988fb185aa12e Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Wed, 17 Jun 2026 19:02:51 +1000
Subject: [PATCH 49/50] Tidy up for ruff

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index 14fdc3f..fb750a9 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -172,7 +172,9 @@ def interpolate_monthly(cube, beg_year, end_year):
 
     # Add month categorisation to extract each month's series
     if not cube.coords("month_number"):
-        iris.coord_categorisation.add_month_number(cube, "time", name="month_number")
+        iris.coord_categorisation.add_month_number(
+            cube, "time", name="month_number"
+        )
 
     # Interpolate each month separately and interleave the data
     for m in range(1, MONTHS_IN_A_YEAR + 1):
@@ -182,8 +184,11 @@ def interpolate_monthly(cube, beg_year, end_year):
             # Select target time points for this month across all years
             m_tpoints = tpoints[m - 1 :: MONTHS_IN_A_YEAR]
             # Interpolate across the years for this month
-            m_interpolated = m_cube.interpolate([("time", m_tpoints)], iris.analysis.Linear())
-            # Place the interpolated data back into the interleaved target indices
+            m_interpolated = m_cube.interpolate(
+                [("time", m_tpoints)], iris.analysis.Linear()
+            )
+            # Place the interpolated data back
+            # into the interleaved target indices
             new_cube.data[m - 1 :: MONTHS_IN_A_YEAR] = m_interpolated.data
 
     # Clean up month_number coordinate if added

From 57a03051d2a4f9954f85e85fe42bf08101b97c5b Mon Sep 17 00:00:00 2001
From: Paul Leopardi <paul.leopardi@anu.edu.au>
Date: Wed, 17 Jun 2026 19:23:28 +1000
Subject: [PATCH 50/50] Refactor interpolate_monthly

---
 CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py | 62 ++++++++++---------
 1 file changed, 34 insertions(+), 28 deletions(-)

diff --git a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
index fb750a9..e1317d9 100644
--- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
+++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py
@@ -137,6 +137,39 @@ def extend_years(cube):
     return cubelist.concatenate_cube()
 
 
+def _interpolate_months_separately(cube, tpoints):
+    # Perform linear time-interpolation
+    new_cube = cube.interpolate([("time", tpoints)], iris.analysis.Linear())
+    new_cube.data = np.ma.asarray(new_cube.data)
+
+    # Add month categorisation to extract each month's series
+    if not cube.coords("month_number"):
+        iris.coord_categorisation.add_month_number(
+            cube, "time", name="month_number"
+        )
+
+    # Interpolate each month separately and interleave the data
+    for m in range(1, MONTHS_IN_A_YEAR + 1):
+        month_constraint = iris.Constraint(month_number=m)
+        m_cube = cube.extract(month_constraint)
+        if m_cube is not None:
+            # Select target time points for this month across all years
+            m_tpoints = tpoints[m - 1 :: MONTHS_IN_A_YEAR]
+            # Interpolate across the years for this month
+            m_interpolated = m_cube.interpolate(
+                [("time", m_tpoints)], iris.analysis.Linear()
+            )
+            # Place the interpolated data back
+            # into the interleaved target indices
+            new_cube.data[m - 1 :: MONTHS_IN_A_YEAR] = m_interpolated.data
+
+    # Clean up month_number coordinate if added
+    if cube.coords("month_number"):
+        cube.remove_coord("month_number")
+
+    return new_cube
+
+
 def interpolate_monthly(cube, beg_year, end_year):
     # Get original time units and calendar
     time_coord = cube.coord("time")
@@ -166,34 +199,8 @@ def interpolate_monthly(cube, beg_year, end_year):
         dtype=np.float64,
     )
 
-    # Perform linear time-interpolation
-    new_cube = cube.interpolate([("time", tpoints)], iris.analysis.Linear())
-    new_cube.data = np.ma.asarray(new_cube.data)
-
-    # Add month categorisation to extract each month's series
-    if not cube.coords("month_number"):
-        iris.coord_categorisation.add_month_number(
-            cube, "time", name="month_number"
-        )
-
     # Interpolate each month separately and interleave the data
-    for m in range(1, MONTHS_IN_A_YEAR + 1):
-        month_constraint = iris.Constraint(month_number=m)
-        m_cube = cube.extract(month_constraint)
-        if m_cube is not None:
-            # Select target time points for this month across all years
-            m_tpoints = tpoints[m - 1 :: MONTHS_IN_A_YEAR]
-            # Interpolate across the years for this month
-            m_interpolated = m_cube.interpolate(
-                [("time", m_tpoints)], iris.analysis.Linear()
-            )
-            # Place the interpolated data back
-            # into the interleaved target indices
-            new_cube.data[m - 1 :: MONTHS_IN_A_YEAR] = m_interpolated.data
-
-    # Clean up month_number coordinate if added
-    if cube.coords("month_number"):
-        cube.remove_coord("month_number")
+    new_cube = _interpolate_months_separately(cube, tpoints)
 
     # Create new DimCoord with contiguous bounds
     new_time_coord = iris.coords.DimCoord(
@@ -232,7 +239,6 @@ def interpolate_monthly(cube, beg_year, end_year):
                 target_idx = tbounds_dict[key]
                 new_cube.data[target_idx] = cube.data[i]
 
-    new_cube.data = np.ma.asarray(new_cube.data)
     return new_cube