import numpy as np
import pandas as pd
import datetime
from ......utils.datastores.empty import init_empty
from ..utils.point_data import read_point_output, read_point_input, calc_and_write_point_dep
# JvP 20210705: added import statements, MODULE_NAME module level logger
import logging
import os
import sys
import netCDF4 as nc4
from ......utils.check.errclass import CifError, CifRuntimeError
MODULE_NAME = __name__[__name__.index(
'TM5'):] if 'TM5' in __name__ else __name__
logger = logging.getLogger(MODULE_NAME)
# Original make_obs function by A. Berchet
[docs]
def make_obs_AB(self, datastore, runsubdir, mode, tracer, do_simu=True):
"""Write the TM5 observation file for one sub-period and one tracer.
Formats observation metadata (station ID, level, coordinates, parameter,
time step) and writes to the TM5 point-data observation file under
*runsubdir*. For adjoint mode, also writes the departure vector.
Args:
self: TM5 model plugin instance.
datastore: CIF data-store for the observation tracer.
runsubdir (str): path to the period run directory.
mode (str): ``'fwd'``, ``'tl'``, or ``'adj'``.
tracer (str): tracer / species name.
do_simu (bool): if ``False``, skip writing and only update bookkeeping.
"""
# If empty datastore, do nothing
if datastore.size == 0:
return
# # JvP 20210602, telecon with AB: only relevant for multiple species
# and/or multiple datastreams (e.g. point and sat obs)
# # Otherwise, crop the datastore to active species
# # Save it to the model in case it is needed later
# if not hasattr(self, "dataobs") or getattr(self, "reset_obs", True):
# self.dataobs = {spec: init_empty()
# for spec in self.chemistry.acspecies.attributes}
#
# self.dataobs[tracer] = pd.concat([self.dataobs[tracer], datastore],
# axis=0, sort=False)
# If do not need to do CHIMERE simulation, just update obs datastore
if not do_simu:
return
# Write only species simulated by the model
mask = (
# AB to JvP (dd.20210602): replace self.dataobs[tracer] with
# datastore from here on...
# self.dataobs[tracer]["parameter"]
datastore["parameter"]
.str.upper()
.isin(["CH4"])
)
# AB to JvP (dd.20210602): replace self.dataobs[tracer] with datastore...
# data2write = self.dataobs[tracer].loc[mask]
data2write = datastore.loc[mask]
# For adjoint, check that there is no NaN values
if mode == "adj":
if not np.all(~np.isnan(data2write["maindata"].loc[:, "adj_out"])):
raise CifError("WARNING: pycif will drive TM5 adjoint "
"with NaNs values! Check prior informations")
# Then write the data
# pyCIF needs observations overlapping several model time steps to be
# unfolded
# That mean that if "dtstep" > 1, you need to create as many virtual
# observations
# in the observation file as dtstep
# Refer to CHIMERE or LMDZ for examples
# New make_obs function by J.C.A. van Peet
[docs]
def make_obs(self, datastore, runsubdir, mode, tracer, ddi, do_simu=True):
"""
PURPOSE
Write point observations to file readable by TM5
ARGS
No idea...
KWARGS
No idea...
NOTE
This function is called from the function outputs2native_adj, which is
defined in pycif/plugins/models/TM5/io/outputs2native_adj.py
VERSION HISTORY
2.3 28-10-2021 by J.C.A. van Peet
*) Updated determination of iter_nr.
*) When mode=adj, the CIF observation increments ('adj_out') are now
written to the point departure file instead of the TM5 forcings
(which are retained for comparison purposes).
2.2 22-10-2021 by J.C.A. van Peet
Antoine Berchet added ddi to the function header. No idea what it does
though...
2.1 20-07-2021 by J.C.A. van Peet
*) Updated the way time is selected from the data2write variable
*) Added on-the-fly generation of stationlist_CH4.dat, which is
then used in params.py for the variable station.list.file.
2.0 05-07-2021 by J.C.A. van Peet
Adapted for TM5.
1.0 ??-??-???? by A. Berchet
See original code above.
"""
# Set the name of this function
PROG_NAME = MODULE_NAME+".make_obs"
# Local logger
logger = logging.getLogger(PROG_NAME)
logger.setLevel(logging.DEBUG)
# Some debug statements...
logger.debug("")
logger.debug("*"*30)
logger.debug(PROG_NAME+" => DEBUG:")
logger.debug(f" datastore = {datastore}")
logger.debug(f" runsubdir = {runsubdir}")
logger.debug(f" mode = {mode}")
logger.debug(f" tracer = {tracer}")
logger.debug(f" kwarg do_simu = {do_simu}")
# If empty datastore, do nothing
if datastore.size == 0:
return
# JvP 20210705: based on a telecon I had with AB (dd. 2-6-2021), I commented
# out some parts of the original code. I did not copy those parts here, but
# they are included for future reference in the original code above.
# No idea what do_simu does. In the original code it's just checked if it
# is True. If not, you return to the calling function. However, I would like
# to know if that ever happens, so I raise an error to stop CIF if that ever
# happens.
if not do_simu:
# Original
# return
# JvP 20211022: Originial code. Judging by the code by Antoine below,
# you should just return if do_simu == False. But I still would like
# to know if that ever happens, so I reverted back to my old code.
logger.critical("*"*30)
logger.critical(PROG_NAME+" => ERROR: do_simu is False.")
logger.critical(" How to proceed???")
logger.critical(" Computer says no...")
logger.critical("*"*30)
raise CifRuntimeError
#
# Code Antoine used to replace my code above.
# logger.critical("")
# logger.critical("*"*30)
# logger.critical(PROG_NAME + " => ERROR: do_simu is False.")
# logger.critical(" Do nothing")
# logger.critical("*"*30)
# logger.critical("")
# return
# end if do_simu
# Get a mask with only species simulated by the model (i.e. CH4). Then use
# the dataframe.loc method/object/property (whatever, see
# https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.loc.html
# for more info) to select only the rows from the dataframe where the
# "parameter" is CH4. I think that the resulting data2write is another
# dataframe...
mask = datastore["metadata"]["parameter"].str.upper().isin(["CH4"])
data2write = datastore.loc[mask]
# >>> JvP 20211022
# The following lines were added by Antoine. No idea what they do...
# Unfold data2write depending on the duration of each observation
inds = [0] + list(np.cumsum(
data2write[("metadata", "dtstep")].values[:-1]))
ref_dtstep = data2write["metadata"]["dtstep"].values
data2write = data2write.reset_index().loc[
data2write.index.repeat(data2write["metadata"]["dtstep"])]
idx = np.zeros(len(data2write), dtype=int)
idx[inds] = inds
np.maximum.accumulate(idx, out=idx)
data2write.loc[:, ("metadata", "tstep")
] += np.arange(len(data2write)) - idx
data2write[("metadata", "date")] = \
self.tstep_dates[ddi][data2write["metadata"]["tstep"].astype(int)] \
+ datetime.timedelta(minutes=30)
# <<< JvP 20211022.
# Number of observations in data2write
n_obs = len(data2write)
# Apparently, in the adjoint run, there can be NaN values in the data.
# So for the adjoint, check that there are no NaN values and raise a
# RuntimeError if there are.
if mode == "adj":
if not np.all(pd.notnull(data2write["maindata"].loc[:, "adj_out"])):
raise CifError("WARNING: pycif will drive TM5 adjoint "
"with NaNs values! Check prior informations")
# end if
# end if
# Generate station data from data2write:
# *) station names can be any length
# *) station codes are strings of 11 characters (e.g. ALT_NOA_000), which
# for now I'll just set to 'STN_XXX_000'
# *) station_type seems to be 'FM' for all stations (Flask, M...?)
# *) station_id is a 1 based index linking the measurements to the
# unique stations
# *) station_tower seems to be 0 for all stations (no tower measurements
# assimilated...)
# *) Use dtype='S' to be able to write to netCDF as character arrays,
# i.e. something that Fortran can understand...
# If you don't use 'S' as dtype, in Python3 the strings will be
# unicode by default, and that messes up the stringtochar function
# somehow.
#
# **************************************************************************
# * According to a mail by Antoine (dd. 15-07-2021), the format of the *
# * "station" column in data2write changes according to the version of *
# * CIF that you're using: *
# * Do you generate observations from random values? or from a monitor? *
# * With the random values, stations just have a number (in your *
# * version it is indeed an integer, which is a problem; in newest *
# * version, it is a number but as a string) *
# * When you lad from a custom monitor, you can give any name to your *
# * station. *
# * In other words, the following may crash after some CIF update... *
# **************************************************************************
#
# This is fixed below:
metadata = data2write["metadata"]
station_full = pd.Series(data=np.asarray(metadata['station']), dtype=str)
station_list = pd.Series(index=metadata['station'].drop_duplicates().values,
data=range(len(metadata['station'].drop_duplicates())))
station_name = np.array("station_" + station_full, dtype="S")
station_code = np.array(
("STN_" + station_full + "_000").str.zfill(11), dtype="S")
station_type = np.array(len(metadata) * ['FM'], dtype="S")
station_id = np.array(station_list[metadata['station']] + 1)
station_tower = np.array(len(metadata) * [0])
logger.debug(f" station_name = {station_name}")
logger.debug(f" station_code = {station_code}")
logger.debug(f" station_type = {station_type}")
logger.debug(f" station_id = {station_id}")
logger.debug(f" station_tower = {station_tower}")
# Unique stations in station_name
unique_stn, unique_stn_idx = np.unique(station_name, return_index=True)
unique_stn_alt = metadata['alt'].iloc[unique_stn_idx].to_numpy()
unique_stn_code = station_code[unique_stn_idx]
unique_stn_lat = metadata['lat'].iloc[unique_stn_idx].to_numpy()
unique_stn_lon = metadata['lon'].iloc[unique_stn_idx].to_numpy()
unique_stn_name = station_name[unique_stn_idx]
unique_stn_tower = station_tower[unique_stn_idx]
unique_stn_type = station_type[unique_stn_idx]
logger.debug(f" unique_stn = {unique_stn}")
logger.debug(f" unique_stn_idx = {unique_stn_idx}")
logger.debug(f" unique_stn_alt = {unique_stn_alt}")
logger.debug(f" unique_stn_code = {unique_stn_code}")
logger.debug(f" unique_stn_lat = {unique_stn_lat}")
logger.debug(f" unique_stn_lon = {unique_stn_lon}")
logger.debug(f" unique_stn_name = {unique_stn_name}")
logger.debug(f" unique_stn_tower = {unique_stn_tower}")
logger.debug(f" unique_stn_type = {unique_stn_type}")
# Original comment by Antoine... I guess I'll check for dtstep > 1 and
# just raise an exception if that occurs...
#
# Then write the data
# pyCIF needs observations overlapping several model time steps to be unfolded
# That mean that if "dtstep" > 1, you need to create as many virtual observations
# in the observation file as dtstep
# Refer to CHIMERE or LMDZ for examples
# Set the name of the file
point_file = runsubdir + '/point_input.nc4'
logger.debug(" point_file = " + point_file)
# Delete any pre-existing files, and open the output file. Usually,
# selecting "w" when opening a file should create a new file overwriting
# any existing file. But if the file is opened in hdfview (which happens
# regularly), this issues an error. Deleting the file explicitly supresses
# that error. But os.remove raises an error if the file does not exist
# (at least in some Python versions I've tried this), so wrap it in a try
# statement to suppress it.
try:
os.remove(point_file)
except:
pass
# end try
# Open file
nc_id = nc4.Dataset(point_file, 'w')
# Global attributes
nc_id.listfile = 'Generated by pyCIF'
nc_id.obsdatadir = 'Generated by pyCIF'
# >>> Create dimensions
# idate = 6: year, month, day, hour, minute, second
nc_id.createDimension('idate', 6)
# id = n_obs: an integer array with indices of the observations...
nc_id.createDimension('id', n_obs)
id_var = nc_id.createVariable('id', np.int32, ('id',))
id_var[()] = np.arange(n_obs, dtype=np.int32) + 1
id_var.longname = 'observation index (1-based)'
# Note that max(...,key=len) gives you the longest string in the array,
# and len returns the length of that string
nc_id.createDimension('len_station_code', len(max(station_code, key=len)))
nc_id.createDimension('len_station_name', len(max(station_name, key=len)))
nc_id.createDimension('len_station_type', len(max(station_type, key=len)))
nc_id.createDimension('len_tracer', len(tracer))
nc_id.createDimension('station', unique_stn.size)
# <<< Create dimensions
# Create variables
alt_var = nc_id.createVariable('alt', np.float64, ('id',))
bias_id_var = nc_id.createVariable(
'bias_id', np.int32, ('id',))
date_components_var = nc_id.createVariable(
'date_components', np.int16, ('id', 'idate'))
lat_var = nc_id.createVariable('lat', np.float64, ('id',))
lon_var = nc_id.createVariable('lon', np.float64, ('id',))
mixing_ratio_var = nc_id.createVariable(
'mixing_ratio', np.float64, ('id',))
mixing_ratio_error_var = nc_id.createVariable(
'mixing_ratio_error', np.float64, ('id',))
mixing_ratio_stdv_var = nc_id.createVariable(
'mixing_ratio_stdv', np.float64, ('id',))
sampling_strategy_var = nc_id.createVariable(
'sampling_strategy', np.int16, ('id',))
station_alt_var = nc_id.createVariable(
'station_alt', np.float32, ('station',))
station_code_var = nc_id.createVariable(
'station_code', 'S1', ('station', 'len_station_code'))
station_id_var = nc_id.createVariable(
'station_id', np.int32, ('id',))
station_lat_var = nc_id.createVariable(
'station_lat', np.float32, ('station',))
station_lon_var = nc_id.createVariable(
'station_lon', np.float32, ('station',))
station_name_var = nc_id.createVariable(
'station_name', 'S1', ('station', 'len_station_name'))
station_tower_var = nc_id.createVariable(
'station_tower', np.int32, ('station',))
station_type_var = nc_id.createVariable(
'station_type', 'S1', ('station', 'len_station_type'))
tower_var = nc_id.createVariable('tower', np.int32, ('id',))
tracer_var = nc_id.createVariable(
'tracer', 'S1', ('len_tracer',))
# Fill variables
alt_var[()] = metadata['alt']
bias_id_var[()] = np.ones((n_obs,), dtype=np.int32)
#
# JvP 20210720
# Apparently, the date column contains np.datetime64[ns] objects, so that
# does include a time value. I don't know why I set the time to 12 to start
# with, but you can just use i.[hour|minute|second] in the command below.
# More info on np.datetime64: https://numpy.org/doc/stable/reference/arrays.datetime.html
# Note that this interface is stable only since the latest numpy version
# (1.21), and that older versions considered this interface as experimental
# https://numpy.org/doc/1.20/reference/arrays.datetime.html.
# date_components_var [()] = np.array( [[i.year, i.month, i.day, 12, 0, 0] for i in data2write['date']] )
date_components_var[()] = np.array(
[[i.year, i.month, i.day, i.hour, i.minute, i.second]
for i in metadata['date']])
# Here, replace lat/lon by the middle of grid cells, as the rest of the CIF does
# the interpolation
lat_var[()] = self.domain.zlat[
metadata["i"].astype(int), metadata["j"].astype(int)]
lon_var[()] = self.domain.zlon[
metadata["i"].astype(int), metadata["j"].astype(int)]
mixing_ratio_var[()] = data2write["maindata"]['obs'] # should be in ppbv
#
# NOTE: MIXING_RATIO_ERROR SEEMS TO BE A CONSTANT 3.0 IN THE POINT_INPUT FILE...
mixing_ratio_error_var[()] = np.full((n_obs,), 3.0)
#
# NOTE: THERE'S NO ERROR IN THE OBSERVATIONS, SO I USE CONSTANT OF 0.0001 FOR NOW...
mixing_ratio_stdv_var[()] = 0.001 * data2write["maindata"]['obs']
#
# NOTE: SAMPLING_STRATEGY SEEMS TO BE A CONSTANT 4 IN THE POINT_INPUT FILE...
sampling_strategy_var[()] = np.full((n_obs,), 4)
#
station_alt_var[()] = unique_stn_alt
station_alt_var.long_name = 'station altitude'
station_code_var[()] = nc4.stringtochar(unique_stn_code)
station_code_var.long_name = 'station code'
station_id_var[()] = station_id
station_id_var.long_name = 'station index (1-based)'
station_lat_var[()] = unique_stn_lat
station_lat_var.long_name = 'station latitude'
station_lon_var[()] = unique_stn_lon
station_lon_var.long_name = 'station longitude'
station_name_var[()] = nc4.stringtochar(unique_stn_name)
station_name_var.long_name = 'station name'
station_tower_var[()] = unique_stn_tower
station_tower_var.long_name = 'tower flag: 1=tower, 0=surface'
station_type_var[()] = nc4.stringtochar(unique_stn_type)
station_type_var.long_name = 'station observation type'
tower_var[()] = station_tower
tracer_var[()] = nc4.stringtochar(np.array(tracer, dtype='S'))
# Close file
nc_id.close()
# *********************************************************
# * QUESTIONS TO ANTOINE: *
# * -) Is there a way to get a name for the stations? *
# * Currently, there are only numbers in the station *
# * column which I convert to the name "station_00X". *
# * -) What is the unit of the 'obs' column in datastore? *
# * -) Is there an error for the 'obs' column? *
# * -) What is the time of the observations? Currently, *
# * only the date is given. I've set the time to 12 *
# * for now. *
# * -) How to access input_dir or other variables from *
# * the yaml file? I now copied the global input_dir *
# to the model plugin input_dir variable *
# *********************************************************
# >>> JvP 20210720
# Write station locations to a stationlist_CH4[...].dat file that is
# also used in params.py for the station.list.file variable.
# Set the name of the stationlist file
stationlist_file = runsubdir + '/stationlist_CH4.dat'
# Write data to the stationlist file
with open(stationlist_file, 'w') as f_id:
# Write header line
print(
'ID_OBS LAT LON ALT TP TO LT_min LT_max STATIONNAME', file=f_id)
# For each of the stations, write the corresponding variables
# Notes:
# *) the columns 'TP' through 'LT_max' have fixed values, corresponding
# to those in the original stationlist file.
# *) the arrays unique_stn_code and unique_stn_name are converted to
# unicode with 'astype(str)', because else their contents in the
# file look like "b'STN_000_000'"
for stn_id, stn_lat, stn_lon, stn_alt, stn_name in \
zip(unique_stn_code.astype(str), unique_stn_lat,
unique_stn_lon, unique_stn_alt, unique_stn_name.astype(str)):
print("%s %7.2f %7.2f %7.1f %s %2i %6.1f %6.1f %s" %
(stn_id, stn_lat, stn_lon, stn_alt, 'FM', 0, 0.0, 24.0, stn_name), file=f_id)
# end for
# end with
# <<< JvP 20210720
# Some debug statements...
logger.debug("*"*30)
logger.debug("")
# Now should write the departure file when in adjoint mode
if mode == "adj":
# >>> JvP 20211022
# After a successful forward run, you have to make a "point_departures.nc4"
# file. In TM5, this file is made in the "dep" job by TMVar.py, which calls
# Observations_JRC.PointObs_JRC.applyPointObs(...) as:
# point.applyPointObs( trackfile=trackfile, trackfile_apri=trackfile_apri,
# mismatchfile=departurefile,
# Params=None, montecarlo=False )
# where:
# trackfile = some/path/to/point_output.nc4 (written by TM5)
# trackile_apri = some/path/to/point_output.nc4 from the apri run (i.e. iter 1)
# departurefile = some/path/to/point_departures.nc4 (which is read by the
# adjoint run of TM5, and which is set in params.py when
# the config file is written).
# Message
logger.info(PROG_NAME+" => writing point departure file.")
# >>> JvP 20211028
# Get the iteration number from the runsubdir, since it is not stored
# somewhere easily accessible. According to Antoine (TM5-CIF telecon
# 27-10-2021) the iteration number is -1 for the apriori run (or first
# iteration), -2 for the aposterior run (or final iteration) and 2,3,...
# for the 'normal' iterations. This means that when you convert the
# iteration number to a string with '%04i'%iter_nr, the apriori iter_nr
# number will become '-001', the aposteriori iter_nr will become '-002',
# and the rest will become '0002', '0003' etc. If you combine the iter_nr
# with a prefix like 'fwd_', the strings become (in chronological order):
# 'fwd_-001', 'fwd_0002', 'fwd_0003', ..., 'fwd_-002'.
# Notes:
# *) Luckily, you only need the string representation of this number, so
# it is sufficient to extract the part after '[fwd|adj]_'.
# *) Similar code is used in params.write_tm5_rc(...)
# Split runsubdir into its consituent parts
runsubdir_parts = runsubdir.split(sep='/')
# Check if runsubdir_parts[-2] starts with either 'fwd_' or 'adj_'.
if (not runsubdir_parts[-2].startswith(('fwd_', 'adj_'))):
logger.critical("")
logger.critical("*"*30)
logger.critical(
PROG_NAME+" => ERROR: path component does not start with 'fwd_' or 'adj_'!")
logger.critical(" runsubdir = "+runsubdir)
logger.critical(f" runsubdir parts = {runsubdir_parts}")
logger.critical(" Computer says no...")
logger.critical("*"*30)
logger.critical("")
raise CifRuntimeError
# end if
# Get the string representing the iteration number, and log a message.
iter_nr = runsubdir_parts[-2][4:]
logger.info(PROG_NAME+f" => iter_nr = {iter_nr}")
# <<< JvP 20211028
# Get the a-priori, forward, and adjoint directories
# Here, use self.adj_refdir defined in run.py, to be sure that it is
# the last computed forward
apri_dir = '/'.join(runsubdir_parts[0:-2]) + \
'/fwd_-001/' + runsubdir_parts[-1]+'/'
# fwd_dir = '/'.join(runsubdir_parts[0:-2])+'/fwd_'+iter_nr+'/'+runsubdir_parts[-1]+'/'
fwd_dir = f'{self.adj_refdir}/{runsubdir_parts[-1]}/'
adj_dir = '/'.join(runsubdir_parts[0:-2]) + \
'/adj_'+iter_nr+'/'+runsubdir_parts[-1]+'/'
#
logger.info(PROG_NAME+f" => apri_dir = {apri_dir}")
logger.info(PROG_NAME+f" => fwd_dir = {fwd_dir}")
logger.info(PROG_NAME+f" => adj_dir = {adj_dir}")
# Read the point_output.nc4 file for the current iteration.
# You can check if your functions work the same as the originals by
# using point files from a TM5 run.
# point_output_file = '/home/jpt930/NOBACKUP/tm5/20210222-tm5_sf-18m_points_only-ERA5-spiv/tmvar/var4d/iter-0001/fwd/output/point_output.nc4'
point_output_file = fwd_dir+'output/point_output.nc4'
logger.info(PROG_NAME+f" => point_output_file = {point_output_file}")
point_output_data = read_point_output(point_output_file)
# Read the point_output.nc4 file from the a-priori (i.e. iteration 1)
# You can check if your functions work the same as the originals by
# using point files from a TM5 run.
# apri_point_output_file = '/home/jpt930/NOBACKUP/tm5/20210222-tm5_sf-18m_points_only-ERA5-spiv/tmvar/var4d/iter-0001/fwd/output/point_output.nc4'
# AB: Is this part needed? we can simply take point_output_data?
# apri_point_output_file = apri_dir+'output/point_output.nc4'
# logger.info(PROG_NAME+" => apri_point_output_file = %s" % apri_point_output_file)
# apri_point_output_data = read_point_output(apri_point_output_file)
#
# # Check if the sample_ids match
# if( not np.array_equal( point_output_data['sample_id'], apri_point_output_data['sample_id'] ) ):
# logger.critical("")
# logger.critical("*"*30)
# logger.critical(PROG_NAME+" => ERROR: Something wrong with sample_id!")
# logger.critical(" point_output_data['sample_id'] = %s" % point_output_data['sample_id'] )
# logger.critical(" apri_point_output_data['sample_id'] = %s" % apri_point_output_data['sample_id'] )
# logger.critical(" mode = %s" % mode )
# logger.critical(" runsubdir = %s" % runsubdir )
# logger.critical(" iter_nr = %s" % iter_nr )
# logger.critical(" Computer says no...")
# logger.critical("*"*30)
# logger.critical("")
# raise RuntimeError(PROG_NAME+" => ERROR: something wrong with sample_id!")
# # end if
# Read the point_input.nc4 file from the last forward
# You can check if your functions work the same as the originals by
# using point files from a TM5 run.
# apri_point_input_file = '/home/jpt930/NOBACKUP/tm5/20210222-tm5_sf-18m_points_only-ERA5-spiv/input/point_input.nc4'
fwd_point_input_file = fwd_dir+'point_input.nc4'
logger.info(PROG_NAME+f" => fwd_point_input_file = {fwd_point_input_file}")
fwd_point_input_data = read_point_input(
fwd_point_input_file, point_output_data['sample_id'])
# Set point departure filename
point_dep_file = runsubdir+'/point_departures.nc4'
logger.info(PROG_NAME+f" => point_dep_file = {point_dep_file}")
# Calculate the departure data
#
# JvP 20211028: Antoine entered the following example code when he made
# some modifications to the TM5 branch:
# print(__file__)
# import code
# code.interact(local=dict(locals(), **globals()))
#
# # At this point the variable data2write["obs_incr"], includes:
# # R^(-1) x (H(xb) - y)
# The first three lines just exit to an interactive prompt, including
# all variables and their values at that point in the code. Nice
# debugging tip. The value of obs_incr is probably the same as the
# TM5 'forcing' calculated in the function calc_and_write_point_dep in
# .../pycif/plugins/models/TM5/io/utils/point_data.py. So I added
# obs_incr as keyword to the function call. Inside the function,
# obs_incr is written to file as 'forcing' while the original TM5
# focings are retained in the 'forcing_TM5' variable for comparison
# purposes.
#
# calc_and_write_point_dep( point_output_data, apri_point_output_data,
# apri_point_input_data, point_dep_file )
#
calc_and_write_point_dep(point_output_data, point_output_data,
fwd_point_input_data, point_dep_file,
adj_out=data2write["maindata"]["adj_out"])
# <<< JvP 20211022
# end if (mode == 'adj')
# logger.debug("")
# logger.debug("*"*30)
# logger.debug(PROG_NAME+" => Computer says no!")
# logger.debug("*"*30)
# logger.debug("")
# try:
# raise RuntimeError
# except RuntimeError as e:
# #logger.exception("OOPS!")
# logger.critical(e, exc_info=True)
# #raise # => Will display the traceback on screen a second time
# sys.exit() # => Just exit.
# end try
# end function make_obs
