Python code for aerosol, greenhouse gas, nitrogen, ozone, solar and volcanic forcings for ScenarioMIP for ACCESS-ESM1.6

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,
@@ -20,16 +21,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)
@@ -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),

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,
@@ -14,18 +15,18 @@
 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=(
             "Generate input files from CMIP7 pre-industrial 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")
@@ -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),

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)

CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_Bio_interpolate.py

@@ -0,0 +1,206 @@
+# Interpolate CMIP7 SM Biomass burning emissions to ESM1.6 grid
+
+from argparse import ArgumentParser
+from datetime import datetime
+from pathlib import Path
+
+import iris
+import netCDF4
+from aerosol.cmip7_aerosol_common import (
+    load_cmip7_aerosol,
+    zero_poles,
+)
+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,
+    interpolate_monthly,
+    save_ancil,
+    set_coord_system,
+)
+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 ScenarioMIP biomass forcings"
+        ),
+        parents=[common_parser()],
+    )
+    parser.add_argument("--scenario")
+    parser.add_argument("--dataset-date-range")
+    parser.add_argument("--save-filename")
+    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)
+        / "ScenarioMIP"
+        / "IIASA-IAMC"
+        / args.dataset_version
+        / "atmos"
+        / "mon"
+        / _biomass_variable(species).replace("-", "_")
+        / "gn"
+        / args.dataset_vdate
+    )
+
+
+def cmip7_sm_aerosol_biomass_filepath(args, species, date_range):
+    dirpath = _biomass_dirpath(args, species)
+    filename = (
+        f"{_biomass_variable(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(args, species):
+    cube = load_cmip7_aerosol_biomass(
+        args,
+        species,
+        args.dataset_date_range,
+        cmip7_date_constraint_from_years(
+            CMIP7_SM_BEG_YEAR,
+            CMIP7_SM_END_YEAR,
+        ),
+    )
+    interpolated = interpolate_monthly(
+        cube, CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR
+    )
+    return extend_years(interpolated)
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    save_cmip7_aerosol_biomass(
+        args,
+        cmip7_sm_aerosol_biomass_filepath,
+        load_cmip7_sm_aerosol_biomass,
+        esm_sm_aerosol_save_dirpath(args),
+    )

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)

CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_SO2_interpolate.py

@@ -0,0 +1,64 @@
+# Interpolate CMIP7 HI SO2 emissions to ESM1.6 grid
+
+from argparse import ArgumentParser
+
+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 esm_sm_aerosol_save_dirpath
+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,
+)
+from cmip7_ancil_argparse import (
+    common_parser,
+    dms_filename_parser,
+)
+from cmip7_PI import fix_esm15_pi_ancil_date
+
+
+def parse_args():
+    parser = ArgumentParser(
+        prog="cmip7_SM_SO2_interpolate",
+        description=(
+            "Generate input files from CMIP7 ScenarioMIP SO2 forcings"
+        ),
+        parents=[
+            common_parser(),
+            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()
+
+
+def load_cmip7_sm_so2_aerosol_anthro(args, species):
+    return load_cmip7_sm_aerosol_anthro(args, species)
+
+
+def load_sm_dms(args):
+    # Use the CMIP6 DMS
+    dms_ancil_dirpath = (
+        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_pi_ancil_date)
+
+
+if __name__ == "__main__":
+    args = parse_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),
+    )

CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol.py

@@ -0,0 +1,21 @@
+from pathlib import Path
+
+from cmip7_ancil_constants import ANCIL_TODAY
+
+
+def esm_sm_aerosol_ancil_dirpath(args):
+    return (
+        Path(args.ancil_target_dirname)
+        / "scenarios"
+        / args.scenario
+        / "atmosphere"
+        / "aerosol"
+    )
+
+
+def esm_sm_aerosol_save_dirpath(args):
+    return (
+        esm_sm_aerosol_ancil_dirpath(args)
+        / args.esm_grid_rel_dirname
+        / ANCIL_TODAY
+    )