Commit 577bf02e authored by Antoine Berchet's avatar Antoine Berchet
Browse files

Implementing background in FLEXPART, including vertical and horizontal...

Implementing background in FLEXPART, including vertical and horizontal interpolation of global concentration fields
parent 7974424c
import os
import datetime
from pycif.utils import path
from logging import info
......@@ -31,6 +32,7 @@ def default_fetch(
dir_dd = dd.strftime(input_dir)
file_dd = dd.strftime(input_file)
tmp_files.append("{}/{}".format(dir_dd, file_dd))
tmp_dates.append(datetime.datetime.strptime(file_dd, input_file))
# Fetching
local_files = []
......
......@@ -71,12 +71,12 @@ def init_components(plugin):
else None)
)
if cmp is None and not hasattr(tracer, "read"):
if cmp is None and not hasattr(tracer, "fetch"):
warning(
"{} in your Yaml is not recognized as"
" a valid input for the model"
" or you did not give enough"
" info to define a read method"
" info to define a fetch method"
.format(comp)
)
tracer.fetch = default_fetch
......
......@@ -92,9 +92,9 @@ def create_domain(domain,
np.ravel_multi_index(in_indexes, zlat_glob.shape).flatten()
domain.raveled_indexes_glob = raveled_indexes
lat_out = np.delete(zlat_glob, raveled_indexes)
lon_out = np.delete(zlon_glob, raveled_indexes)
lat_out = np.delete(zlat_glob, [])
lon_out = np.delete(zlon_glob, [])
# Putting together longitudes and latitudes
domain.zlon = np.append(domain.zlon_in.flatten(), lon_out)[:, np.newaxis]
domain.zlat = np.append(domain.zlat_in.flatten(), lat_out)[:, np.newaxis]
......
from .read import read
from .get_domain import get_domain
from .write import write
import datetime
import glob
import os
import xarray as xr
from pycif.utils.classes.setup import Setup
def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
# Looking for a reference file to read lon/lat in
list_file = glob.glob("{}/*nc".format(ref_dir))
domain_file = None
# Either a file is specified in the Yaml
if ref_file in list_file:
domain_file = "{}/{}".format(ref_dir, ref_file)
# Or loop over available file regarding file pattern
else:
for flx_file in list_file:
try:
date = datetime.datetime.strptime(
os.path.basename(flx_file), ref_file
)
domain_file = flx_file
break
except ValueError:
continue
if domain_file is None:
raise Exception(
"EDGARv5 domain could not be initialized as no file was found"
)
# Read lon/lat in
nc = xr.open_dataset(domain_file, decode_times=False)
lon = nc["longitude"]
lat = nc["latitude"]
lon_min = lon.min()
lon_max = lon.max() + (lon[-1] - lon[-2])
lat_min = lat.min() / 2
lat_max = lat.max() + (lat[-1] - lat[-2])
nlon = lon.size
nlat = lat.size
# Initializes domain
setup = Setup.from_dict(
{
"domain": {
"plugin": {
"name": "dummy",
"version": "std",
"type": "domain",
},
"xmin": lon_min,
"xmax": lon_max,
"ymin": lat_min,
"ymax": lat_max,
"nlon": nlon,
"nlat": nlat,
}
}
)
Setup.load_setup(setup, level=1)
# Initializes vertical resolution
pres = nc["pressure"]
setup.domain.nlev = len(pres)
setup.domain.sigma_a = pres.values / 1e2
setup.domain.sigma_b = 0. * pres.values
return setup.domain
import datetime
import os
import numpy as np
import xarray as xr
from netCDF4 import Dataset
from pycif.utils.netcdf import readnc
def read(
self,
name,
tracdir,
tracfile,
varnames,
dates,
interpol_flx=False,
comp_type=None,
tracer=None,
model=None,
**kwargs
):
"""Get fluxes from pre-computed fluxes and load them into a pyCIF
variables
Args:
self: the fluxes Plugin
name: the name of the component
tracdir, tracfile: flux directory and file format
dates: list of dates to extract
interpol_flx (bool): if True, interpolates fluxes at time t from
values of surrounding available files
"""
list_files = tracfile
if type(tracfile) != list:
list_files = [dd.strftime(tracfile) for dd in dates]
data = []
for dd, ff in zip(dates, list_files):
tmp = readnc(os.path.join(tracdir, ff),
[self.varname_init])
data.append(tmp[0])
# Putting in DataArray
xmod = xr.DataArray(
data, coords={"time": dates}, dims=("time", "lev", "lat", "lon")
)
# Converting ppb to ppm for pycif
xmod /= 1e3
return xmod
import os
import copy
import numpy as np
import pandas as pd
from netCDF4 import Dataset
import xarray as xr
from pycif.utils.netcdf import save_nc
def write(self, name, lbc_file, data, mode="a", comp_type=None):
"""Write flux to INICOND or BOUN_CONC CHIMERE compatible files.
"""
# If mode is 'a' but file does not exit, switch to mode 'w'
if mode == "a" and not os.path.isfile(lbc_file):
mode = "w"
# Loading longitudes and latitudes
lon = self.domain.zlon
lat = self.domain.zlat
nlev = self.domain.nlev
# Write INI_CONCS
if comp_type == "inicond":
write_iniconcs(name, lbc_file, data, lon, lat, mode)
else:
lon_side = self.domain.zlon_side
lat_side = self.domain.zlat_side
write_bounconcs(
name,
lbc_file,
data,
lon,
lat,
lon_side,
lat_side,
nlev,
mode,
comp_type,
)
def write_iniconcs(name, lbc_file, data, lon, lat, mode="a"):
nhours, nlev, nmerid, nzonal = np.shape(data)
# Dimensions
spstrlen = 23
dimnames = [
"Time",
"south_north",
"west_east",
"bottom_top",
"SpStrLen",
"DateStrLen",
"Species",
]
dimlens = [None, nmerid, nzonal, nlev, spstrlen, 19, 1]
# Variables names, dimension and attributes
varnames = ["Times", "species", "lon", "lat", name]
vardims = [
("Time", "DateStrLen"),
("Species", "SpStrLen"),
("south_north", "west_east"),
("south_north", "west_east"),
("bottom_top", "south_north", "west_east"),
]
dtypes = ["c", "c", "f", "f", "d"]
units = ["", "", "degrees_east", "degrees_north", "ppb"]
attributes = [
{},
{},
{"long_name": "Longitude"},
{"long_name": "Latitude"},
{"long_name": "{}".format(name)},
]
# Variables to save
times = [
list(pd.to_datetime(d).strftime("%Y-%m-%d_%H:00:00"))
for d in data.time.values
]
specs_var = [list(name.ljust(spstrlen))]
# Append species to existing file
if os.path.isfile(lbc_file) and mode == "a":
with Dataset(lbc_file, "a") as f:
ljust_specs_in = f.variables["species"][:].astype(str)
specs_in = ["".join(p).strip() for p in ljust_specs_in]
if name not in specs_in:
specs = [list(name.ljust(spstrlen))]
specs = xr.Dataset(
{"species": (("Species", "SpStrLen"), specs)}
)
ds = xr.open_dataset(lbc_file)
ds_spec = xr.Dataset(
{"species": (("Species", "SpStrLen"), ljust_specs_in)}
)
del ds["species"]
ds_spec = xr.concat([ds_spec["species"], specs["species"]], "Species")
ds["species"] = ds_spec
ds.to_netcdf(lbc_file, format="NETCDF3_CLASSIC")
specs_in += [name]
specs_var = [list(s.ljust(spstrlen)) for s in specs_in]
variables = [times, specs_var, lon, lat, data[0].values]
save_nc(
lbc_file,
variables,
varnames,
vardims,
dimnames,
dimlens,
units,
dtypes,
mode=mode,
attributes=attributes,
format="NETCDF3_CLASSIC",
)
def write_bounconcs(
name,
lbc_file,
data,
lon,
lat,
lon_side,
lat_side,
nlev,
mode="a",
comp_type=None,
):
nhours = len(data)
nsides = len(lon_side)
nmerid, nzonal = np.shape(lon)
# Dimensions
spstrlen = 23
datestrlen = 19
dimnames = [
"Time",
"south_north",
"west_east",
"bottom_top",
"h_boundary",
"SpStrLen",
"DateStrLen",
"Species",
]
dimlens = [None, nmerid, nzonal, nlev, nsides, spstrlen, datestrlen, 1]
# Variables names, dimension and attributes
varnames = ["Times", "species", "lon", "lat", "top_conc", "lat_conc"]
vardims = [
("Time", "DateStrLen"),
("Species", "SpStrLen"),
("south_north", "west_east"),
("south_north", "west_east"),
("Time", "south_north", "west_east", "Species"),
("Time", "bottom_top", "h_boundary", "Species"),
]
dtypes = ["c", "c", "f", "f", "d", "d"]
units = ["", "", "degrees_east", "degrees_north", "ppb", "ppb"]
attributes = [
{},
{},
{"long_name": "Longitude"},
{"long_name": "Latitude"},
{"long_name": "{}".format(varnames[-2])},
{"long_name": "{}".format(varnames[-1])},
]
# Variables to save
times = [
list(pd.to_datetime(d).strftime("%Y-%m-%d_%H:00:00"))
for d in data.time.values
]
# Output data
top_in = (
data[:, 0]
if comp_type == "topcond"
else np.zeros((nhours, nmerid, nzonal, 1))
)
lat_in = (
data.data[..., 0, :, np.newaxis]
if comp_type == "latcond"
else np.zeros((nhours, nlev, nsides, 1))
)
ljust_specs_in = [list(name.ljust(spstrlen))]
specs_in = [name]
# Append species to existing file
if os.path.isfile(lbc_file) and mode == "a":
with Dataset(lbc_file, "a") as f:
ljust_specs_in = f.variables["species"][:].astype(str)
specs_in = ["".join(p).strip() for p in ljust_specs_in]
top_in = f.variables["top_conc"][:]
lat_in = f.variables["lat_conc"][:]
if name in specs_in:
ispec = specs_in.index(name)
if comp_type == "topcond":
top_in[..., ispec] = data[:, 0]
else:
lat_in[..., ispec] = data[..., 0, :]
else:
specs_in = specs_in + [name]
ljust_specs_in = [list(s.ljust(spstrlen)) for s in specs_in]
if comp_type == "topcond":
top_in = np.concatenate(top_in, data, axis=3)
lat_in = np.concatenate(
lat_in, np.zeros((nhours, nlev, nsides, 1)), axis=3
)
else:
top_in = np.concatenate(
top_in, np.zeros((nhours, nmerid, nzonal, 1)), axis=3
)
lat_in = np.concatenate(lat_in, data, axis=3)
variables = [times, ljust_specs_in, lon, lat, top_in, lat_in]
save_nc(
lbc_file,
variables,
varnames,
vardims,
dimnames,
dimlens,
units,
dtypes,
mode=mode,
attributes=attributes,
format="NETCDF3_CLASSIC",
)
......@@ -6,8 +6,10 @@ from . import chimere_icbc
from . import grib2_ecmwf
from . import lmdz_ic
from . import lmdz_prescrconcs
from . import noaa_glob_avg
Fields.register_plugin("CHIMERE", "icbc", chimere_icbc)
Fields.register_plugin("LMDZ", "ic", lmdz_ic)
Fields.register_plugin("LMDZ", "prescrconcs", lmdz_prescrconcs)
Fields.register_plugin("ECMWF", "grib2", grib2_ecmwf)
Fields.register_plugin("NOAA", "glob_avg", noaa_glob_avg)
......@@ -68,7 +68,7 @@ def read(self, name, tracdir, tracfile, varnames, dates,
# Extract data outside nest domain
flx_reg_out = \
np.delete(data.reshape(data.shape[0], -1),
self.domain.raveled_indexes_glob, 1)
[], 1)
else:
flx_reg_out = \
......@@ -80,7 +80,7 @@ def read(self, name, tracdir, tracfile, varnames, dates,
np.append(flx_reg_in, flx_reg_out, axis=1)[:, :, np.newaxis],
axis=0)
# Put data into dataarry
# Put data into dataarray
xmod = xr.DataArray(trcr_flx[:, np.newaxis],
coords={'time': np.array(times)},
dims=('time', 'lev', 'lat', 'lon'))
......
import os
from .ini_periods import ini_periods
from .run import run
from .ini_mapper import ini_mapper
......@@ -5,6 +7,7 @@ from .io.native2inputs import native2inputs
from .io.native2inputs_adj import native2inputs_adj
from .io.outputs2native import outputs2native
from .io.outputs2native_adj import outputs2native_adj
from .utils.flexpart_header import read_header
requirements = {'domain': {'name': 'FLEXPART', 'version': 'std',
'empty': False, 'any': False},
......@@ -15,7 +18,8 @@ requirements = {'domain': {'name': 'FLEXPART', 'version': 'std',
required_inputs = ['fluxes']
default_values = {
"read_background": False
"read_background": False,
"periods": "1MS"
}
......@@ -31,12 +35,17 @@ def ini_data(plugin, **kwargs):
"""
# Initializes default values
# Period of sub-simulations: default = 1 month
if not hasattr(plugin.plugin, 'periods'):
plugin.periods = '1MS'
else:
plugin.periods = plugin.plugin.periods
# Update required inputs if reading background
if plugin.read_background:
plugin.required_inputs.append("inicond")
# Read header for vertical resolution
fp_header_nest = read_header(
os.path.join(plugin.run_dir_nest,
plugin.outheight_header))
plugin.domain.heights = \
fp_header_nest.outheight[fp_header_nest.outheight > 0]
plugin.domain.nlev = len(plugin.domain.heights)
return plugin
......@@ -11,11 +11,14 @@ def ini_mapper(model, transform_type, inputs={}, outputs={}, backup_comps={}):
# Executable
mapper = {'inputs':
{('fluxes', s): dict_surface
for s in ['CH4']},
{('fluxes', s): dict_surface for s in ['CH4']},
'outputs': {('concs', s): {}
for s in ['CH4']}
}
# Updating mapper if using background
if model.read_background:
mapper["inputs"].update({('inicond', s): dict_surface for s in ['CH4']})
return mapper
import os
from logging import debug
import numpy as np
import pandas as pd
from .read import read_flexpart_grid
def flux_contribution(self, mode, subdir, dataobs,
fp_header_nest, fp_header_glob):
for obs_i, row in enumerate(dataobs.itertuples()):
row = dataobs.iloc[obs_i]
station = row.station
fluxes = {"sim": list(self.dataflx.values())[0]["spec"]}
if mode == "tl" and "incr" in list(self.dataflx.values())[0]:
fluxes["sim_tl"] = list(self.dataflx.values())[0]["incr"]
# Read nested grids
runsubdir_nest = os.path.join(
self.run_dir_nest, station.upper(), subdir)
file_date = row.date.strftime('%Y%m%d%H%M%S')
file_name = 'grid_time_nest_{}_001'.format(file_date)
debug("Reading {} for station {}".format(file_name, station))
if not os.path.isfile(os.path.join(runsubdir_nest, file_name)):
return
grid_nest, gtime, ngrid = read_flexpart_grid(
runsubdir_nest, file_name, fp_header_nest,
numscale=self.numscale)
# Conversion of footprints
grid_nest *= self.coeff * self.mmair / self.molarmass
# Multiply footprints with fluxes for fwd and tl
for data_id in fluxes:
dataflx = fluxes[data_id]
flx_dates = pd.DatetimeIndex(
dataflx.time.values).to_pydatetime()
inds_flx = \
np.argmin(np.abs(np.array(gtime)[:, np.newaxis]
- flx_dates[np.newaxis, :]), axis=1)
nest_sim = (
grid_nest.T[:ngrid].reshape(ngrid, -1)
* dataflx[inds_flx, 0, :self.domain.zlon_in.size, 0]) \
.sum()
# Filling simulation
sim_col = dataobs.columns.get_indexer([data_id])
dataobs.iloc[obs_i, sim_col] = nest_sim
# Read global footprints
# TODO: read correction factor dry air!
runsubdir_glob = os.path.join(
self.run_dir_glob, station.upper(), subdir)
file_name = 'grid_time_{}_001'.format(file_date)
if not os.path.isfile(os.path.join(runsubdir_glob, file_name)):
return
grid_glob, gtime_glob, ngrid_glob = \
read_flexpart_grid(
runsubdir_glob, file_name, fp_header_glob,
numscale=self.numscale)