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
[docs]
def get_domain(ref_dir, ref_file, input_interval, 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(
"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
