Source code for pycif.plugins.transforms.complex.diagmet.utils.checkcfl

import numpy as np
from ......utils.geometry.dist_matrix import dist_matrix




def checkcfl(transf, inout_datastore, ddi, mapper):
    """Check and enforce CFL stability: limit wind speeds so that the Courant number stays below 1.

    Sets the number of CHIMERE physical time steps per hour (``nphourm``) so
    that the horizontal Courant number stays below 0.5 everywhere in the
    domain, given the grid-cell sizes (great-circle distance between corners)
    and the strongest horizontal winds. The vertical CFL contribution is not
    implemented (no explicit vertical wind is available from ECMWF). See
    :doc:`/documentation/doc-models/chimere/diagmet` (section 8) for the full
    derivation.

    Args:
        transf (Plugin): diagmet transform instance.
        inout_datastore (dict): mutable datastore.
        ddi (datetime): current sub-simulation date.
        mapper (dict): transform mapper.
    """
    # Compute size of grid cells
    rearth = 6371.03e3
    domain = mapper["inputs"][("meteo", "winz")]["domain"]
    zlon = domain.zlon
    zlat = domain.zlat
    zlonc = domain.zlonc
    zlatc = domain.zlatc

    # Size in meridional direction
    radlat = np.radians(0.5 * (zlatc[:, 1:] + zlatc[:, :-1]))
    dradlon = np.diff(np.radians(0.5 * (zlonc[:, 1:] + zlonc[:, :-1])), axis=0)
    ysize = rearth * 2 * np.arcsin(np.sqrt(
        0.5 * (1 - np.sin(radlat[1:]) * np.sin(radlat[:-1])
               - np.cos(radlat[1:]) * np.cos(radlat[:-1]) * np.cos(dradlon)
               )
    ))

    # Size in meridional direction
    radlat = np.radians(0.5 * (zlatc[1:] + zlatc[:-1]))
    dradlon = np.diff(np.radians(0.5 * (zlonc[1:] + zlonc[:-1])), axis=1)
    xsize = rearth * 2 * np.arcsin(np.sqrt(
        0.5 * (1 - np.sin(radlat[:, 1:]) * np.sin(radlat[:, :-1])
               - np.cos(radlat[:, 1:]) * np.cos(radlat[:, :-1]) * np.cos(
            dradlon)
        )
    ))

    # estimation of time step using CFL criteria
    # conversion from entrainment flux to vertical speed
    # = *1.3 for air density in kg/m3
    cflmax = 0.5

    # Variables from other parameters
    alti = transf.diag_misc["alti"]
    winm = inout_datastore["outputs"][("meteo", "winm")][ddi]["spec"]
    winz = inout_datastore["outputs"][("meteo", "winz")][ddi]["spec"]
    # missing: winw = inout_datastore["outputs"][("meteo", "winw")][ddi]["spec"]
    # but where does it come from when DP is read?

    layerdepth = np.diff(alti, axis=1)

    dxoveru = cflmax * np.abs(xsize / winz).min(axis=(1, 2, 3))
    dxoverv = cflmax * np.abs(ysize / winm).min(axis=(1, 2, 3))
    # missing: dxoverw = cflmax * np.abs(layerdepth / winw).min(axis=(1, 2, 3))
    nphourm = (3600. / np.minimum(dxoveru, dxoverv)).round()  # + dxoverw

    inout_datastore["outputs"][("meteo", "nphourm")][ddi]["spec"] = nphourm