Skip to content

Instantly share code, notes, and snippets.

@atteggiani
Last active November 30, 2023 08:33
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save atteggiani/797585a6ff215a590174dadd3526bcef to your computer and use it in GitHub Desktop.
Example of how to create the Nitrogen deposition ancillary file for ACCESS-ESM1.5 from a set of netCDF source files.
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment