import numpy as np
import glob
import datetime
import os
import xarray as xr
from .....utils.classes.setup import Setup
from .....utils.netcdf import readnc
import copy
from .....utils.check.errclass import CifError
[docs]
def get_domain(ref_dir, ref_file, input_interval, target_dir, tracer=None):
"""Build the legacy LMDZ4 72x96x19 Domain from a reference concentration
file.
Looks for a reference NetCDF file in ``ref_dir``, either an exact
match for ``ref_file`` or the first file whose name matches
``ref_file`` used as a date-format pattern, reads ``lat``/``lon``
(flipping latitude to increasing order if needed), computes
grid-cell corners manually (with an ad hoc pole adjustment), and
hardcodes the 19-level hybrid-sigma coefficients (``sigma_a``/
``sigma_b``) when the file's ``presnivs`` dimension has size 19,
since they are not available in the input files.
.. warning::
If ``presnivs`` size differs from 19, ``sigma_a``/``sigma_b``
are never assigned in this function, which raises an
unbound-variable error rather than a domain-specific one.
Args:
ref_dir: directory to search for a reference file.
ref_file: exact file name, or a date-format pattern used to
match one of the files in ``ref_dir``.
input_interval: unused, accepted for interface compatibility.
target_dir: unused, accepted for interface compatibility.
tracer: unused, accepted for interface compatibility.
Returns:
Domain: the legacy LMDZ4 grid domain.
Raises:
CifError: if no reference file can be found in ``ref_dir``.
"""
# Looking for a reference file to read lon/lat in
list_file = glob.glob(f"{ref_dir}/*nc")
domain_file = None
# Either a file is specified in the Yaml
if ref_file in list_file:
domain_file = f"{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 CifError(
"Old LMDz domain could not be initialized as no file was found"
)
print('Domain file for old LMDZ presrcribed concs:', domain_file)
lat, lon = readnc(domain_file, ['lat', 'lon'])
# must be increasing
if lat[1] < lat[0]:
lat2 = np.flip(lat)
lat = copy.copy(lat2)
nlon = lon.size
nlat = lat.size
# avoid pole issues: #ATTENTION: check for all grids, here OK for 73 lat
lat[0] = -90. + (lat[0] - lat[1]) / 4.
lat[nlat - 1] = 90. - (lat[0] - lat[1]) / 4.
lon_min = lon.min()
lon_max = lon.max()
lat_min = lat.min()
lat_max = lat.max()
# Compute corners
lonc = np.zeros(nlon + 1)
latc = np.zeros(nlat + 1)
for i in range(nlon + 1):
if i == nlon - 1:
lonc[i] = lon[i] - (lon[i]-lon[i-1]) / 2.
elif i == nlon:
lonc[i] = lon[i-1] + (lon[i-1] - lonc[i-1])
else:
lonc[i] = lon[i] - (lon[i+1]-lon[i]) / 2.
for i in range(nlat + 1):
if i == nlat - 1:
latc[i] = lat[i] - (lat[i]-lat[i-1]) / 2.
elif i == nlat:
latc[i] = lat[i-1] + (lat[i-1] - latc[i-1])
else:
latc[i] = lat[i] - (lat[i+1]-lat[i]) / 2.
latc[latc > 90.] = 90.
latc[latc < -90.] = -90.
# Read vertical information in domain_file
# a et b pas dispos -> en dur
list_pres = readnc(domain_file, ['presnivs'])
nlevs = list_pres.size
if nlevs == 19:
sigma_a = [5.55111512e-12, 5.93108832e+02, 1.39918307e+03, 2.62124872e+03,
4.45505592e+03, 7.06566628e+03, 1.05237284e+04, 1.46678361e+04,
1.88878426e+04, 2.19670207e+04, 2.24441964e+04, 1.99314918e+04,
1.59450159e+04, 1.20320029e+04, 8.60958234e+03, 5.77252298e+03,
3.55095334e+03, 1.91474921e+03, 7.76486716e+02, 0.00000000e+00]
sigma_b = [1.00000000e+000, 9.76412025e-001, 9.44760806e-001, 8.97554243e-001,
8.28154016e-001, 7.31510416e-001, 6.05607987e-001, 4.54277335e-001,
2.91469876e-001, 1.44243468e-001, 4.52531725e-002, 6.06048280e-003,
1.47155134e-004, 8.60044395e-008, 6.12220132e-015, 7.09824791e-033,
2.12863677e-086, 1.95858259e-296, 0.00000000e+000, 0.00000000e+000]
# 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" # TO CHECK FOR LMDZ
}
}
)
Setup.load_setup(setup, level=1)
zlon, zlat = np.meshgrid(lon, lat)
zlonc, zlatc = np.meshgrid(lonc, latc)
setup.domain.zlon = zlon
setup.domain.zlat = zlat
setup.domain.zlonc = zlonc
setup.domain.zlatc = zlatc
return setup.domain
