#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Author: Dominik Brunner (DB), Empa, Switzerland
# DB, 15 May 2017: first implementation on daint.cscs.ch
# DB, 06 Jun 2017: modified to add surface pressure fields
# hjm, 22 Jun 2018: Translated to Python
######################
## Problems:
## - when comparing to the "processed2" icbc form Gerrit, PSURF is slightly different
## - when checking the website, the b levels are not the same as in CAMS
## - It is not possible to delete a variable from a netcdf4 file, so resort to calling ncks at the end
import argparse
import os
import logging
import shutil
import netCDF4 as nc
import numpy as np
from subprocess import call
import sys
import importlib
###############################
## Constants : CAMS levels ##
###############################
# hybrid coefficients of the 137 model level version
# http://www.ecmwf.int/en/forecasts/documentation-and-support/137-model-levels
a_half_137 = [
0.0000000000e+00, 2.0003650188e+00, 3.1022410393e+00, 4.6660838127e+00,
6.8279771805e+00, 9.7469663620e+00, 1.3605423927e+01, 1.8608930588e+01,
2.4985717773e+01, 3.2985710144e+01, 4.2879241943e+01, 5.4955463409e+01,
6.9520576477e+01, 8.6895881653e+01, 1.0741574097e+02, 1.3142550659e+02,
1.5927940369e+02, 1.9133856201e+02, 2.2796894836e+02, 2.6953958130e+02,
3.1642074585e+02, 3.6898236084e+02, 4.2759249878e+02, 4.9261602783e+02,
5.6441345215e+02, 6.4333990479e+02, 7.2974414062e+02, 8.2396783447e+02,
9.2634490967e+02, 1.0372011719e+03, 1.1568536377e+03, 1.2856103516e+03,
1.4237701416e+03, 1.5716229248e+03, 1.7294489746e+03, 1.8975192871e+03,
2.0760959473e+03, 2.2654316406e+03, 2.4657705078e+03, 2.6773481445e+03,
2.9003913574e+03, 3.1351193848e+03, 3.3817436523e+03, 3.6404682617e+03,
3.9114904785e+03, 4.1949306641e+03, 4.4908173828e+03, 4.7991494141e+03,
5.1198950195e+03, 5.4529907227e+03, 5.7983447266e+03, 6.1560742188e+03,
6.5269467773e+03, 6.9118706055e+03, 7.3118691406e+03, 7.7274121094e+03,
8.1593540039e+03, 8.6085253906e+03, 9.0764003906e+03, 9.5626826172e+03,
1.0065978516e+04, 1.0584631836e+04, 1.1116662109e+04, 1.1660067383e+04,
1.2211547852e+04, 1.2766873047e+04, 1.3324668945e+04, 1.3881331055e+04,
1.4432139648e+04, 1.4975615234e+04, 1.5508256836e+04, 1.6026115234e+04,
1.6527322266e+04, 1.7008789062e+04, 1.7467613281e+04, 1.7901621094e+04,
1.8308433594e+04, 1.8685718750e+04, 1.9031289062e+04, 1.9343511719e+04,
1.9620042969e+04, 1.9859390625e+04, 2.0059931641e+04, 2.0219664062e+04,
2.0337863281e+04, 2.0412308594e+04, 2.0442078125e+04, 2.0425718750e+04,
2.0361816406e+04, 2.0249511719e+04, 2.0087085938e+04, 1.9874025391e+04,
1.9608572266e+04, 1.9290226562e+04, 1.8917460938e+04, 1.8489707031e+04,
1.8006925781e+04, 1.7471839844e+04, 1.6888687500e+04, 1.6262046875e+04,
1.5596695312e+04, 1.4898453125e+04, 1.4173324219e+04, 1.3427769531e+04,
1.2668257812e+04, 1.1901339844e+04, 1.1133304688e+04, 1.0370175781e+04,
9.6175156250e+03, 8.8804531250e+03, 8.1633750000e+03, 7.4703437500e+03,
6.8044218750e+03, 6.1685312500e+03, 5.5643828125e+03, 4.9937968750e+03,
4.4573750000e+03, 3.9559609375e+03, 3.4892343750e+03, 3.0572656250e+03,
2.6591406250e+03, 2.2942421875e+03, 1.9615000000e+03, 1.6594765625e+03,
1.3875468750e+03, 1.1432500000e+03, 9.2650781250e+02, 7.3499218750e+02,
5.6806250000e+02, 4.2441406250e+02, 3.0247656250e+02, 2.0248437500e+02,
1.2210156250e+02, 6.2781250000e+01, 2.2835937500e+01, 3.7578129768e+00,
0.0000000000e+00, 0.0000000000e+00
]
b_half_137 = [
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 3.8199999608e-08, 6.7607002165e-06,
2.4348000807e-05, 5.8921999880e-05, 1.1191429803e-04, 1.9857739971e-04,
3.4037968726e-04, 5.6155532366e-04, 8.8969792705e-04, 1.3528055279e-03,
1.9918379840e-03, 2.8571242001e-03, 3.9709536359e-03, 5.3778146394e-03,
7.1333767846e-03, 9.2614600435e-03, 1.1806022376e-02, 1.4815628529e-02,
1.8318451941e-02, 2.2354844958e-02, 2.6963520795e-02, 3.2176095992e-02,
3.8026399910e-02, 4.4547960162e-02, 5.1773015410e-02, 5.9728413820e-02,
6.8448252976e-02, 7.7958308160e-02, 8.8285736740e-02, 9.9461667240e-02,
1.1150465161e-01, 1.2444812804e-01, 1.3831289113e-01, 1.5312503278e-01,
1.6891041398e-01, 1.8568944931e-01, 2.0349121094e-01, 2.2233286500e-01,
2.4224400520e-01, 2.6324188709e-01, 2.8535401821e-01, 3.0859845877e-01,
3.3293908834e-01, 3.5825419426e-01, 3.8436332345e-01, 4.1112476587e-01,
4.3839120865e-01, 4.6600329876e-01, 4.9380031228e-01, 5.2161920071e-01,
5.4930114746e-01, 5.7669216394e-01, 6.0364806652e-01, 6.3003581762e-01,
6.5573596954e-01, 6.8064302206e-01, 7.0466899872e-01, 7.2773873806e-01,
7.4979656935e-01, 7.7079755068e-01, 7.9071676731e-01, 8.0953603983e-01,
8.2725608349e-01, 8.4388113022e-01, 8.5943180323e-01, 8.7392926216e-01,
8.8740754128e-01, 8.9990049601e-01, 9.1144818068e-01, 9.2209565639e-01,
9.3188077211e-01, 9.4085955620e-01, 9.4906443357e-01, 9.5654952526e-01,
9.6335172653e-01, 9.6951341629e-01, 9.7507840395e-01, 9.8007160425e-01,
9.8454189301e-01, 9.8849952221e-01, 9.9198400974e-01, 9.9500250816e-01,
9.9763011932e-01, 1.0000000000e+00
]
# hybrid coefficients of the 60 model level version
# http://www.ecmwf.int/en/forecasts/documentation-and-support/60-model-levels
a_half_60 = [
0.0000000000e+00, 2.0000000000e+01, 3.8425338745e+01, 6.3647796631e+01,
9.5636962891e+01, 1.3448330688e+02, 1.8058435059e+02, 2.3477905273e+02,
2.9849584961e+02, 3.7397192383e+02, 4.6461816406e+02, 5.7565112305e+02,
7.1321801758e+02, 8.8366040039e+02, 1.0948347168e+03, 1.3564746094e+03,
1.6806403809e+03, 2.0822739258e+03, 2.5798886719e+03, 3.1964216309e+03,
3.9602915039e+03, 4.9067070312e+03, 6.0180195312e+03, 7.3066328125e+03,
8.7650546875e+03, 1.0376125000e+04, 1.2077445312e+04, 1.3775324219e+04,
1.5379804688e+04, 1.6819472656e+04, 1.8045183594e+04, 1.9027695312e+04,
1.9755109375e+04, 2.0222203125e+04, 2.0429863281e+04, 2.0384480469e+04,
2.0097402344e+04, 1.9584328125e+04, 1.8864750000e+04, 1.7961359375e+04,
1.6899468750e+04, 1.5706449219e+04, 1.4411125000e+04, 1.3043218750e+04,
1.1632757812e+04, 1.0209500000e+04, 8.8023554688e+03, 7.4388046875e+03,
6.1443164062e+03, 4.9417773438e+03, 3.8509133301e+03, 2.8876965332e+03,
2.0637797852e+03, 1.3859125977e+03, 8.5536181641e+02, 4.6733349609e+02,
2.1039389038e+02, 6.5889236450e+01, 7.3677425385e+00, 0.0000000000e+00,
0.0000000000e+00
]
b_half_60 = [
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00,
7.5823496445e-05, 4.6139489859e-04, 1.8151560798e-03, 5.0811171532e-03,
1.1142909527e-02, 2.0677875727e-02, 3.4121163189e-02, 5.1690407097e-02,
7.3533833027e-02, 9.9674701691e-02, 1.3002252579e-01, 1.6438430548e-01,
2.0247590542e-01, 2.4393314123e-01, 2.8832298517e-01, 3.3515489101e-01,
3.8389211893e-01, 4.3396294117e-01, 4.8477154970e-01, 5.3570991755e-01,
5.8616840839e-01, 6.3554745913e-01, 6.8326860666e-01, 7.2878581285e-01,
7.7159661055e-01, 8.1125342846e-01, 8.4737491608e-01, 8.7965691090e-01,
9.0788388252e-01, 9.3194031715e-01, 9.5182150602e-01, 9.6764522791e-01,
9.7966271639e-01, 9.8827010393e-01, 9.9401944876e-01, 9.9763011932e-01,
1.0000000000e+00
]
[docs]def main(date, inpath, outpath, param):
"""Prepare CAMS CO2, CO and NOx boundary conditions for
**int2lm/int2cosmo** for the project SMARTCARB.
The CAMS data consists of
- CO2 and CO fields of experiment gf39, class RD (approx. 15 km
resolution, 137 levels). See
https://atmosphere.copernicus.eu/change-log-gf39
- NO and NO2 from the CAMS operational product, exp 0001, class MC
(60 km resolution)
The data sets are retrieved as individual 3-hourly files from the MARS
archive at ECMWF with names::
cams_0001_2015010500.nc # for NO and NO2
sfc_0001_2015010500.nc # for log of surface pressure
cams_gf39_2015010500.nc # for CO and CO2
sfc_gf39_2015010500.nc # for log of surface pressure
The path to the directory of the CAMS data can optionally be supplied,
otherwise the script should be invoked in the directory of the CAMS data.
The script generates 8 individual 3-hourly IC/BC files::
cams_nox_yyyymmddhh.nc
cams_co2_yyyymmddhh.nc
Usage::
python cams4int2cosmo.py date [-i inpath -o outpath]
Output:
``cams_NOX_YYYYMMDD00.nc`` to ``cams_NOX_YYYYMMDD21.nc``
``cams_CO2_YYYYMMDD00.nc`` to ``cams_CO2_YYYYMMDD21.nc``
Parameters
----------
date : str
date in format YYYYMMDD
inpath : str
path of original CAMS files (default is current path)
outpath : str
path where output files should be generated (default is current path)
param : dict
"""
try:
os.makedirs(outpath, exist_ok=True)
except (OSError, PermissionError):
logging.error("Creating output directory failed")
raise
tracer2dict = {
'CO2': dict(name="CO2_BG",
short_name="co2",
long_name="carbon_dioxide"),
'CO': dict(name="CO_BG", short_name="co", long_name="carbon_monoxide"),
'CH4': dict(name="CH4_BG", short_name="ch4", long_name="methane"),
'NOX': dict(name="NOX_BG",
short_name="NOx",
long_name="nitrogen_oxide"),
}
species = []
for s in param["species"]:
logging.info(s)
logging.info(tracer2dict[s])
try:
species.append(tracer2dict[s])
except KeyError:
logging.error("Variable " + s +
" is not part of the available list of variables.")
main_process(date, inpath, outpath, species, param)
def main_process(date, inpath, outpath, species, param):
if param["lev"] == 137:
a_half = a_half_137 #defined globally
b_half = b_half_137 #defined globally
offset = 37 # only levels 38 to 137 are retrieved from MARS
elif param["lev"] == 60:
a_half = a_half_60 #defined globally
b_half = b_half_60 #defined globally
offset = 0 # all levels
else:
logging.error(
"The list of hybrid parameters for the amount of levels " +
param["lev"] + "is not defined.")
raise
hyam = [(a_half[i] + a_half[i + 1]) / 2.
for i in range(offset,
len(a_half) - 1)]
hybm = [(b_half[i] + b_half[i + 1]) / 2.
for i in range(offset,
len(a_half) - 1)]
to_print = ",".join(
[species[i]["short_name"] for i in range(len(species))])
logging.info('Processing ' + to_print + ' for time ' +
date.strftime("%Y%m%d%H"))
infile = os.path.join(
inpath, param["prefix1"] + "_" + date.strftime("%Y%m%d%H") + ".nc")
sfcfile = os.path.join(
inpath, param["prefix2"] + "_" + date.strftime("%Y%m%d%H") + ".nc")
outfile = os.path.join(
outpath, param["suffix"] + "_" + date.strftime("%Y%m%d%H") + ".nc")
try:
shutil.copy(infile, outfile)
except FileNotFoundError:
logging.error("file " + infile + " not found")
raise
except (PermissionError, OSError):
logging.error("Copying file " + infile + " failed")
raise
with nc.Dataset(outfile, "a", format="NETCDF4") as outf:
# copy surface pressure from surface file
with nc.Dataset(sfcfile, "r") as inf:
lnsp = inf.variables["lnsp"]
outf.createVariable("PSURF", "f8", lnsp.dimensions)
outf.variables["PSURF"][:] = np.exp(lnsp[:])
outf["PSURF"].setncattr("long_name", "surface pressure")
outf["PSURF"].setncattr("units", "Pa")
# change the calender attribute to proleptic_gregorian
outf["time"].setncattr("calendar", "proleptic_gregorian")
# rename fields according to Carbosense4D definitions
for s in species:
if s["name"] != "NOX_BG":
outf.renameVariable(s["short_name"], s["name"])
else:
# Create NOX from NO and NO2
outf.createVariable("NOX_BG", "f8", outf["no"].dimensions)
outf["NOX_BG"][:] = outf["no"][:] + outf["no2"][:]
outf.renameDimension("latitude", "lat")
outf.renameDimension("longitude", "lon")
outf.renameVariable("latitude", "lat")
outf.renameVariable("longitude", "lon")
# Reverse latitude for all fields so that startlat < endlat
for v in outf.variables:
dims = outf[v].dimensions
if 'lat' in dims:
lat_ind = dims.index('lat')
if lat_ind == 0:
outf[v][:] = outf[v][::-1]
elif lat_ind == 1:
outf[v][:] = outf[v][:, ::-1]
elif lat_ind == 2:
outf[v][:] = outf[v][:, :, ::-1]
elif lat_ind == 3:
outf[v][:] = outf[v][:, :, :, ::-1]
elif lat_ind == 4:
outf[v][:] = outf[v][:, :, :, :, ::-1]
# Add the units; these must match those in $int2cosmo/src/trcr_gribtabs.f90
for s in species:
chem = s["name"]
long = s["long_name"]
outf[chem].setncattr("units", "kg kg-1")
outf[chem].setncattr("units_desc",
"kg of substance per kg of dry air")
outf[chem].setncattr(
"long_name",
"mass mixing ratio of " + chem[:-3] + " from outside domain")
outf[chem].setncattr("standard_name",
"mass_fraction_of_" + long + "_in_air")
# Add hybrid coefficients
outf.createVariable("hyam", np.float64, "level")
outf.createVariable("hybm", np.float64, "level")
outf["hyam"][:] = hyam
outf["hybm"][:] = hybm
outf["hyam"].setncattr("units", "Pa")
outf["hybm"].setncattr("units", "1")
outf["hyam"].setncattr("long_name",
"hybrid A coefficient at layer midpoints")
outf["hybm"].setncattr("long_name",
"hybrid B coefficient at layer midpoints")
# create new dimension lev1
outf.createDimension("lev1", 1)
outf.createVariable("P0", float, "lev1")
outf["P0"][:] = [1]
outf["P0"].setncattr("units", "Pa")
outf["P0"].setncattr("long_name", "reference pressure")
# if NOX_BG, delete no and no2
for s in species:
if s["name"] == "NOX_BG":
if "no" in outf.variables and "no2" in outf.variables:
call([
"ncks", "--overwrite", "-x", "-v", "no,no2", outfile,
outfile
])
if __name__ == "__main__":
parser = argparse.ArgumentParser(
description=
"Prepare CAMS CO2, CO, CH4 and NOx boundary conditions for int2lm/int2cosmo",
formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument('date', type=str, help='date in format YYYYMMDD')
parser.add_argument('param', type=str, help='dictionary of the parameters')
parser.add_argument(
'-i',
type=str,
metavar="inpath",
help='path of original CAMS files (default is current path)',
default=os.getcwd())
parser.add_argument(
'-o',
type=str,
metavar="outpath",
help=
"path where output files should be generated (default is inpath/processed)",
)
indir = parser.get_default("i")
parser.set_defaults(o=os.path.join(indir, "processed"))
args = parser.parse_args()
main(args.date, args.i, args.o, args.param)