Commit 36ef7c3c authored by Antoine Berchet's avatar Antoine Berchet
Browse files

Merge branch 'devel' into 'LSCE'

Devel

See merge request satinv/cif!169
parents d39d0d01 c4753402
......@@ -289,9 +289,10 @@ def doc_from_arginputs(input_arguments, write_headers=True):
towrite.append(" **accepted values**: {}".format(accepted))
elif type(accepted) is dict:
towrite.append(" **accepted values**:")
towrite.append("")
for key in accepted:
towrite.append(" {}: {}".format(key, accepted[key]))
towrite.append(" - {}: {}".format(key, accepted[key]))
else:
towrite.append(" **accepted type**: {}".format(accepted))
......
......@@ -19,6 +19,10 @@ mode:
name: gausscost
version: std
reload_from_previous: true
montecarlo:
nsample: 10
perturb_x: true
perturb_y: false
plugin:
name: 4dvar
version: std
......
......@@ -6,7 +6,7 @@ workdir: !join [*outdir, /inversion_long_bands_ensrf_]
datei: 2010-01-01
datef: 2010-01-05 00:00:00
mode:
nsample: 10
nsample: 50
plugin:
name: EnSRF
version: std
......
......@@ -219,11 +219,15 @@ NCERR(__LINE__)
NCERR(__LINE__)
idr=mm5date2numeric(datebuf)
idex = idate(ihourrun)
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS FILE'
print *,'IHOURRUN=',ihourrun,' EXPECTED= ' &
,idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
if (ignore_chck_dates == 0) then
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS FILE'
print *,'IHOURRUN=',ihourrun,' EXPECTED= ' &
,idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
endif
else
print *,'*** NO CHECK ON DATES in BOUN_CONCS'
endif
allocate(buf3(nspecboun, nhbound_domain, nverti))
stvec4 =(/1 , 1 , 1 , ihourrun+1 /)
......
......@@ -144,11 +144,15 @@ NCERR(__LINE__)
bounconc_ncid,bounconc_times_varid,datebuf,(/1,1/),(/dlen,1/))
NCERR(__LINE__)
idex = idate(0)
idr=mm5date2numeric(datebuf)
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS FILE'
print *,'IHOURRUN=0, EXPECTED= ',idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
idr=mm5date2numeric(datebuf)
if (ignore_chck_dates == 0) then
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS FILE'
print *,'IHOURRUN=0, EXPECTED= ',idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
endif
else
print *,'*** NO CHECK ON DATES in BOUN_CONCS'
endif
! Boundary concentration increment file is already open
......@@ -204,11 +208,15 @@ NCERR(__LINE__)
bounconcincr_ncid,bounconcincr_times_varid,datebuf,(/1,1/),(/dlen,1/))
NCERR(__LINE__)
idex = idate(0)
idr=mm5date2numeric(datebuf)
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS INCREMENT FILE'
print *,'IHOURRUN=0, EXPECTED= ',idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
idr=mm5date2numeric(datebuf)
if (ignore_chck_dates == 0) then
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS INCREMENT FILE'
print *,'IHOURRUN=0, EXPECTED= ',idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
endif
else
print *,'*** NO CHECK ON DATES in BOUN_CONCS'
endif
!*** Lateral Boundary Concentrations
......@@ -247,11 +255,15 @@ NCERR(__LINE__)
bounconc_ncid, bounconc_times_varid, datebuf,(/1,2/),(/dlen,1/))
NCERR(__LINE__)
idex = idate(1)
idr=mm5date2numeric(datebuf)
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS FILE'
print *,'IHOURRUN=1, EXPECTED= ',idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
idr=mm5date2numeric(datebuf)
if (ignore_chck_dates == 0) then
if(idr.ne.idex) then
print *,'*** ERROR: WRONG EXPECTED DATE IN BOUN_CONCS FILE'
print *,'IHOURRUN=1, EXPECTED= ',idex,' BOUN_CONCS= ',idr
call exit1('Exiting')
endif
else
print *,'*** NO CHECK ON DATES in BOUN_CONCS'
endif
stvec4 =(/1 , 1 , 1, 2/)
cntvec4=(/nspecboun, nhbound_domain, nverti, 1/)
......
"""
Read grib files from ECMWF
Options :
1) Surface fields only, 3d field by default
2) Decumulation : in this case 2 dates are saved in the fetch list and the decumulation is next_date - date
"""
from .fetch import fetch
from .get_domain import get_domain
from .read import read
......@@ -13,3 +23,21 @@ _fullname = "ECMWF grib2 data files"
# "accepted":
# },
#}
input_arguments={
"decumul": {
"doc": "In case of cumulated variable"
"If True, decumulation done",
"default": False,
"accepted": bool
},
"surface": {
"doc": "In case of surface variable"
"If True, 2d field",
"default": False,
"accepted": bool
}
}
......@@ -5,6 +5,7 @@ import pandas as pd
from .....utils import path
from .utils import find_valid_file
from logging import info
def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
......@@ -14,12 +15,22 @@ def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
list_files = {}
for dd in list_period_dates:
dir_dd = dd.strftime(ref_dir)
dir_dd_next = (dd + datetime.timedelta(days=1)).strftime(ref_dir)
files_3d, dates_3d = find_valid_file(dir_dd, ref_file, dd, dir_dd_next)
dir_dd_next = (dd + datetime.timedelta(hours=3)).strftime(ref_dir)
dir_dd_previous = (dd - datetime.timedelta(hours=3)).strftime(ref_dir)
files_3d, dates_3d = find_valid_file(dir_dd, ref_file, dd, dir_dd_next, ref_dir_previous=dir_dd_previous)
if os.path.isfile(files_3d[0]):
list_dates[dd] = [[dd, dd + datetime.timedelta(hours=3)]]
list_files[dd] = [files_3d[0]]
if hasattr(tracer, "decumul"):
if tracer.decumul:
if list_files:
list_files[dd] = [files_3d]
else :
list_files[dd] = [files_3d]
else :
list_files[dd] = [[files_3d[0]]]
else :
list_files[dd] = [[files_3d[0]]]
# Fetching
# Renaming target files according to date in case
......@@ -28,5 +39,6 @@ def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
target_file = "{}/{}".format(target_dir, dd.strftime(ref_file))
path.link(files_3d[0], target_file)
local_files.append(target_file)
#info(list_files)
#info(list_dates)
return list_files, list_dates
......@@ -7,7 +7,6 @@ import copy
def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
date_ref = list(input_dates.values())[0][0]
dir_dd = date_ref.strftime(target_dir)
files_3d, dates_3d = find_valid_file(dir_dd, ref_file, date_ref, dir_dd)
lon, lat = grib_file_reader(files_3d[0], ["longitude", "latitude"])
......@@ -20,10 +19,12 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
if jscan == 0:
lat2 = np.flip(lat)
lat = copy.deepcopy(lat2)
pv = grib_file_reader(files_3d[0], [], attribute="pv")
hybrid = np.array(grib_file_reader(
files_3d[0], ["hybrid"])).flatten().astype(int)
if hasattr(tracer, "surface"):
if tracer.surface==False:
pv = grib_file_reader(files_3d[0], [], attribute="pv")
hybrid = np.array(grib_file_reader(
files_3d[0], ["hybrid"])).flatten().astype(int)
lon_min = lon.min() # - (lon[1] - lon[0]) / 2
lon_max = lon.max() # + (lon[-1] - lon[-2]) / 2
......@@ -38,22 +39,27 @@ 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)
# Reconstruct alpha and beta
ecm_nlevs = int(len(pv) / 2 - 1)
sigma_a = np.empty(ecm_nlevs)
sigma_b = np.empty(ecm_nlevs)
for ii in range(ecm_nlevs - 1):
sigma_a[ii] = (pv[ecm_nlevs - ii] + pv[ecm_nlevs - ii - 1]) / 2
sigma_b[ii] = \
(pv[1 + 2 * ecm_nlevs - ii] + pv[1 + 2 * ecm_nlevs - ii - 1]) / 2
# Vertical crop if all hybrid levels are not available
sigma_a = sigma_a[ecm_nlevs - hybrid[::-1]]
sigma_b = sigma_b[ecm_nlevs - hybrid[::-1]]
ecm_nlevs = hybrid.size
# Forcing non zero top level
if sigma_a[-1] == 0:
sigma_a[-1] = sigma_a[-2] / 100
if hasattr(tracer, "surface"):
if tracer.surface:
ecm_nlevs=1
sigma_a = np.array([0])
sigma_b = np.array([1])
else:
ecm_nlevs = int(len(pv) / 2 - 1)
sigma_a = np.empty(ecm_nlevs)
sigma_b = np.empty(ecm_nlevs)
for ii in range(ecm_nlevs - 1):
sigma_a[ii] = (pv[ecm_nlevs - ii] + pv[ecm_nlevs - ii - 1]) / 2
sigma_b[ii] = \
(pv[1 + 2 * ecm_nlevs - ii] + pv[1 + 2 * ecm_nlevs - ii - 1]) / 2
# Vertical crop if all hybrid levels are not available
sigma_a = sigma_a[ecm_nlevs - hybrid[::-1]]
sigma_b = sigma_b[ecm_nlevs - hybrid[::-1]]
ecm_nlevs = hybrid.size
# Forcing non zero top level
if sigma_a[-1] == 0:
sigma_a[-1] = sigma_a[-2] / 100
# Initializes domain
setup = Setup.from_dict(
......
......@@ -2,7 +2,7 @@ import numpy as np
import xarray as xr
from .utils import grib_file_reader
from logging import info
def read(
self,
......@@ -12,27 +12,51 @@ def read(
files,
interpol_flx=False,
comp_type=None,
tracer=None,
**kwargs
):
xout = []
outdates = []
for dd, dd_file in zip(dates, files):
jscan = grib_file_reader(dd_file, [], 'jScansPositively')
spec = grib_file_reader(dd_file, varnames)
spec_cum=[]
for ddi_file in dd_file:
jscan = grib_file_reader(ddi_file, [], 'jScansPositively')
speci = grib_file_reader(ddi_file, varnames)
# Flip latitudes if jscan = 0
if jscan == 0 and tracer.surface==False :
speci = speci[:, ::-1, :]
elif jscan == 0 and tracer.surface:
speci = speci[::-1, :]
# Flip levels if not surface field
if tracer.surface==False:
speci=speci[::-1, ...]
spec_cum.append(speci)
# decumulation if decumul
if hasattr(tracer, "decumul"):
if tracer.decumul:
info("DECUMULATION DONE")
if len(spec_cum)>1:
spec = spec_cum[1]-spec_cum[0]
else:
info('impossible DECUMULATION ?')
spec = spec_cum[0]
else : spec = spec_cum[0]
else: spec = spec_cum[0]
# Flip latitudes if jscan = 0
if jscan == 0:
spec = spec[:, ::-1, :]
# Flip levels anyway
xout.append(spec[::-1, ...])
xout.append(spec)
outdates.append(dd[0])
if hasattr(tracer, "surface"):
if tracer.surface:
xout = np.array(xout)[:,np.newaxis,: ,:]
else :
xout = np.array(xout)
xmod = xr.DataArray(
np.array(xout),
xout,
coords={"time": outdates},
dims=("time", "lev", "lat", "lon"),
)
return xmod
......@@ -56,19 +56,41 @@ def grib_file_reader(filepath, varname, attribute=None):
return varout
def find_valid_file(ref_dir, file_format, dd, ref_dir_next):
def find_valid_file(ref_dir, file_format, dd, ref_dir_next, ref_dir_previous=False):
# Get all files and dates matching the file and format
list_files_orig = os.listdir(ref_dir)
# Convert ref date
ref_date = datetime.datetime.strptime(dd.strftime(file_format), file_format)
previous_date = ref_date - datetime.timedelta(hours=3)
if previous_date.month < ref_date.month and ref_dir_previous:
try : list_files_orig += os.listdir(ref_dir_previous)
except: info ("Did not find any valid GRIB files in {} "
"with format {}"
.format(ref_dir_previous, file_format))
next_date = ref_date + datetime.timedelta(hours=3)
if next_date.month>ref_date.month:
try : list_files_orig += os.listdir(ref_dir_next)
except: info ("Did not find any valid GRIB files in {} "
"with format {}"
.format(ref_dir_previous, file_format))
list_dates_cur = []
list_files_cur = []
for f in list_files_orig:
try:
if f.find('idx') < 0:
list_dates_cur.append(
datetime.datetime.strptime(f, file_format))
list_files_cur.append(f)
list_dates_cur.append(
datetime.datetime.strptime(f, file_format))
list_files_cur.append(f)
except:
continue
#info ("exception 24.grb I for {}".format(f))
try:
if f.find('idx') < 0 :
list_dates_cur.append(datetime.datetime.strptime(f.replace('24.grb','23.grb'), file_format)+datetime.timedelta(minutes=59))
list_files_cur.append(f)
except:
#info ("exception 24.grb II for {}".format(f))
continue
list_files = np.array(list_files_cur)
list_dates = np.array(list_dates_cur)
......@@ -77,58 +99,21 @@ def find_valid_file(ref_dir, file_format, dd, ref_dir_next):
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 GRIB files in {} "
"with format {}. Please check your yml file"
.format(ref_dir, file_format))
# Convert ref date
ref_date = datetime.datetime.strptime(dd.strftime(file_format), file_format)
# Compute deltas
mask = (list_dates - ref_date) <= datetime.timedelta(0)
# find nearest previous date
file_ref1 = ref_dir + list_files[mask][np.argmax(list_dates[mask])]
date_ref1 = list_dates[mask].max()
mask = (list_dates - ref_date) >= datetime.timedelta(0)
# If empty, try to look in the directory for the next month
if len(list_dates[mask]) == 0:
list_files_next_orig = os.listdir(ref_dir_next)
list_dates_next = []
list_files_next = []
for f in list_files_next_orig:
try:
list_dates_next.append(
datetime.datetime.strptime(f, file_format))
list_files_next.append(f)
except:
continue
list_files_next = np.array(list_files_cur + list_files_next)
list_dates_next = np.array(list_dates_cur + list_dates_next)
# Sorting along dates
isort = np.argsort(list_dates_next)
list_dates_next = list_dates_next[isort]
list_files_next = list_files_next[isort]
if list_files_next == []:
raise Exception("Did not find any valid GRIB files in {} "
"with format {}. Please check your yml file"
.format(ref_dir_next, file_format))
mask = (list_dates_next - ref_date) >= datetime.timedelta(0)
file_ref2 = ref_dir_next + list_files_next[mask][
np.argmin(list_dates_next[mask])]
date_ref2 = list_dates_next[mask].min()
else:
file_ref2 = ref_dir + list_files[mask][np.argmin(list_dates[mask])]
date_ref2 = list_dates[mask].min()
# find nearest next date
file_ref2 = ref_dir + list_files[mask][np.argmin(list_dates[mask])+1]
date_ref2 = list_dates[mask].min()
# Reconvert to original date
dd1 = dd + (date_ref1 - ref_date)
......
"""
Read CAMS products
Option to define pressure coordinate names depending on CAMS version
"""
from .get_domain import get_domain
from .fetch import fetch
from .read import read
......@@ -6,3 +12,14 @@ _name = "CAMS"
_version = "netcdf"
_fullname = "CAMS netcdf files"
input_arguments={
"aibi_name": {
"doc": "To choose ai/bi vertical coordinate names"
" instead of hyam and hybm; this depends on the version of the CAMS "
"product to be read",
"default": False,
"accepted": bool
},
}
......@@ -6,6 +6,7 @@ import xarray as xr
from .....utils.classes.setup import Setup
import copy
def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
# Looking for a reference file to read lon/lat in
......@@ -61,8 +62,13 @@ 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["hyam"].values
sigma_b = nc["hybm"].values
if tracer.aibi_name:
sigma_a = nc["ap"].values
sigma_b = nc["bp"].values
else :
sigma_a = nc["hyam"].values
sigma_b = nc["hybm"].values
nlevs = sigma_a.size
# Initializes domain
......
import datetime
import calendar
import os
import pandas as pd
import numpy as np
......@@ -30,16 +31,16 @@ def read(
comp_type: type of boundary conditions to generate
"""
if ddi is None:
raise Exception("CAMS netCDF read function was called "
"without specifying ddi")
# Variable name to extract
var2extract = name
if varnames != "":
var2extract = varnames
# Reading fields for periods within the simulation window
xout = []
opened_file = ""
......@@ -48,16 +49,20 @@ def read(
if dd_file != opened_file:
ds = xr.open_dataset(dd_file)
opened_file = dd_file
ntimes = ds.dims["time"]
freq = dd[0].days_in_month * 24 / ntimes
freq = pd.DatetimeIndex([dd[0]]).days_in_month[0] * 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[var2extract].values[date_index]
if lat[1] < lat[0]:
if lat[1] < lat[0] and conc.ndim == 4:
conc = conc[:, :, ::-1, :]
elif lat[1] < lat[0] and conc.ndim == 3:
conc = conc[:, ::-1, :]
xout.append(conc)
......@@ -66,5 +71,4 @@ def read(
coords={"time": np.array(dates)[:, 0]},
dims=("time", "lev", "lat", "lon"),
)
return xmod
......@@ -4,6 +4,7 @@ import os
import calendar
import numpy as np
def find_valid_file(ref_dir, file_format, dd):
# Get all files and dates matching the file and format
......
"""
Read TNO yearly fluxes and apply time profil
Time profile is considered in UTC. Time zone does not taken into account.
All the profil files are mandatory even the vertical one.
If point_sources: True, vertical profil is also applied
WARNING : Currently PS are put in the TNO grid
"""
from .fetch import fetch
from .get_domain import get_domain
from .read import read
from .write import write
_name = "TNO"
_version = "netcdf"
input_arguments={
"point_sources": {
"doc": "Point Soucre type"