Create pre-industrial and historical ozone forcing ancillaries from UK ozone CMIP7 NetCDF files (#111)

* Code to obtain ESM1.6 ozone from UKESM NetCDF files

* Update ozone code after first tests

* Do not regrid

* Do not use the Gregorian calendar

* Fix a typo in a comment

* Add the PI_mean ozone script to the collection of ozone scripts (#122)

* Add a PI_mean script for 1850-1870 mean as PI control

* Remove the added 'month' coordinate before saving

* Guess time bounds to ensure that time is contiguous

* Set the calendar to Gregorian but replace the bounds (#125)

* Save using the Gregorian calendar

* Try replacing time bounds

* Indent correctly

Author: penguian

CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py

@@ -45,7 +45,7 @@ def esm_grid_mask_cube(args):
     return cube
 
 
-def set_gregorian(var):
+def set_gregorian(var, replace_bounds=False):
     # Change the calendar to Gregorian for the model
     time = var.coord("time")
     origin = time.units.origin
@@ -59,17 +59,43 @@ def set_gregorian(var):
             date.year, date.month, date.day, date.hour, date.minute, date.second
         )
         tvals[i] = newunits.date2num(newdate)
-        for j in range(2):
-            date = time.units.num2date(tbnds[i][j])
-            newdate = cftime.DatetimeProlepticGregorian(
+        if replace_bounds:
+            beg_date = cftime.DatetimeProlepticGregorian(
                 date.year,
                 date.month,
-                date.day,
+                1,
                 date.hour,
                 date.minute,
                 date.second,
             )
-            tbnds[i][j] = newunits.date2num(newdate)
+            tbnds[i][0] = newunits.date2num(beg_date)
+            if date.month == 12:
+                end_year = date.year + 1
+                end_month = 1
+            else:
+                end_year = date.year
+                end_month = date.month + 1
+            end_date = cftime.DatetimeProlepticGregorian(
+                end_year,
+                end_month,
+                1,
+                date.hour,
+                date.minute,
+                date.second,
+            )
+            tbnds[i][1] = newunits.date2num(end_date)
+        else:
+            for j in range(2):
+                date = time.units.num2date(tbnds[i][j])
+                newdate = cftime.DatetimeProlepticGregorian(
+                    date.year,
+                    date.month,
+                    date.day,
+                    date.hour,
+                    date.minute,
+                    date.second,
+                )
+                tbnds[i][j] = newunits.date2num(newdate)
     time.points = tvals
     time.bounds = tbnds
     time.units = newunits
@@ -134,7 +160,9 @@ def zero_poles(cube):
     cube.data[:, -1] = 0.0
 
 
-def save_ancil(cubes, save_dirpath, save_filename):
+def save_ancil(
+    cubes, save_dirpath, save_filename, gregorian=True, replace_bounds=False
+):
     """
     Handle both a list and a single cube
     """
@@ -146,8 +174,9 @@ def save_ancil(cubes, save_dirpath, save_filename):
     """
     for cube in cubes:
         cube.attributes["grid_staggering"] = 3  # New dynamics
-        cube.attributes["time_type"] = 1  # Gregorian
-        set_gregorian(cube)
+        if gregorian:
+            cube.attributes["time_type"] = 1  # Gregorian
+            set_gregorian(cube, replace_bounds=replace_bounds)
     """
     ANTS doesn't set the calendar header for monthly fields
     See fileformats/ancil/time_headers.py

CMIP7/esm1p6/atmosphere/ozone/cmip7_HI_ozone_generate.py

@@ -0,0 +1,53 @@
+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_HI_ozone_generate",
+        description=(
+            "Generate input files from UK CMIP7 historical ozone forcings"
+        ),
+    )
+    return parser.parse_args()
+
+
+def esm_hi_ozone_save_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "modern"
+        / "historical"
+        / "forcing"
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )
+
+
+def save_cmip7_hi_ozone(args, cube):
+    # Save as an ancillary file
+    save_dirpath = esm_hi_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_hi_ozone(args, esm_cube)

CMIP7/esm1p6/atmosphere/ozone/cmip7_PI_mean_ozone_generate.py

@@ -0,0 +1,67 @@
+from argparse import ArgumentParser
+from pathlib import Path
+
+import iris.coord_categorisation
+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_PI_mean_ozone_generate",
+        description=(
+            "Generate input files from 21 year mean of UK CMIP7 historical"
+            "ozone forcings"
+        ),
+    )
+    return parser.parse_args()
+
+
+def esm_pi_ozone_save_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "modern"
+        / "pre-industrial-mean"
+        / "forcing"
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )
+
+
+def save_cmip7_pi_mean_ozone(args, cube):
+    # Add a "month_number" variable.
+    iris.coord_categorisation.add_month_number(
+        cube, "time", name="month_number"
+    )
+    # Aggregate by "month_number" to obtain a mean for each month.
+    cube = cube.aggregated_by("month_number", iris.analysis.MEAN)
+    # Re-create the time bounds to ensure that time is contiguous.
+    time = cube.coord("time")
+    time.bounds = None
+    time.guess_bounds(bound_position=0.5)
+    # Remove the added "month_number" coordinate before saving.
+    cube.remove_coord("month_number")
+    # Save as an ancillary file
+    save_dirpath = esm_pi_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_pi_mean_ozone(args, esm_cube)

CMIP7/esm1p6/atmosphere/ozone/cmip7_PI_ozone_generate.py

@@ -0,0 +1,53 @@
+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_PI_ozone_generate",
+        description=(
+            "Generate input files from UK CMIP7 pre-industrial ozone forcings"
+        ),
+    )
+    return parser.parse_args()
+
+
+def esm_pi_ozone_save_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "modern"
+        / "pre-industrial"
+        / "forcing"
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )
+
+
+def save_cmip7_pi_ozone(args, cube):
+    # Save as an ancillary file
+    save_dirpath = esm_pi_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_pi_ozone(args, esm_cube)

CMIP7/esm1p6/atmosphere/ozone/cmip7_ozone.py

@@ -0,0 +1,32 @@
+from argparse import ArgumentParser
+from pathlib import Path
+
+import iris
+from cmip7_ancil_common import fix_coords
+
+
+def ozone_parser():
+    parser = ArgumentParser(add_help=False)
+    parser.add_argument("--ukesm-ancil-dirpath")
+    parser.add_argument("--ukesm-netcdf-filename")
+    parser.add_argument("--save-filename")
+    return parser
+
+
+def cmip7_ozone_filepath(args):
+    dirpath = Path(args.ukesm_ancil_dirpath)
+    filename = args.ukesm_netcdf_filename
+    return dirpath / filename
+
+
+def load_cmip7_ozone(args):
+    filepath = cmip7_ozone_filepath(args)
+    cube = iris.load_cube(filepath)
+    return cube
+
+
+def fix_cmip7_ozone(args, cube):
+    # Make the coordinates compatible with the ESM1.5 grid mask
+    fix_coords(args, cube)
+    cube.data = cube.data.filled(0.0)
+    return cube