import datetime
import glob
import os
from logging import debug
import numpy as np
import pandas as pd
from ......utils.parallel import thread
from ...utils.read import read_flexpart_gridinit
from ......utils.check.errclass import CifError
[docs]
def inicond_contribution_ad(
self, mode, dataobs, fp_header_init, spec, ddi
):
"""Compute the adjoint initial-condition sensitivity from FLEXPART fields.
Reads initial-condition sensitivity grids and weights them by ``adj_out``
departures from *dataobs*, accumulating the result into the in-memory
initial-condition sensitivity store indexed by period and grid cell.
Args:
self: Lagrangian model plugin instance.
mode (str): ``'adj'``.
dataobs: CIF observation data-store (contains ``'adj_out'``).
fp_header_init: FLEXPART initial-condition sensitivity header.
spec (str): species name.
ddi (datetime): sub-simulation period start.
"""
inicond = {"spec": self.datainicond[("inicond", spec)][ddi]["spec"]}
nlat, nlon = self.domain.zlat.shape
nz = (fp_header_init.outheight != 0.).sum()
ini_dates = pd.DatetimeIndex(
self.datainicond[("inicond", spec)][ddi]["spec"].time.values
).to_pydatetime()
ndates =len(ini_dates)
# Execute parallel threads
nthreads = self.nthreads
thread_intervals = np.linspace(0, len(dataobs), nthreads + 1).astype(int)
@thread
def thread_function(ithread):
inicond_sensit_tmp = np.zeros((ndates, nz, nlat, nlon))
for obs_i in range(thread_intervals[ithread], thread_intervals[ithread + 1]):
process_obs_row(
self, dataobs, ithread,
fp_header_init,
ini_dates, inicond, inicond_sensit_tmp,
obs_i
)
return inicond_sensit_tmp
inicond_sensit = sum(thread_function(range(nthreads)))
# Flush inicond
self.datainicond[("inicond", spec)][ddi]["spec"] = None
return inicond_sensit
[docs]
def process_obs_row(
self, dataobs, ithread,
fp_header_init,
ini_dates,
inicond, inicond_sensit_tmp, obs_i
):
"""Compute the adjoint initial-condition sensitivity for a single observation row.
Reads the initial-condition sensitivity field for *obs_i* and weights it
by ``adj_out`` from *dataobs*, accumulating the result into
*inicond_sensit_tmp*.
Args:
self: Lagrangian model plugin instance.
dataobs: CIF observation data-store for the period.
ithread (int): thread index.
fp_header_init: FLEXPART initial-condition sensitivity header.
ini_dates: list of initial-condition sensitivity dates.
inicond (dict): initial-condition arrays.
inicond_sensit_tmp: thread-local accumulator for the sensitivity.
obs_i (int): row index into the observation data-store.
"""
row = dataobs.iloc[obs_i]["metadata"]
row_main = dataobs.iloc[obs_i]["maindata"]
station = row.station
network = row.network
subdir = row.date.strftime(self.footprint_dir_format)
# Translate station name if needed
if hasattr(self, "station_name_dict"):
station = self.dict_station_name[station.upper()].upper()
# Infer folder structure
runsubdir_init = os.path.join(
self.run_dir_bg,
self.footprint_stat_subdir_format.format(
stat=station.upper(), network=network),
subdir
)
release_date = row.date - pd.to_timedelta(self.release_shift)
file_date = release_date.strftime('%Y%m%d%H%M%S')
file_name = self.file_ini_format.format(
date=file_date, stat=station.upper(), network=network)
list_valid = glob.glob(os.path.join(runsubdir_init, file_name))
if list_valid == []:
debug(f"WARNING: file not found: {os.path.join(runsubdir_init, file_name)}")
return
elif len(list_valid) > 1:
raise CifError(
f"Multiple files fit the specified format {self.file_ini_format}. "
f"This can be related to the use of a wildcard... "
f"Please check your yml"
)
file_name = os.path.basename(list_valid[0])
debug(f"Thread #{ithread}: Reading {file_name} for station {station}")
grid_init = read_flexpart_gridinit(
runsubdir_init, file_name, fp_header_init)
# Normalize grid_init to make sure that total sensitivity is 1
grid_init /= grid_init.sum()
# Multiply 3-D sensitivity to background concentrations
# WARNING: do not deal with temporal and vertical dimension yet
nz = (fp_header_init.outheight != 0.).sum()
ini_sensit = grid_init.T.reshape(nz, -1)
inds_inicond = np.argmin(
(np.array(
[row.date -
pd.to_timedelta(self.backward_trajdays)]
)[:, np.newaxis] - ini_dates[np.newaxis, :]
) >= datetime.timedelta(0), axis=1) - 1
istartsensit = (
0 if self.domain.zlon.size == self.domain.zlon_in.size
else self.domain.zlon_in.size
)
inicond_sensit_tmp[inds_inicond, :, 0, istartsensit:] \
+= ini_sensit * row_main.adj_out
