diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_HI_SO2_interpolate.py index e3aea5d..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, @@ -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), diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_PI_SO2_interpolate.py index ffe8ee8..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, @@ -14,10 +15,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 +26,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") @@ -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_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..b3b3abb --- /dev/null +++ b/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), + ) 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..00a83ad --- /dev/null +++ b/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), + ) 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..03d3cfb --- /dev/null +++ b/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 + ) 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..6b3aafb --- /dev/null +++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SM_aerosol_anthro.py @@ -0,0 +1,77 @@ +from argparse import ArgumentParser +from pathlib import Path + +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_ancil_argparse import common_parser +from cmip7_ancil_common import ( + cmip7_date_constraint_from_years, + extend_years, + fix_coords, + interpolate_monthly, +) +from cmip7_SM import CMIP7_SM_BEG_YEAR, CMIP7_SM_END_YEAR + + +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") + 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_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_anthro(args, species): + cube = load_cmip7_aerosol( + args, + cmip7_sm_aerosol_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 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), + ) diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_SO2_interpolate.py index 993ab20..044b734 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"]] @@ -120,5 +119,8 @@ def save_cmip7_so2_aerosol_anthro( (dms_lon, 2), ], ) - - save_ancil([so2_low, so2_high, dms], save_dirpath, args.save_filename) + save_ancil( + [so2_low, so2_high, dms], + save_dirpath, + args.save_filename, + ) 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): 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): diff --git a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py index 5425e1a..09cc010 100644 --- a/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py +++ b/CMIP7/esm1p6/atmosphere/aerosol/cmip7_aerosol_common.py @@ -32,5 +32,7 @@ def zero_poles(cube): # For aerosol emissions they should be zero latdim = cube.coord_dims("latitude") assert latdim == (1,) - 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 diff --git a/CMIP7/esm1p6/atmosphere/cmip7_SM.py b/CMIP7/esm1p6/atmosphere/cmip7_SM.py new file mode 100644 index 0000000..060215c --- /dev/null +++ b/CMIP7/esm1p6/atmosphere/cmip7_SM.py @@ -0,0 +1,18 @@ +from pathlib import Path + +from cmip7_ancil_constants import ANCIL_TODAY + +CMIP7_SM_BEG_YEAR = 2022 +CMIP7_SM_END_YEAR = 2100 + + +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/cmip7_ancil_common.py b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py index f76bae1..e1317d9 100644 --- a/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py +++ b/CMIP7/esm1p6/atmosphere/cmip7_ancil_common.py @@ -8,6 +8,8 @@ import cf_units import cftime import iris +import iris.analysis +import iris.coord_categorisation import mule import numpy as np from cmip7_ancil_constants import ( @@ -135,6 +137,111 @@ 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") + 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 = 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, + ) + + # Interpolate each month separately and interleave the data + new_cube = _interpolate_months_separately(cube, tpoints) + + # Create new DimCoord with contiguous bounds + new_time_coord = iris.coords.DimCoord( + tpoints, + 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") + 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.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]) + 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) + 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 @@ -151,15 +258,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,) - cube.data[:, 0] = 0.0 - cube.data[:, -1] = 0.0 - - def save_ancil( cubes, save_dirpath, save_filename, gregorian=True, replace_bounds=False ): diff --git a/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py b/CMIP7/esm1p6/atmosphere/ghg/cmip7_HI_ghg_generate.py index d1fb43f..ab8f83b 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 ( + cmip7_ghg_update_namelists_file, + load_cmip7_ghg_series_mmr, ) @@ -33,176 +24,16 @@ 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..e2c6bac --- /dev/null +++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_SM_ghg_generate.py @@ -0,0 +1,41 @@ +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 ( + cmip7_ghg_update_namelists_file, + load_cmip7_ghg_series_mmr, +) + + +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() + + 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_GHG_BEG_YEAR, CMIP7_SM_END_YEAR + ) + + # Patch the greenhouse gas namelist. + cmip7_ghg_update_namelists_file( + ghg_mmr_dict, CMIP7_SM_GHG_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..e3bfe23 --- /dev/null +++ b/CMIP7/esm1p6/atmosphere/ghg/cmip7_ghg_series.py @@ -0,0 +1,178 @@ +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) 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_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) 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" 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) 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..d518556 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, "CMIP", "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..3ba7903 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,50 @@ 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) + 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) + 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_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) 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..ae6e8c9 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,149 @@ 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( + 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 + )