import os
import shutil
import filecmp
import copy
from netCDF4 import Dataset
import numpy as np
import xarray as xr
import pandas as pd
import pathlib
from logging import warning, debug
from ......utils import path
from ......utils.hdf5 import _hdf5_lock
from .utils import replace_dates
[docs]
def make_boundcond(self, datastore, runsubdir, ddi, mode, input_type):
"""
Generates boundary conditions files for CHIMERE
:param self:
:param datastore:
:type datastore: dict
:param runsubdir:
:type runsubdir: str
:param sdc:
:type sdc: str
:param hour_dates:
:param mode:
:param input_type:
:return:
"""
# List of dates for which emissions are needed
list_dates = pd.date_range(ddi, periods=self.nhours + 1, freq="h")
# Name of variable in netCDF
nc_varname = "top_conc" if input_type == "topcond" else "lat_conc"
# Fixed name for BC files
fileout = f"{runsubdir}/BOUN_CONCS.nc"
fileoutincr = f"{runsubdir}/BOUN_CONCS.increment.nc"
# Loop on all active species
# If in datastore, take data, otherwise, link to original BOUN_CONCS
data2dump = {}
data2dump_tl = {}
for spec in self.chemistry.acspecies.attributes:
trid = (input_type, spec)
ref_trid = (input_type, spec)
# Update ref_trid if in ensemble mode
perturbed_from_reference = False
if hasattr(self, "perturbed_species"):
if spec in self.perturbed_species:
ref_trcr = self.perturbed_species[spec]
ref_trid = (input_type, ref_trcr)
# If spec not explicitly defined in datastore,
# fetch general component information if available
if trid in datastore:
pass
elif trid in [(k[0], k[1].replace("__sample#", "")) for k in datastore]:
trid_ind = [
(k[0], k[1].replace("__sample#", "")) for k in datastore
].index(trid)
trid = list(datastore.keys())[trid_ind]
# If trid comes from an unperturbed input, but in ensemble mode
# needs to fetch info from the reference species
elif ref_trid in datastore:
trid = ref_trid
perturbed_from_reference = True
# If spec not explicitly defined in datastore,
# fetch general component information if available
elif trid not in datastore and (input_type, "") in datastore:
trid = (input_type, "")
else:
continue
tracer = datastore[trid]
tracer_data = tracer["data"][ddi]
# If no data is provided, just copy from original file
if "spec" not in tracer_data:
dirorig = tracer["dirorig"]
fileorig = tracer["fileorig"]
fileini = ddi.strftime(f"{dirorig}/{fileorig}")
# If does not exist, just link
linked = False
if not os.path.isfile(fileout):
path.link(fileini, fileout)
linked = True
# Otherwise, check for difference
if not linked or perturbed_from_reference:
if not os.path.isfile(fileini):
warning(
f"The boundary condition file {fileini} does not exist to initialize the species {spec}. The concentrations will be initialized to zero. Please check your yaml if you expect a different behaviour")
continue
if not filecmp.cmp(fileini, fileout) or perturbed_from_reference:
# Fetch values from fileini
spec2read = spec if not perturbed_from_reference else ref_trid[1]
with _hdf5_lock:
with Dataset(fileini, "r") as ds:
ljust_specs_in = ds.variables["species"][:].astype(str)
specs_in = ["".join(p).strip() for p in ljust_specs_in]
# Skip if spec not in fileini
if spec2read not in specs_in:
warning(
f"{spec2read} is not accounted for in the boundary conditions file {fileini}. Please check your LBC whether it should be")
continue
ispec = specs_in.index(spec2read)
if input_type == "latcond":
lbcin = ds.variables["lat_conc"][..., ispec]
lbcin = lbcin[..., np.newaxis, :]
else:
lbcin = ds.variables["top_conc"][..., ispec]
lbcin = lbcin[:, np.newaxis, ...]
lbcin = xr.DataArray(
lbcin,
coords={"time": list_dates},
dims=("time", "lev", "lat", "lon"),
)
# Check that value is not the same in fileout
with Dataset(fileout, "r") as ds_out:
ljust_specs_out = ds_out.variables["species"][:].astype(
str)
specs_out = ["".join(p).strip()
for p in ljust_specs_out]
# Skip if spec not in fileini
if spec2read not in specs_out:
debug(
f"{spec2read} is not available in {fileout}. "
f"Just fetching from {fileini}"
)
else:
ispec = specs_out.index(spec2read)
if input_type == "latcond":
lbcout = ds_out.variables["lat_conc"][..., ispec]
lbcout = lbcout[..., np.newaxis, :]
else:
lbcout = ds_out.variables["top_conc"][..., ispec]
lbcout = lbcout[:, np.newaxis, ...]
lbcout = xr.DataArray(
lbcout,
coords={"time": list_dates},
dims=("time", "lev", "lat", "lon"),
)
if np.all(lbcout == lbcin):
continue
# If lbc file is still a link, should be copied
# to be able to modify it locally
if os.path.islink(fileout):
file_orig = pathlib.Path(fileout).resolve()
os.unlink(fileout)
shutil.copy(file_orig, fileout)
# Now saves to buffer variable before saving
data2dump[spec] = copy.deepcopy(lbcin)
# Repeat operations for tl
if mode == "tl":
# If does not exist, just copy and put everything to 0
if not os.path.isfile(fileoutincr):
shutil.copy(fileout, fileoutincr)
with _hdf5_lock:
with Dataset(fileoutincr, "a") as ds:
ds.variables["lat_conc"][:] = \
0. * ds.variables["lat_conc"][:]
ds.variables["top_conc"][:] = \
0. * ds.variables["top_conc"][:]
# If spec not in reference file, just pass
with _hdf5_lock:
with Dataset(fileini, "r") as ds:
ljust_specs_in = ds.variables["species"][:].astype(str)
specs_in = ["".join(p).strip().lower()
for p in ljust_specs_in]
if spec.lower() not in specs_in:
continue
# Load list of species in BOUN_CONCS.nc
with _hdf5_lock:
with Dataset(fileoutincr, "r") as ds:
ljust_specs_in = ds.variables["species"][:].astype(str)
specs_in = ["".join(p).strip().lower()
for p in ljust_specs_in]
# If species already in increment file, it is already at 0
if spec.lower() in specs_in:
continue
# Fetch correct structure from BOUN_CONCS.nc and put to 0
with _hdf5_lock:
with xr.open_dataset(fileoutincr) as ds:
if spec.lower() in specs_in:
ispec = specs_in.index(spec.lower())
if input_type == "latcond":
lbcin = ds["lat_conc"][..., ispec].values
lbcin = lbcin[..., np.newaxis, :]
else:
lbcin = ds["top_conc"][..., ispec].values
lbcin = lbcin[:, np.newaxis, ...]
# Stop here if already all zeros
if np.all(lbcin == 0):
continue
else:
if input_type == "latcond":
lbcin = np.zeros(
(len(list_dates), self.domain.nlev,
1, self.domain.nlon_side)
)
else:
lbcin = np.zeros(
(len(list_dates), 1, self.domain.nlat, self.domain.nlon))
lbcin = xr.DataArray(
0. * lbcin,
coords={"time": list_dates},
dims=("time", "lev", "lat", "lon"),
)
# Dump to BOUN_CONCS.increments.nc
debug(f"Replacing increments values to 0 for {spec} in {fileoutincr}")
data2dump_tl[spec] = copy.deepcopy(lbcin)
else:
# Replace existing link by copy
# of original file to modify it
path.copyfromlink(fileout)
# Write boundary conditions
lbc_fwd = tracer_data["spec"]
data2dump[spec] = copy.deepcopy(lbc_fwd)
if mode == "tl":
path.copyfromlink(fileoutincr)
lbc_tl = tracer_data.get("incr", 0.0 * lbc_fwd)
self.latcond.write(
spec, fileoutincr, lbc_tl, comp_type=input_type
)
data2dump_tl[spec] = copy.deepcopy(lbc_tl)
# Now dump buffer data if any
if len(data2dump) > 0:
self.latcond.write(
list(data2dump.keys()), fileout, data2dump, comp_type=input_type
)
# Now dump buffer data if any
if len(data2dump_tl) > 0:
self.latcond.write(
list(data2dump_tl.keys()), fileoutincr, data2dump_tl, comp_type=input_type
)
# Check that the dates are consistent with what CHIMERE expects
replace_dates(fileout, list_dates, self.ignore_input_dates)
if mode == "tl":
replace_dates(fileoutincr, list_dates, self.ignore_input_dates)
