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

import numpy as np




def altipres(transf, inout_datastore, ddi, mapper):
    """Compute mid-layer pressures, altitudes, and air density.

    Derives:

    * ``pres_mid`` (Pa) — pressure at the middle of each model layer,
      stored in ``transf.diag_misc['pres']`` for use by downstream steps.
    * ``alti`` (m) — cumulative layer thickness above the surface,
      output as ``('meteo', 'alti')``.
    * ``airm`` (molecules cm⁻³) — air number density,
      output as ``('meteo', 'airm')`` using :math:`n = 7.29 \\times 10^{16} \\cdot P/T`.
    * ``oro`` (m) — orography from geopotential (divided by g = 9.81),
      output as ``('meteo', 'oro')``.

    See :doc:`/documentation/doc-models/chimere/diagmet` (section 1) 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 (provides domain sigma coefficients
            via ``mapper["inputs"][("meteo", "winz")]["domain"]``).
    """
    # Compute pressure (Pa) "at the middle of the layers"
    domain = mapper["inputs"][("meteo", "winz")]["domain"]
    sigma_a = domain.sigma_a
    sigma_b = domain.sigma_b
    sigma_a_mid = domain.sigma_a_mid
    sigma_b_mid = domain.sigma_b_mid
    if domain.pressure_unit == "hPa":
        sigma_a = 100 * domain.sigma_a
        sigma_a_mid = 100 * domain.sigma_a_mid

    psurf = inout_datastore["inputs"][("meteo", "psurf")][ddi]["spec"].values
    
    pres_mid = sigma_b_mid[np.newaxis, :, np.newaxis, np.newaxis] * psurf \
        + sigma_a_mid[np.newaxis, :, np.newaxis, np.newaxis]
    pres = sigma_b[np.newaxis, :, np.newaxis, np.newaxis] * psurf \
        + sigma_a[np.newaxis, :, np.newaxis, np.newaxis]
    transf.diag_misc["pres"] = pres_mid

    # Temperature (K)
    temp = inout_datastore["inputs"][("meteo", "temp")][ddi]["spec"]

    # Altitudes (m) "at the interface of the layers"
    rsurg = 29.27

    # Issue with log when top of atmosphere at 0 Pa
    pres[pres == 0] = 0.001 # Pa

    dpres = np.diff(np.log(pres), axis=1)
    dalt = rsurg * temp * dpres
    alti = np.cumsum(dalt, axis=1)
    inout_datastore["outputs"][("meteo", "alti")][ddi]["spec"] = np.abs(alti)

    # Air density (molec/cm3)
    inout_datastore["outputs"][("meteo", "airm")][ddi]["spec"] = \
        7.2868e16 * pres_mid / temp  # put constant in double?

    # Geopotential to orography
    inout_datastore["outputs"][("meteo", "oro")][ddi]["spec"] = \
        inout_datastore["inputs"][("meteo", "oro")][ddi]["spec"] / 9.81