Source code for pycif.plugins.datastreams.fields.oldlmdz_ic.get_domain

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




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 BCs:', domain_file)

    # Define horizontal grid :reference from PYVAR for 96 x 72
    # nlon = 97
    # zlon = np.zeros(nlon)
    # for i in range(nlon): zlon[i] = -180. + i*3.75
    # nlat = 73
    # zlat = np.zeros(nlat)
    # for i in range(nlat): zlat[i] = 90.0 - i*2.5
    # zlat[0] = 89.375
    # zlat[nlat-1] = -89.375
    # print('LATLATLATLAT',zlat)
    # print('LONLONLONLON',zlon)
    # print('cccccccccccc',zlat*np.pi/180.)

    rlatu, rlonv = readnc(domain_file, ['rlatu', 'rlonv'])
    lon = rlonv[:-1]*180./np.pi
    lat = rlatu*180./np.pi
    nlon = lon.size
    nlat = lat.size
    # avoid pole issues: #ATTENTION: check for all grids, here OK for 96 x 72
    lat[nlat - 1] = -90. + (lat[0] - lat[1]) / 4.
    lat[0] = 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
    sigma_a = readnc(domain_file, ['ap'])
    sigma_b = readnc(domain_file, ['bp'])
#    print('aaaaaaaaa',sigma_a)
#    print('bbbbbbbbb',sigma_b)
    nlevs = sigma_a.size - 1

    # 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