DIAGMET meteorological pre-processing

The diagmet transform (API reference) 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 (METEO.nc).

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 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 pycif.plugins.transforms.complex.diagmet.adjoint.adjoint() is a no-op.

Notations used below:

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

Variable glossary (CHIMERE output names, see Meteorological input):

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 (\(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 (\(K_{zz}\))	m²/s
topc	Total precipitation	kg/m²
atte	Cloud attenuation	0-1

1. Pressure, altitude, air density — altipres

altipres

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

\[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 (\(dz = -\frac{R T}{g}\, d(\ln P)\), using \(R/g = 29.27\) K⁻¹m for the diagmet convention):

\[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⁻³):

\[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 \(\Phi_\mathrm{sfc}\) (read as ('meteo', 'oro')):

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

2. Interpolation to the surface — defcolumn

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:

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

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

\[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 \(10^{-10}\) kg/kg. Temperature at the surface uses the ECMWF 2 m temperature \(T_{2m}\) directly, extrapolated consistently with the lowest-level temperature \(T(1)\):

\[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

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 \((\cos\theta, \sin\theta)\) is derived from the geographic coordinates of cell corners (zlonc, zlatc):

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

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

\[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

mean_z0_shf_extra_urban_temp

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

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

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

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

\[\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 configuration options. The lowest-level winds are corrected by the urban wind factor:

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

5. Sensible and virtual heat fluxes — sv_heat_flux

sv_heat_flux

Potential temperature at every level:

\[\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 \(\mathrm{sshf}\), \(\mathrm{slhf}\) (W/m²):

\[\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):

\[\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 \(h_{cs} = 1500\) m and the clipped, positive flux \(\mathrm{potf}_+ = \max(\mathrm{potf}, 10^{-6})\):

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

6. Friction velocity — friction_velocity

friction_velocity

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

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

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

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

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

\[\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 \(f_m^{(0)} = 75\, C_{Dn}^2 \sqrt{z_*/z_0}\):

\[\begin{split}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}\end{split}\]

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

\[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

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 \(\theta_v\) (dry static energy proxy, called po internally, computed in step 6) to a fixed reference height \(z_\mathrm{therm} = 25\) m, for later use by the Obukhov-length calculation (step 9):

\[\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 \(l_0\) is the highest layer index with \(z(l_0) \le z_\mathrm{therm}\).

8. CFL check — checkcfl

checkcfl

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

\[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 \(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

low_cloud_top

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

\[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: \(\mathrm{RH} = q / q_\mathrm{sat}\). Surface relative humidity (sreh) is \(\mathrm{RH}\) at the surface level, clipped to \([0, 1]\). A maximum relative humidity rhmaxx is also derived below a fixed altitude \(z_\mathrm{cldmax} = 1000\) m, clipped to \([\mathrm{crhx}, 1]\) with \(\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

obukov_length

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

\[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 \(z_2 = 0.5\, z(1)\) and at the roughness length \(z_0 = \mathrm{az0}\):

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

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

\[\begin{split}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}\end{split}\]

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

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

11. Vertical turbulent diffusivity — vertical_turbulent_diffusivity

vertical_turbulent_diffusivity

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

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

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

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

\[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 \(K_\mathrm{max} = 500\) m²/s and \(K_\mathrm{min,BL}\) linearly enhanced between the dry value \(K_\mathrm{min,dry} = 0.1\) m²/s and the cloudy value \(K_\mathrm{min,wet} = 5.0\) m²/s according to the maximum relative humidity rhmaxx of step 9:

\[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 \(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 \(\alpha = L_v q / (R T)\), \(\chi = L_v^2 q / (C_p R_v T^2)\)) applied in layers where relative humidity exceeds \(\mathrm{crhx} = 0.90\):

\[\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 \(\lambda = 150\) m:

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

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

\[\begin{split}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)\end{split}\]

12. Convective mass fluxes — convection

convection

When deep_convection is enabled, ECMWF provides updraught mass flux/detrainment (\(M_u\), \(D_u\)) and downdraught mass flux/detrainment (\(M_d\), \(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:

\[\begin{split}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}\end{split}\]

\[\begin{split}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}\end{split}\]

Outputs dpdu, dpdd, dpeu, dped are \(D_u\), \(D_d\), \(E_u\), \(E_d\), each clipped to non-negative values.

13. Precipitation — precipitations

precipitations

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

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

14. Cloud optical thickness — cloud_optical_thickness

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 \([0,1]\)) using fixed reference optical depths for a fully-covered sky (\(\mathrm{opd}_{l,0}=50\), \(\mathrm{opd}_{m,0}=10\), \(\mathrm{opd}_{h,0}=2\)):

\[\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:

\[\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 low_cloud_top() are not yet ported from the original diagmet.f90 — only the “read cloud cover fraction” option (value 0) is implemented.