Example of how to create the Nitrogen deposition ancillary file for ACCESS-ESM1.5 from a set of netCDF source files.

make_nitrogen_deposition_ancil.py

#!/usr/bin/env python3
"""
Create the ACCESS-ESM1.5 ancillary file needed for nitrogen deposition in UKESM1 for the CMIP6 SSP370 scenarioMIP
00447    0 NITROGEN DEPOSITION (kgN/m2/s) (ACCESS-CM2)
or 
00884    0 NITROGEN DEPOSITION (gN/m2/d) (ACCESS-ESM1.5)
"""
import os
import ants
from ants.fileformats.ancil import _cubes_to_ancilfile, _mule_set_lbuser2
import iris
import xarray as xr
import mule
import iris.fileformats
import iris.analysis
import numpy as np

ancdir = "/scratch/public/dm5220/tilo"

# Read in the resolution-specific land-sea mask
gridFile_esm = "/g/data/access/projects/access/data/ancil/access_v2/qrparm.mask"
# gridFile_cm = "/g/data/access/projects/access/data/ancil/access_cm2_n96e/O1/qrparm.mask"

# Source variable(s)
sourcevar = ["drynhx","wetnhx","drynoy","wetnoy"]
# Path to source file(s)
infile = []
for var in sourcevar:
    infile.append(os.path.join(ancdir,"source",f"{var}.nc"))

# Path to the ancillary historical file
hist_file = "/g/data/p66/txz599/data/ancil/CMIP6/Ndep_1849_2015.anc"

# Specify output ancillary file(s) to produce
ancil_name = "ndep_2014_2101_ssp370_esm1.5.anc"
stash = 'm01s00i884'
description = "Nitrogen deposition field for driving UKESM1 (ssp370)"
outfile = os.path.join(ancdir,ancil_name)

def regrid(inCube,gridfile):
    """
    Subroutine to regrid source data from its native grid to
    the model grid at the required resolution.
    """
    gridCube = ants.load_cube(gridfile,iris.AttributeConstraint(STASH='m01s00i030'))
    inCube.coord('latitude').coord_system=gridCube.coord('latitude').coord_system
    inCube.coord('longitude').coord_system=gridCube.coord('longitude').coord_system
    inCube.coord('longitude').circular=True
    gridCube.coord('longitude').circular=True
    ants.utils.cube.guess_horizontal_bounds(inCube)
    ants.utils.cube.guess_horizontal_bounds(gridCube)
    outCube = inCube.regrid(gridCube, iris.analysis.AreaWeighted())
    return outCube

def convert_units(data):
    '''
    Convert units from kgN/m2/s to gN/m2/d
    '''
    coeff = 1e3*60*60*24 
    return data*coeff

def preprocess():
    """
    Preprocess source file(s) to create a single file with
    the sum of all source variables, and duplicate last year of data
    to cover until end of 2100
    """
    xr.set_options(keep_attrs=True)
    # Set up the data as sum of all source variables
    data = xr.open_mfdataset(infile)
    sum_data = data[sourcevar].to_array().sum('variable')
    sum_data.name = 'ndep'
    # Duplicate last year twice to cover until end of 2101
    end_data = xr.apply_ufunc(lambda x: np.tile(x,(2,1,1)),
        sum_data[-12:],
        dask='allowed',
        exclude_dims=set(['time']),
        input_core_dims=[['time','lat','lon']],
        output_core_dims=[['time','lat','lon']]
    )
    # Create new time coord for last 2 years, by adding 2 to
    # last 2 years of historical data
    new_time_coord = [dt.replace(year=dt.year+2) for dt in sum_data.time[-12*2:].values]
    # set new time coord as time coordinate for last_year 
    new = xr.concat([sum_data,end_data.assign_coords(time=new_time_coord)],dim='time')
    return new.to_iris()

    
def main():
    '''
    Create ancillary file
    '''
    cube = preprocess()
    cube.var_name = "n_dep"
    cube.long_name = "NITROGEN DEPOSITION (gN/m2/d)"
    cube.attributes["STASH"]=iris.fileformats.pp.STASH.from_msi(stash)
    cube.attributes["grid_staggering"] = 3 #New Dynamics
    cube = regrid(cube,gridFile_esm)
    cube.data = convert_units(np.clip(cube.data,0.0,float('Inf')))
    cubes = ants.utils.cube.as_cubelist(cube)
    ancilfile = _cubes_to_ancilfile(cubes)
    # Get year '2014' of historical data and prepend it to ancilfile to cover from 2014
    hdata_2014 = list(filter(lambda x: x.lbyr == 2014, mule.load_umfile(hist_file).fields))
    for field in hdata_2014:
        field.set_data_provider(mule.ArrayDataProvider(field.get_data().astype(np.float64)))
    ancilfile.fields = hdata_2014+ancilfile.fields
    # Set lbuser2
    _mule_set_lbuser2(ancilfile)
    # Set some other header values
    ancilfile.fixed_length_header.raw[8] = 1 # proleptic gregorian calendar
    ancilfile.fixed_length_header.raw[12] = 703 # UM model version
    ancilfile.fixed_length_header.raw[21] = 2014 # T1 year
    ancilfile.fixed_length_header.raw[22] = 1 # T1 month
    ancilfile.fixed_length_header.raw[23] = 16 # T1 day
    ancilfile.fixed_length_header.raw[24] = 0 # T1 hour
    ancilfile.fixed_length_header.raw[28] = 2101 # T2 year
    ancilfile.fixed_length_header.raw[29] = 12 # T2 month
    ancilfile.fixed_length_header.raw[30] = 16 # T2 day
    ancilfile.fixed_length_header.raw[31] = 0 # T2 hour
    ancilfile.integer_constants.raw[3] = cube.coord('time').shape[0]+12 # number of timesteps
    for field in ancilfile.fields:
        field.lbhr = 0 # Hour
        field.raw[7:13] = 0 # LBYRD to LBDAYD 
        field.lbtim = 1 # proleptic gregorian calendar
        field.lbpack = 0 # Not packed
        field.lbrel = 3 # Header release for UM > 8.1
        field.lbsrce = 7031111 # UM version 703
        field.lbvc = 129 # Surface level type
        field.lblev = 9999 # Surface
        field.lbuser3 = 0 # Rimwidth and Halo size
        field.lbuser6 = 0
        field.lbproc = 0
        field.bmks = 1.0
        field.lbnrec = 2048-field.lblrec%2048 + field.lblrec
    # Write the ancillary file
    ancilfile.to_file(outfile)
    print(f"Written {outfile}")

if __name__ == '__main__':
    main()