Commit eaefb06a authored by Isabelle Pison's avatar Isabelle Pison
Browse files

lmdz kinetic beginning

parent e664dc70
import numpy as np
from pycif.utils.classes.setup import Setup
def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
date_ref = list(input_dates.values())[0][0]
print('Here, read the horizontal grid e.g. longitudes and latitudes')
#lon, lat = centers of the grid cells?
#nlon = lon.size
#nlat = lat.size
# minimum/maximum = including extreme borders of the grid cells?
#lon_min = lon.min() - (lon[1] - lon[0]) / 2 ?
#lon_max = lon.max() + (lon[-1] - lon[-2]) / 2
#lat_min = lat.min() - (lat[1] - lat[0]) / 2
#lat_max = lat.max() + (lat[-1] - lat[-2]) / 2
print('Here, read the vertical information, such as the number of vertical levels')
#nlevs =
# Reconstruct alpha and beta
print('Put here how to get at the vertical grid .e.g. sigma coefficients')
#sigma_a = np.empty(nlevs)
#sigma_b = np.empty(nlevs)
# 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,
"nlev": nlevs,
"sigma_a": sigma_a,
"sigma_b": sigma_b,
"pressure_unit": "Pa"
}
}
)
Setup.load_setup(setup, level=1)
return setup.domain
from .get_domain import get_domain
from .fetch import fetch
from .read import read
from .write import write
_name = "CAMS"
_version = "netcdf"
#requirements = {
# "domain": {
# "name": "CHIMERE",
# "version": "std",
# "empty": False,
# "any": False,
# },
#}
import os
import numpy as np
from pycif.utils import path
from .utils import find_valid_file
def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
list_files = {}
list_dates = {}
#print('fffffffffffff',input_dates)
for datei in input_dates:
tmp_files = []
tmp_dates = []
for dd in input_dates[datei]:
#print('ppppppppppp',dd)
files_cams, dates_cams = find_valid_file(ref_dir, ref_file, dd)
print(files_cams,dates_cams)
tmp_files.extend(files_cams)
tmp_dates.extend(dates_cams)
# print('FFFFFFFFFFFF',tmp_dates,tmp_files)
# Fetching
local_files = []
#print('RRRRRRRRRR',ref_file)
for f, dd in zip(tmp_files, tmp_dates):
target_file = "{}/{}".format(target_dir, dd.strftime(ref_file))
path.link(ref_dir+'/'+f, target_file)
local_files.append(target_file)
unique_dates, unique_index = np.unique(tmp_dates, return_index=True)
list_files[datei] = np.array(tmp_files)[unique_index]
list_dates[datei] = unique_dates
#print('oooooooooooooooo',list_dates)
return list_files, list_dates
import numpy as np
import glob
import datetime
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 bc_file in list_file:
try:
date = datetime.datetime.strptime(
os.path.basename(bc_file), ref_file
)
domain_file = bc_file
break
except ValueError:
continue
if domain_file is None:
raise Exception(
"CAMS domain could not be initialized as no file was found"
)
print('Domain file for CAMS BCs:',domain_file)
nc = xr.open_dataset(domain_file, decode_times=False)
# Read lon/lat in domain_file
lon = nc["longitude"]
lat = nc["latitude"]
nlon = lon.size
nlat = lat.size
#print(nlon,nlat)
#print(lon[0])
lon_min = lon.min()
lon_max = lon.max()
lat_min = lat.min()
lat_max = lat.max()
print(lon_min)
# Read vertical information in domain_file
sigma_a = nc["ap"].values
sigma_b = nc["bp"].values
#print('SSSSSSSSSSSSSSSS',sigma_a)
#print('ssssssssssssssss',sigma_b)
nlevs = sigma_a.size - 1
#print('LLLLLLLLLLLLLLLLLLLLLL',nlevs)
# 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,
"nlev": nlevs,
"sigma_a": sigma_a[1:],
"sigma_b": sigma_b[1:],
"pressure_unit": "Pa"
}
}
)
Setup.load_setup(setup, level=1)
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,
comp_type=None,
**kwargs
):
"""Get BCs from raw files and load them into a pyCIF
variables
Args:
self: the BC Plugin
name: the name of the component
tracdir, tracfile: raw data directory and file format
dates: list of dates to extract
comp_type: type of boundary conditions to generate
"""
print('ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ')
if type(tracfile) == str:
tracfile = [dd.strftime(tracfile) for dd in dates]
else:
if len(tracfile) != len(dates) :
raise Exception(
"Try read CAMS files from a list of dates and a "
"list of files, but not of same length:\n{}\n{}".format(
tracfile, dates
)
)
list_files = tracfile[:]
print('LLLLLLLLLLLLLL',list_files)
# Reading fields for periods within the simulation window
trcr_conc = []
times = []
for dd, dd_file in zip(dates, list_files):
file_conc = dd.strftime(dd_file)
dir_conc = dd.strftime(tracdir)
if not os.path.isfile("{}/{}".format(dir_conc, file_conc)) and getattr(
self, "closest_year", False
):
info(
"Warning: could not find correct year for CAMS; "
"using closest available one"
)
list_dates = [
datetime.datetime.strptime(os.path.basename(f), tracfile)
for f in glob.glob("{}/z_cams_l_nilu_*nc".format(dir_conc))
]
delta_dates = np.abs(dd - np.array(list_dates))
file_conc = list_dates[np.argmin(delta_dates)].strftime(tracfile)
nc = xr.open_dataset(
"{}/{}".format(dir_conc, file_conc), decode_times=False
)
#print('DDDDDDDDDDD',dd,nc['time'])
# Convert number of hours since first day of the month to datetime
times.extend(
[datetime.datetime(dd.year, dd.month, 1) + datetime.timedelta(hours=t)
for t in nc['time'].values]
)
#print('TTTTTTTTTT',times)
trcr_conc.append(nc[varnames]['time' == dates].values)
# Putting in DataArray for the right period = combien d'heures??? Juste dates???
# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
print('TTTTTTTTT',np.array(trcr_conc).shape)
print('ddddddddddd',dates,times[0])
xmod = xr.DataArray(
np.array(trcr_conc),
coords={"time": dates},
dims=("time", "lev", "lat", "lon"),
)
return xmod
import datetime
import glob
import os
import numpy as np
def find_valid_file(ref_dir, file_format, dd):
# Get all files and dates matching the file and format
list_files_orig = os.listdir(ref_dir)
list_dates = []
list_files = []
for f in list_files_orig:
try:
list_dates.append(datetime.datetime.strptime(f, file_format))
list_files.append(f)
except:
continue
list_files = np.array(list_files)
list_dates = np.array(list_dates)
# Sorting along dates
isort = np.argsort(list_dates)
list_dates = list_dates[isort]
list_files = list_files[isort]
if list_files == []:
raise Exception("Did not find any valid CAMS files in {} "
"with format {}. Please check your yml file"
.format(ref_dir, file_format))
#print('yyyyyyyyyyyyyyy',dd,list_dates,list_files)
# si date = premiere heure du mois -> fichier = celui du mois d'avant
# sinon, fichier = celui du mois
if dd.hour == 0 and dd.day == 1:
file_ref1 = list_files[0]
dd1 = list_dates[0]
else:
file_ref1 = list_files[1]
dd1 = list_dates[1]
file_ref2 = list_files[1]
dd2 = list_dates[1]
return [file_ref1, file_ref2], [dd1, dd2]
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",
)
......@@ -56,6 +56,14 @@ requirements = {
"type": "fields",
"newplg": True,
},
"kinetic": {
"name": "LMDZ",
"version": "photochem",
"empty": True,
"any": False,
"type": "fields",
"newplg": True,
},
"prodloss3d": {
"name": "LMDZ",
"version": "prodloss3d",
......@@ -78,6 +86,7 @@ required_inputs = [
"def",
"chem_fields",
"prescrconcs",
"kinetic",