.. _doc-models-chimere-diagmet:

########################################
DIAGMET meteorological pre-processing
########################################

The :doc:`diagmet transform </documentation/plugins/transforms/complex/diagmet>`
(:doc:`API reference </documentation/api/plugins/transforms/complex/diagmet>`)
replicates the historical ``diagmet.f90`` pre-processor of CHIMERE. It converts raw
ECMWF fields (read on ECMWF model levels, in an ECMWF-native grid) into the derived
meteorological variables expected by CHIMERE on its own vertical grid (:doc:`METEO.nc
<meteo>`).

This page documents the physical formulas carried out by each processing step. The
implementation is split into one module per step under
``pycif/plugins/transforms/complex/diagmet/utils/``, executed in the order below by
:func:`pycif.plugins.transforms.complex.diagmet.forward.forward`. Steps communicate
through a scratch dictionary, ``transf.diag_misc``, that is created at the start of
the forward run and discarded at the end.

diagmet has no adjoint: meteorological fields are not control variables, so
:func:`pycif.plugins.transforms.complex.diagmet.adjoint.adjoint` is a no-op.

Notations used below:

.. list-table::
    :header-rows: 1
    :widths: 30 70

    * - Symbol
      - Meaning
    * - :math:`P`
      - pressure (Pa)
    * - :math:`T`
      - temperature (K)
    * - :math:`q`
      - specific humidity (kg/kg)
    * - :math:`z`
      - altitude above ground level (m)
    * - :math:`g = 9.81`
      - gravitational acceleration (m/s²)
    * - :math:`\kappa = R/C_p = 0.2857`
      - Poisson exponent
    * - :math:`R = 287.04`
      - specific gas constant of dry air (J/kg/K)
    * - :math:`C_p = 1005`
      - specific heat of dry air at constant pressure (J/kg/K)
    * - :math:`R_v = 461.5`
      - specific gas constant of water vapour (J/kg/K)
    * - :math:`L_v = 2.45\times 10^6`
      - latent heat of vaporisation (J/kg)
    * - :math:`\kappa_v = 0.4`
      - von Kármán constant

Variable glossary (CHIMERE output names, see :doc:`meteo`):

.. list-table::
    :header-rows: 1
    :widths: 15 65 20

    * - Name
      - Meaning
      - Unit
    * - ``alti``
      - Altitude of layer top
      - m
    * - ``airm``
      - Air number density
      - molec/cm³
    * - ``oro``
      - Orography
      - m
    * - ``winz``
      - Zonal wind
      - m/s
    * - ``winm``
      - Meridional wind
      - m/s
    * - ``tem2``
      - 2 m air temperature
      - K
    * - ``usta``
      - Friction velocity
      - m/s
    * - ``aerr``
      - Aerodynamic resistance
      - s/m
    * - ``obuk``
      - Obukhov length
      - m
    * - ``wsta``
      - Convective velocity scale (:math:`w_*`)
      - m/s
    * - ``sreh``
      - Surface relative humidity
      - 0-1
    * - ``nphourm``
      - Time steps per hour
      - hour⁻¹
    * - ``dpeu``
      - Entrainment in updraught
      - kg/m²/s
    * - ``dped``
      - Entrainment in downdraught
      - kg/m²/s
    * - ``dpdu``
      - Detrainment in updraught
      - kg/m²/s
    * - ``dpdd``
      - Detrainment in downdraught
      - kg/m²/s
    * - ``kzzz``
      - Vertical turbulent diffusivity (:math:`K_{zz}`)
      - m²/s
    * - ``topc``
      - Total precipitation
      - kg/m²
    * - ``atte``
      - Cloud attenuation
      - 0-1

.. contents::
    :local:
    :depth: 1

1. Pressure, altitude, air density — ``altipres``
==================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.altipres`

Pressures are diagnosed from the hybrid sigma coefficients of the input domain
(:math:`a`, :math:`b`, at layer interfaces and mid-points) and the surface pressure
:math:`P_{surf}`:

.. math::

    P_\mathrm{mid}(l) = b_\mathrm{mid}(l)\, P_\mathrm{surf} + a_\mathrm{mid}(l)
    \qquad
    P(l) = b(l)\, P_\mathrm{surf} + a(l)

Layer-interface altitudes above the surface are obtained by vertically integrating
the hydrostatic equation (:math:`dz = -\frac{R T}{g}\, d(\ln P)`, using
:math:`R/g = 29.27` K⁻¹m for the diagmet convention):

.. math::

    z(l) = \sum_{k \le l} 29.27 \; T(k) \left[\ln P(k-1) - \ln P(k)\right]

Air number density follows from the ideal gas law (with the constant folding in
Avogadro's number and the conversion to cm⁻³):

.. math::

    n_\mathrm{air} = 7.2868\times 10^{16} \; \frac{P_\mathrm{mid}}{T}
    \quad \text{[molec/cm}^3\text{]}

Orography is recovered from the input surface geopotential :math:`\Phi_\mathrm{sfc}`
(read as ``('meteo', 'oro')``):

.. math::

    z_\mathrm{oro} = \Phi_\mathrm{sfc} / g

2. Interpolation to the surface — ``defcolumn``
================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.defcolumn`

An extra "surface" level (index 0) is prepended to every 3-D field so that the
column includes both the model's lowest full level and the diagnostic surface
state. The extrapolation weight is set by the relative altitude of the two lowest
levels:

.. math::

    w = \frac{z(1)}{z(2) - z(1)}

and every field :math:`X` (pressure, specific humidity, total water) is
linearly extrapolated to the surface as:

.. math::

    X(0) = X(1) + \big(X(1) - X(2)\big)\, w

Specific humidity and total condensed water (liquid, optionally + ice + rain,
depending on the ``cice``/``rain`` options) are floored at :math:`10^{-10}` kg/kg.
Temperature at the surface uses the ECMWF 2 m temperature :math:`T_{2m}` directly,
extrapolated consistently with the lowest-level temperature :math:`T(1)`:

.. math::

    T(0) = T_{2m} - \big(T(1) - T_{2m}\big) \frac{2}{z(0) - 2}

All downstream steps operate on this augmented, ``nlev + 1``-level column.

3. Wind rotation — ``uv_rotation``
==================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.uv_rotation`

ECMWF winds are given on a zonal/meridional (lat/lon) basis, while CHIMERE grids can
be rotated/curvilinear. The local grid-cell orientation :math:`(\cos\theta,
\sin\theta)` is derived from the geographic coordinates of cell corners
(``zlonc``, ``zlatc``):

.. math::

    \cos\theta = \frac{dx}{\sqrt{dx^2+dy^2}}
    \qquad
    \sin\theta = \frac{dy}{\sqrt{dx^2+dy^2}}

with :math:`dx` the corner-to-corner longitude difference (converted to a distance
via :math:`\cos(\mathrm{lat})`) and :math:`dy` the corresponding latitude
difference. Winds are then rotated into the model grid basis:

.. math::

    u_\mathrm{model} = \cos\theta \; u - \sin\theta \; v
    \qquad
    v_\mathrm{model} = \sin\theta \; u + \cos\theta \; v

4. Urban corrections — ``mean_z0_shf_extra_urban_temp``
========================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.mean_z0_shf_extra_urban_temp`

CHIMERE's ``LANDUSE`` file gives, for each grid cell, the fraction :math:`f_k` of 9
land-use classes (class index 4 is "urban"). The ``LANDPAR`` file gives a
per-class, per-month roughness length :math:`z_{0,k}`. Four grid-cell diagnostics
are derived and stored in ``transf.diag_misc`` for use by later steps:

.. math::

    \mathrm{awf} = \sum_{k \ne 4} f_k + (1 - \mathrm{uwinfac}) \; f_4
    \qquad \text{(urban wind-reduction factor)}

.. math::

    \mathrm{az0} = \sum_k f_k \; z_{0,k}(\mathrm{month})
    \qquad \text{(area-weighted roughness length)}

.. math::

    \mathrm{auf} = \mathrm{uflxadd} \; f_4
    \qquad \text{(additional urban heat-flux fraction)}

.. math::

    \mathrm{pm} = \mathrm{pblmin} \sum_{k \ne 4} f_k
                  + \max(\mathrm{upblmin}, \mathrm{pblmin}) \; f_4

where ``uwinfac``, ``uflxadd``, ``upblmin`` and ``pblmin`` are the transform's
:doc:`configuration options </documentation/plugins/transforms/complex/diagmet>`.
The lowest-level winds are corrected by the urban wind factor:

.. math::

    u(0) \leftarrow u(0) \; \mathrm{awf}
    \qquad
    v(0) \leftarrow v(0) \; \mathrm{awf}

5. Sensible and virtual heat fluxes — ``sv_heat_flux``
=======================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.sv_heat_flux`

Potential temperature at every level:

.. math::

    \theta = T \left(\frac{P_0}{P}\right)^\kappa, \qquad P_0 = 10^5\ \mathrm{Pa}

Surface sensible (``heat``) and latent (``humf``) kinematic heat fluxes are
obtained from the ECMWF accumulated surface fluxes
:math:`\mathrm{sshf}`, :math:`\mathrm{slhf}` (W/m²):

.. math::

    \mathrm{heat} = \frac{R}{C_p} \; T_{2m} \; \frac{\mathrm{sshf}}{P(0)}
    \qquad
    \mathrm{humf} = \frac{R}{L_v} \; T_{2m} \; \frac{\mathrm{slhf}}{P(0)}

combined into a virtual potential-temperature flux, with an optional urban
contribution ``auf`` (see step 4):

.. math::

    \mathrm{potf} = \mathrm{heat} \left(1 + 0.61\, q(0)\right)
                    + 0.61\, \mathrm{humf}\, \theta(0)
                    + \frac{R}{C_p}\, T_{2m}\, \frac{\mathrm{auf}}{P(0)}

A convective velocity scale is derived using a fixed convective height scale
:math:`h_{cs} = 1500` m and the clipped, positive flux
:math:`\mathrm{potf}_+ = \max(\mathrm{potf}, 10^{-6})`:

.. math::

    w_{*0} = \left(g\, \mathrm{potf}_+\, h_{cs} \,/\, T_{2m}\right)^{1/3}

6. Friction velocity — ``friction_velocity``
=============================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.friction_velocity`

Only computed when ``usta: recompute`` (default); otherwise :math:`u_*` is read
directly from ECMWF. Virtual potential temperature at every level:

.. math::

    \theta_v = \theta \left(1 + 0.61\, q - w\right)

with :math:`w` the total condensed water mixing ratio (step 2). The 10 m wind
speed :math:`u_{10}` is reduced by the urban wind factor
(:math:`u_{10,s} = \mathrm{awf}\, u_{10}`) and floored by a smoothing offset
:math:`u_\mathrm{off} = 0.5` m/s:

.. math::

    u^* _\mathrm{ref} = \max(u_{10,s}, u_\mathrm{off})

A bulk Richardson number is built between the surface and the reference height
:math:`z_* = 10` m using the neutral drag coefficient
:math:`C_{Dn} = \left(\kappa_v / \ln(z_*/z_0)\right)^2`
(:math:`z_0 = \mathrm{az0}`, step 4):

.. math::

    \mathrm{Ri} = \frac{g\, z_*\, \left(\theta_v(1) - \theta_v(0)\right)}
                       {\theta(0)\, {u^*_\mathrm{ref}}^2}

and the Louis (1982) stability function, with
:math:`f_m^{(0)} = 75\, C_{Dn}^2 \sqrt{z_*/z_0}`:

.. math::

    f_m = \begin{cases}
      1 - \dfrac{10\,\mathrm{Ri}}{1 + f_m^{(0)}\sqrt{-\mathrm{Ri}}}
        & \mathrm{Ri} < 0 \\[6pt]
      \dfrac{1}{1 + \dfrac{10\,\mathrm{Ri}}{\sqrt{1 + 5\,\mathrm{Ri}}}}
        & \mathrm{Ri} \ge 0
    \end{cases}

The drag coefficient and friction velocity, combining mechanical and convective
(:math:`w_{*0}`, step 5) contributions:

.. math::

    C_D = C_{Dn}\sqrt{f_m}
    \qquad
    u_* = C_D \sqrt{{u^*_\mathrm{ref}}^2 + (1.2\, w_{*0})^2}

7. Near-surface potential temperature — ``boundary_layer_height``
===================================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.boundary_layer_height`

Despite its name, this step does not compute the PBL height itself — CHIMERE
reuses the ECMWF boundary-layer height (``hght``) directly (it is copied through
unchanged as both an input and output of the transform). It instead interpolates
the virtual potential temperature :math:`\theta_v` (dry static energy proxy,
called ``po`` internally, computed in step 6) to a fixed reference height
:math:`z_\mathrm{therm} = 25` m, for later use by the Obukhov-length calculation
(step 9):

.. math::

    \theta_{v,\mathrm{therm}} = \theta_v(l_0) +
        \big(\theta_v(l_0+1) - \theta_v(l_0)\big)\;
        \frac{z_\mathrm{therm} - z(l_0)}{z(l_0+1) - z(l_0)}

where :math:`l_0` is the highest layer index with :math:`z(l_0) \le z_\mathrm{therm}`.

8. CFL check — ``checkcfl``
============================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.checkcfl`

Grid cell sizes in the zonal (:math:`\Delta x`) and meridional (:math:`\Delta y`)
directions are computed from cell-corner coordinates using the great-circle
(haversine) distance, with Earth radius :math:`R_\mathrm{earth} = 6371.03` km. The
number of CHIMERE physical time steps per hour is set so that the Courant number
stays below :math:`C_\mathrm{max} = 0.5` for the strongest horizontal winds
anywhere in the domain:

.. math::

    n_\mathrm{steps/hour} = \mathrm{round}\left(
        \frac{3600}{\min\left(
            C_\mathrm{max}\, \min\left|\dfrac{\Delta x}{u}\right|,\;
            C_\mathrm{max}\, \min\left|\dfrac{\Delta y}{v}\right|
        \right)}
    \right)

.. note::

    The vertical CFL contribution (from the vertical wind :math:`w`) is not
    implemented: CHIMERE's diagmet does not receive an explicit vertical
    velocity field from ECMWF.

9. Low-level relative humidity and cloud top — ``low_cloud_top``
===================================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.low_cloud_top`

Saturation vapour pressure (Pa) via the Tetens formula, and specific humidity at
saturation:

.. math::

    e_s(T) = 611 \exp\left(\frac{17.27\,(T - 273.15)}{T - 35.86}\right)
    \qquad
    q_\mathrm{sat} = \frac{0.622\, e_s}{P - e_s}

Relative humidity at every level: :math:`\mathrm{RH} = q / q_\mathrm{sat}`.
Surface relative humidity (``sreh``) is :math:`\mathrm{RH}` at the surface level,
clipped to :math:`[0, 1]`. A maximum relative humidity ``rhmaxx`` is also derived
below a fixed altitude :math:`z_\mathrm{cldmax} = 1000` m, clipped to
:math:`[\mathrm{crhx}, 1]` with :math:`\mathrm{crhx} = 0.90`; it is used by the
vertical-diffusivity step (step 11) to enhance mixing in cloudy layers.

10. Obukhov length — ``obukov_length``
=========================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.obukov_length`

Using the near-surface potential temperature :math:`\theta_{v,\mathrm{therm}}`
(step 7), the virtual heat flux ``potf`` (step 5), and friction velocity
:math:`u_*` (step 6), the Obukhov length is:

.. math::

    L = - \frac{\theta_{v,\mathrm{therm}}\; u_*^3}
               {\kappa_v\, g\, \mathrm{potf}}

Non-dimensional heights are built at half the lowest CHIMERE layer thickness
:math:`z_2 = 0.5\, z(1)` and at the roughness length :math:`z_0 = \mathrm{az0}`:

.. math::

    \zeta = z_2 / L \qquad \zeta_0 = z_0 / L

The aerodynamic resistance :math:`r_a` (``aerr``) follows the standard
Monin-Obukhov integrated flux-profile relationships, stable
(:math:`L \ge 0`, log-linear) or unstable (:math:`L < 0`, Businger-Dyer with
:math:`\eta = (1 - 15\zeta)^{1/4}`):

.. math::

    r_a = \begin{cases}
    \dfrac{\ln(z_2/z_0) + 4.7\,(\zeta - \zeta_0)}{\kappa_v\, u_*}
        & L \ge 0 \\[8pt]
    \dfrac{\ln(z_2/z_0) + \ln\!\dfrac{(\eta_0^2+1)(\eta_0+1)^2}
                                       {(\eta^2+1)(\eta+1)^2}
           + 2\left(\arctan\eta - \arctan\eta_0\right)}{\kappa_v\, u_*}
        & L < 0
    \end{cases}

The convective velocity scale ``wsta`` (:math:`w_*`) is non-zero only in unstable
conditions, using the PBL height :math:`h` (ECMWF ``hght``) and the flux
``potfm`` :math:`= \max(\mathrm{potf}, 10^{-6})` (step 5):

.. math::

    w_* = \begin{cases} 0 & L \ge 0 \\
      \left(g\, \mathrm{potfm}\, h \,/\, \theta_{v,\mathrm{therm}}\right)^{1/3}
        & L < 0 \end{cases}

11. Vertical turbulent diffusivity — ``vertical_turbulent_diffusivity``
==========================================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.vertical_turbulent_diffusivity`

*Inside the PBL* (:math:`z_n = z(l)/h \le 1`, with :math:`h` the PBL height): a
convective velocity scale :math:`w_c` is built from the friction/convective
velocities, and combined with an O'Brien-type parabolic vertical profile. In
stable conditions (:math:`L>0`):

.. math::

    w_c = \frac{u_*}{1 + 4.7\, z(l)/L}

in unstable conditions, with :math:`\epsilon_p = \min(0.1, z_n)`:

.. math::

    w_c = \left(u_*^3 + 2.8\, \epsilon_p\, w_*^3\right)^{1/3}

.. math::

    K_{zz} = \mathrm{clip}\left(
        \kappa_v\, w_c\, h\, z_n\,(1-z_n)^2,\;
        \left[K_\mathrm{min,BL},\, K_\mathrm{max}\right]
    \right)

with :math:`K_\mathrm{max} = 500` m²/s and :math:`K_\mathrm{min,BL}` linearly
enhanced between the dry value :math:`K_\mathrm{min,dry} = 0.1` m²/s and the
cloudy value :math:`K_\mathrm{min,wet} = 5.0` m²/s according to the maximum
relative humidity ``rhmaxx`` of step 9:

.. math::

    K_\mathrm{min,BL} = K_\mathrm{min,dry} +
        \left(K_\mathrm{min,wet} - K_\mathrm{min,dry}\right)
        \frac{\mathrm{rhmaxx} - \mathrm{crhx}}{1 - \mathrm{crhx}}

Layers straddling the PBL top are linearly blended towards the free-troposphere
minimum :math:`K_\mathrm{min,up} = 0.1` m²/s.

*Above the PBL*: a local gradient Richardson number is computed from wind shear
and virtual potential temperature gradients between levels, with a moist
correction (using :math:`\alpha = L_v q / (R T)`,
:math:`\chi = L_v^2 q / (C_p R_v T^2)`) applied in layers where relative humidity
exceeds :math:`\mathrm{crhx} = 0.90`:

.. math::

    \mathrm{Ri} = \frac{g\, \Delta\theta_v}{\Delta z \, S^2\, \theta(0)}
    \, ,\qquad
    S^2 = \left(\frac{\Delta u}{\Delta z}\right)^2
        + \left(\frac{\Delta v}{\Delta z}\right)^2 + 10^{-6}

A mixing-length diffusivity scale, with upper-air mixing length
:math:`\lambda = 150` m:

.. math::

    D_k = \frac{S}{\left(\dfrac{1}{\kappa_v z} + \dfrac{1}{\lambda}\right)^2}

and the Louis (1982) free-atmosphere stability correction:

.. math::

    K_{zz} = \mathrm{clip}\left(
    D_k \times \begin{cases}
      1 - \dfrac{8\,\mathrm{Ri}}{1 + 1.286\sqrt{-\mathrm{Ri}}} & \mathrm{Ri} < 0 \\[6pt]
      \dfrac{1}{(1 + 5\,\mathrm{Ri})^2} & \mathrm{Ri} \ge 0
    \end{cases}
    ,\;\left[K_\mathrm{min,up},\, K_\mathrm{max}\right]\right)

12. Convective mass fluxes — ``convection``
==============================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.convection`

When ``deep_convection`` is enabled, ECMWF provides updraught mass flux/detrainment
(:math:`M_u`, :math:`D_u`) and downdraught mass flux/detrainment
(:math:`M_d`, :math:`D_d`) profiles; otherwise all four are set to zero.
Entrainment is recovered as the level-to-level convergence of the mass flux, plus
the detrainment term:

.. math::

    E_u(l) = \begin{cases}
        M_u(1) + D_u(1) & l = 1 \\
        M_u(l) - M_u(l-1) + D_u(l) & l > 1
    \end{cases}

.. math::

    E_d(l) = \begin{cases}
        M_d(l+1) - M_d(l) + D_d(l) & l < n_\mathrm{lev} \\
        -M_d(n_\mathrm{lev}) + D_d(n_\mathrm{lev}) & l = n_\mathrm{lev}
    \end{cases}

Outputs ``dpdu``, ``dpdd``, ``dpeu``, ``dped`` are :math:`D_u`, :math:`D_d`,
:math:`E_u`, :math:`E_d`, each clipped to non-negative values.

13. Precipitation — ``precipitations``
=========================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.precipitations`

Total precipitation is simply the sum of the ECMWF convective (``copc``) and
large-scale (``lspc``) accumulated precipitation:

.. math::

    \mathrm{topc} = \mathrm{copc} + \mathrm{lspc}

14. Cloud optical thickness — ``cloud_optical_thickness``
=============================================================

:mod:`~pycif.plugins.transforms.complex.diagmet.utils.cloud_optical_thickness`

In the current implementation, low/medium/high cloud optical depths are derived
directly from the ECMWF cloud-cover fractions (``clol``, ``clom``, ``cloh``,
each in :math:`[0,1]`) using fixed reference optical depths for a fully-covered
sky (:math:`\mathrm{opd}_{l,0}=50`, :math:`\mathrm{opd}_{m,0}=10`,
:math:`\mathrm{opd}_{h,0}=2`):

.. math::

    \mathrm{opd}_l = \mathrm{opd}_{l,0}\; \mathrm{clol}
    \qquad
    \mathrm{opd}_m = \mathrm{opd}_{m,0}\; \mathrm{clom}
    \qquad
    \mathrm{opd}_h = \mathrm{opd}_{h,0}\; \mathrm{cloh}

and combined into a broadband cloud attenuation factor:

.. math::

    \mathrm{atte} = \exp\left(-0.11\;
        \left(\mathrm{opd}_l + \mathrm{opd}_m + \mathrm{opd}_h\right)^{0.67}\right)

.. note::

    The liquid/ice-water-path-based option (using ``clwc``/``cice`` instead of the
    cloud-cover fraction, selected by the ``clol``/``clom``/``cloh`` = 1 or 2
    options) and the relative-humidity-based diagnosis of cloud fraction from
    :func:`~pycif.plugins.transforms.complex.diagmet.utils.low_cloud_top.low_cloud_top` are not
    yet ported from the original ``diagmet.f90`` — only the "read cloud cover
    fraction" option (value ``0``) is implemented.