Commit da2229a3 authored by Antoine Berchet's avatar Antoine Berchet
Browse files

Fix FLEXPART inicond when missing particles; background from CAMS in FLEXPART

parent e080d741
......@@ -2,6 +2,7 @@ import os
import numpy as np
import pandas as pd
import datetime
from netCDF4 import Dataset
from pycif.utils import path
......@@ -38,16 +39,20 @@ def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None,
list_files = {}
list_dates = {}
for dd in list_period_dates:
file = dd.strftime("{}/{}".format(ref_dir, ref_file))
# Fetch date frequency
with Dataset(file, "r") as f:
ntimes = f.dimensions["time"].size
# to_timedelta does not work with all frequencies!
datef = dd + \
datetime.timedelta(
days=int(pd.DatetimeIndex([dd]).days_in_month[0]))
list_hours = pd.date_range(
dd, datef, freq="6H", closed="left")
list_dates[dd] = [[hh, hh + datetime.timedelta(hours=6)]
for hh in list_hours]
file = dd.strftime("{}/{}".format(ref_dir, ref_file))
list_hours = pd.date_range(dd, datef, periods=ntimes + 1)
list_dates[dd] = [[hh0, hh1 + datetime.timedelta(hours=6)]
for hh0, hh1 in zip(list_hours[:-1], list_hours[1:])]
list_files[dd] = (len(list_hours) * [file])
target_file = "{}/{}".format(target_dir, os.path.basename(file))
......
......@@ -61,11 +61,9 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
latc = np.linspace(lat_min - dy/2., lat_max + dy/2., nlat + 1)
# Read vertical information in domain_file
sigma_a = nc["hyai"].values
sigma_b = nc["hybi"].values
print('XXXXX check if hybrid coefficients must be for the interfaces or the middles of the levelsXXXXX')
nlevs = sigma_a.size - 1
sigma_a = nc["hyam"].values
sigma_b = nc["hybm"].values
nlevs = sigma_a.size
# Initializes domain
setup = Setup.from_dict(
......@@ -98,6 +96,8 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
setup.domain.zlat = zlat
setup.domain.zlonc = zlonc
setup.domain.zlatc = zlatc
# Cyclic domain in longitude
setup.domain.lon_cyclic = True
return setup.domain
import datetime
import os
import pandas as pd
import numpy as np
import xarray as xr
from netCDF4 import Dataset
......@@ -31,6 +31,11 @@ def read(
"""
# Variable name to extract
var2extract = name
if varnames != "":
var2extract = varnames
# Reading fields for periods within the simulation window
xout = []
opened_file = ""
......@@ -40,11 +45,13 @@ def read(
ds = xr.open_dataset(dd_file)
opened_file = dd_file
date_index = int((dd[0] - ddi) / datetime.timedelta(hours=6))
ntimes = ds.dims["time"]
freq = dd[0].days_in_month * 24 / ntimes
date_index = int((dd[0] - ddi) / datetime.timedelta(hours=freq))
# bottom of the atmosphere = at the beginning of the table
lat = ds['latitude']
conc = ds[varnames].values[date_index]
conc = ds[var2extract].values[date_index]
if lat[1] < lat[0]:
conc = conc[:, :, ::-1, :]
......@@ -55,4 +62,5 @@ def read(
coords={"time": np.array(dates)[:, 0]},
dims=("time", "lev", "lat", "lon"),
)
return xmod
......@@ -40,6 +40,9 @@ def inicond_contribution(self, mode, dataobs,
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()
......@@ -61,7 +64,7 @@ def inicond_contribution(self, mode, dataobs,
ini_sim = (
dataini[inds_inicond, :, istartsensit:, 0].values * ini_sensit
).sum()
# Filling simulation
sim_col = dataobs.columns.get_indexer([data_id])
dataobs.iloc[obs_i, sim_col] += ini_sim
......@@ -74,3 +77,6 @@ def inicond_contribution(self, mode, dataobs,
sim_col = dataobs.columns.get_indexer(["background_tl"])
dataobs.iloc[obs_i, sim_col] = ini_sim
......@@ -214,6 +214,9 @@ def native2inputs_adj(
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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment