Spectral and Newtonian Nudging

applications/lfric_atm/metadata/field_def_diags.xml

@@ -45,6 +45,11 @@
   <field id="height_wth" name="height_wth" long_name="height_of_cell_interfaces_above_earth_radius" standard_name="height_above_mean_sea_level" unit="m" grid_ref="full_level_face_grid" />
   <field id="height_w2h" name="height_w2h" long_name="height_of_cell_faces_above_earth_radius" standard_name="height_above_mean_sea_level" unit="m" grid_ref="half_level_edge_grid" />
   <field id="height_w0" name="height_w0" long_name="height_of_cell_vertices_above_earth_radius" standard_name="height_above_mean_sea_level" unit="m" grid_ref="node_grid" />
+  <!-- Spectral nudging fields -->
+  <field id="u_nudging_ext_ref" name="u_nudging_ext_ref" long_name="Eastward wind component regridded from AIFS" unit="m s-1" grid_ref="nudging_ref_half_grid" />
+  <field id="v_nudging_ext_ref" name="v_nudging_ext_ref" long_name="Northward wind component regridded from AIFS" unit="m s-1" grid_ref="nudging_ref_half_grid" />
+  <field id="temperature_nudging_ext_ref" name="temperature_nudging_ext_ref" long_name="Air temperature regridded from AIFS" unit="K" grid_ref="nudging_ref_half_grid" />
+  <field id="surface_pressure_nudging_ext_ref" name="surface_pressure_nudging_ext_ref" long_name="Surface air pressure regridded from AIFS" unit="Pa" grid_ref="nudging_ref_2D_grid" />
   <field id="theta_nudging_ref" name="theta_nudging_ref" long_name="nudging_reference_air_potential_temperature" unit="K" grid_ref="nudging_ref_full_grid" />
   <field id="u_in_w3_nudging_ref" name="u_in_w3_nudging_ref" long_name="nudging_reference_eastward_wind_at_cell_centres" unit="ms-1" grid_ref="nudging_ref_half_grid" />
   <field id="v_in_w3_nudging_ref" name="v_in_w3_nudging_ref" long_name="nudging_reference_northward_wind_at_cell_centres" unit="ms-1" grid_ref="nudging_ref_half_grid" />

applications/lfric_atm/metadata/lfric_dictionary.xml

@@ -430,6 +430,5 @@
   <field id="lf_pond_depth" name="lf_pond_depth" long_name="lf_pond_depth" unit="m" domain_ref="face" axis_ref="sea_ice_categories" />
   <field id="lf_sunocean" name="lf_sunocean" long_name="lf_sunocean" unit="ms-1" domain_ref="face" />
   <field id="lf_svnocean" name="lf_svnocean" long_name="lf_svnocean" unit="ms-1" domain_ref="face" />
-
 </field_group>
 </field_definition>

documentation/source/science_guide/nudging/convolution.rst

@@ -0,0 +1,196 @@
+.. -----------------------------------------------------------------------------
+    (c) Crown copyright Met Office. All rights reserved.
+    The file LICENCE, distributed with this code, contains details of the terms
+    under which the code may be used.
+
+    Some of the content of this file has been produced with the assistance of
+    Met Office GitHub Copilot Enterprise.
+   -----------------------------------------------------------------------------
+.. _nudging_science_convolution:
+
+Filtering by Convolution
+========================
+
+The convolution method (see :ref:`nudging_science_formulation_convolution`)
+applies a scale-selective filter :math:`\mathcal{L}` to the difference
+between the model field and the reference field, so that only the large
+scales are nudged. This section explains why the filter is applied as a
+convolution in physical space, and derives the convolution kernel used on
+the sphere.
+
+Consider filtering a field :math:`f` so that only its large scales are
+retained. Writing the spectrum of :math:`f` as :math:`\mathcal{F}[f]`,
+where :math:`\mathcal{F}` denotes a transform to spectral space, the filter
+:math:`\mathcal{L}` acts as a multiplication of the spectrum by a function
+:math:`\ell` that varies with wavenumber:
+
+.. math:: :label: eq:nudging_spectral_filter
+
+   \mathcal{L}[f] = \mathcal{F}^{-1}\left[ \ell \; \mathcal{F}[f] \right] .
+
+It is vital that :math:`\ell` does not amplify any wavenumbers, which will
+cause the numerical model to go unstable. However, it is not always
+necessary for an analytic form of :math:`\ell` to be preserved
+exactly by the discretisation. For example, an idealised "top-hat" may be
+specified for :math:`\ell`, but it may be acceptable for the discretisation to
+only approximate this.
+
+Why a convolution?
+------------------
+
+The most direct way to apply :eq:`eq:nudging_spectral_filter` would be to
+transform :math:`f` to spectral space, multiply by :math:`\ell`, and
+transform back. However, for a parallel implementation in which :math:`f`
+is distributed geographically across many processors, this is very
+inefficient: the transform requires global gather and scatter
+communications and/or global summations, whose cost grows exponentially
+with the number of processors.
+
+An alternative is to note the convolution theorem: the spectrum of a
+convolved field is the product of the spectra of the field and of the
+convolving function. Defining the convolution of :math:`f` with a function
+:math:`c` as
+
+.. math:: :label: eq:nudging_convolution_def
+
+   \mathcal{C}[f, c] := \int_\Omega f(\boldsymbol{y}) \; c(\boldsymbol{y} - \boldsymbol{x}) \; \mathrm{d}\boldsymbol{y} \equiv f \ast c ,
+
+the convolution theorem gives
+:math:`\mathcal{F}\!\left[\mathcal{C}[f, c]\right] = \mathcal{F}[f] \,
+\mathcal{F}[c]`. Since the filter :eq:`eq:nudging_spectral_filter` acts as a
+multiplication in spectral space, it can therefore be applied as a
+convolution in physical space, provided the convolving function is chosen so
+that :math:`\mathcal{F}[c] = \ell`, i.e. :math:`c = \mathcal{F}^{-1}[\ell]`.
+
+The key advantage is that, if :math:`c` has compact support, the integral in
+:eq:`eq:nudging_convolution_def` reduces to an integral over a local
+neighbourhood of each point. On a distributed mesh this requires only
+communication of neighbouring data through *halo exchanges*, rather than the
+global communications demanded by a spectral transform.
+
+In practice the convolution is truncated by restricting it to a local
+stencil of finite extent. The truncated kernel is no longer exactly equal to
+:math:`\mathcal{F}^{-1}[\ell]`, but can remain a good approximation provided
+the truncation occurs where :math:`c` is already small. When the nudging of
+the spectrum does not need to be exact, this can be acceptable.
+
+The filter on the sphere
+------------------------
+
+To apply this approach on the cubed-sphere, the convolution kernel must be
+derived for the sphere. The cubed-sphere is non-orthogonal and has
+discontinuities at panel edges, which makes a separable
+longitude--latitude Fourier filter challenging. Instead, the field is
+expanded in spherical harmonics, which are appropriate for quasi-uniform
+meshes of the sphere. This use of a convolution for scale-selective
+filtering on a cubed-sphere follows the approach of Thatcher and McGregor
+(2009).
+
+The spherical harmonics, as functions of co-latitude :math:`\theta` and
+longitude :math:`\varphi`, may be written as
+
+.. math:: :label: eq:nudging_spherical_harmonic
+
+   Y_{lm}(\theta, \varphi) = N_{lm} \, e^{i m \varphi} \, P_{lm}(\cos\theta) ,
+
+where :math:`P_{lm}` are the associated Legendre polynomials and
+:math:`N_{lm}` is a normalisation chosen so that the :math:`Y_{lm}` are
+orthonormal over the sphere. A field :math:`f` is expanded as
+:math:`f(\theta,\varphi) = \sum_{l,m} a_{lm} Y_{lm}(\theta,\varphi)`.
+
+We seek the convolution kernel :math:`c` that retains only the spherical
+harmonic components with total wavenumber between :math:`k_{min}` and
+:math:`k_{max}`. Since the spherical harmonics form an orthonormal basis,
+this filter kernel can be written
+
+.. math:: :label: eq:nudging_general_kernel
+
+   c(\theta,\varphi,\theta',\varphi') =
+   \sum_{l=k_{min}}^{k_{max}} \sum_{m=-l}^{l}
+   F_{lm} \, Y_{lm}(\theta,\varphi) \, Y_{lm}^{\ast}(\theta',\varphi') .
+
+The kernel should depend only on the distance between two points, and not on
+their orientation. This requires that there is no dependence on :math:`m` or
+:math:`\varphi`, so that :math:`F_{lm} = F_{l0}`. The sum over :math:`m` can
+then be carried out using the spherical harmonic addition theorem,
+
+.. math:: :label: eq:nudging_addition_theorem
+
+   \sum_{m=-l}^{l} Y_{lm}(\theta,\varphi) \, Y_{lm}^{\ast}(\theta',\varphi')
+   = \frac{2l+1}{4\pi} \, P_l\!\left(\cos\gamma\right) ,
+
+where :math:`P_l` is the Legendre polynomial of degree :math:`l` and
+:math:`\gamma` is the great-circle (central) angle between the two points.
+Taking an idealised top-hat filter in spectral space, :math:`F_{l0} = 1`,
+the kernel becomes a weighted sum of Legendre polynomials of the central
+angle:
+
+.. math:: :label: eq:nudging_kernel
+
+   c(\gamma) = \sum_{l=k_{min}}^{k_{max}} \frac{2l+1}{4\pi} \, P_l(\cos\gamma) .
+
+The Legendre polynomials may be evaluated efficiently using the recurrence
+relation
+
+.. math:: :label: eq:nudging_legendre_recurrence
+
+   \begin{aligned}
+     P_0(\cos\gamma) &= 1 , \\
+     P_1(\cos\gamma) &= \cos\gamma , \\
+     P_l(\cos\gamma) &= \frac{2l-1}{l} \cos\gamma \, P_{l-1}(\cos\gamma)
+       - \frac{l-1}{l} \, P_{l-2}(\cos\gamma) .
+   \end{aligned}
+
+Equation :eq:`eq:nudging_kernel` is the kernel that corresponds exactly to
+the idealised top-hat spectral filter. It resembles a :math:`\mathrm{sinc}`
+function, (mainly!) decaying away from :math:`\gamma = 0`, and is evaluated at the
+points of the stencil surrounding each grid point to give the discrete
+convolution weights.
+
+Convolution Envelope
+--------------------
+
+However the series of Legendre polynomials does not strictly have compact support
+and extends over the whole sphere. An accurate convolution would therefore
+require a stencil covering the whole domain.
+
+Instead we truncate the filter to a local stencil. To avoid the spurious
+amplification of any wavenumbers, the filter is multiplied by a Gaussian
+envelope that decays to zero towards the edge of the stencil. The width of
+this envelope is a configurable parameter of the scheme, and must be chosen
+together with the stencil extent (see :ref:`nudging_user_index` for
+guidance).
+
+Thus the actual filter used is
+
+.. math:: :label: eq:nudging_enveloped_kernel
+
+   c(\gamma) = \left[ \sum_{l=k_{min}}^{k_{max}} \frac{2l+1}{4\pi}
+   P_l(\cos\gamma) \right]
+   \exp\left[ -\frac{1}{2} \left( \frac{\gamma}{\sigma} \right)^2 \right] ,
+
+where :math:`\gamma` is the great-circle angle from the central point and
+:math:`\sigma` is the width of the Gaussian envelope.
+
+:numref:`fig_nudging_envelope` illustrates the effect of the envelope. The
+underlying Legendre filter :eq:`eq:nudging_kernel` oscillates and decays only
+slowly with distance, so that it does not have compact support. Multiplying
+by the Gaussian envelope gives the enveloped filter, which decays smoothly to
+zero well within the stencil, allowing the convolution to be truncated to a
+local region with minimal error, provided the stencil is wide enough for the
+envelope to have decayed sufficiently by its edge.
+
+.. _fig_nudging_envelope:
+
+.. figure:: legendre_envelope_convolution_filter.svg
+   :width: 90%
+   :align: center
+
+   The spherical-harmonic convolution filter for a top-hat retaining total
+   wavenumbers from :math:`k_{min} = 2` to :math:`k_{max} = 20`, plotted
+   against great-circle distance. The underlying Legendre filter
+   (black dashed, :eq:`eq:nudging_kernel`) oscillates and decays only slowly,
+   whereas the enveloped filter (blue, :eq:`eq:nudging_enveloped_kernel`)
+   decays smoothly to zero. The red circles show the discrete filter values
+   sampled at the stencil points for a stencil of depth 12 on a C56 mesh.
+

documentation/source/science_guide/nudging/formulation.rst

@@ -0,0 +1,104 @@
+.. -----------------------------------------------------------------------------
+    (c) Crown copyright Met Office. All rights reserved.
+    The file LICENCE, distributed with this code, contains details of the terms
+    under which the code may be used.
+
+    Some of the content of this file has been produced with the assistance of
+    Met Office GitHub Copilot Enterprise.
+   -----------------------------------------------------------------------------
+.. _nudging_science_formulation:
+
+Formulation
+===========
+
+Let :math:`X` denote a field being nudged (potential temperature
+:math:`\theta`, or one of the horizontal wind components :math:`u`, :math:`v`),
+and let :math:`X_{ref}` denote the corresponding reference field, having
+already been vertically interpolated onto the model's levels (see
+:ref:`nudging_science_vertical_treatment`). At each call to the nudging
+scheme, an increment :math:`\Delta X` is calculated and added to the
+model's tendency for :math:`X`.
+
+Two methods are available for calculating :math:`\Delta X`.
+
+.. _nudging_science_formulation_newtonian:
+
+Newtonian relaxation
+---------------------
+
+In Newtonian relaxation, the increment is simply proportional to the local
+difference between the field and its reference value:
+
+.. math:: :label: eq:nudging_newtonian
+
+   \Delta X = w \left( X_{ref} - X \right)
+
+where :math:`w` is a weight described in :ref:`nudging_science_weights`
+below. This relaxes :math:`X` directly towards :math:`X_{ref}` at every
+horizontal scale present in the reference data, with no spectral
+selectivity.
+
+.. _nudging_science_formulation_convolution:
+
+Convolution
+-----------
+
+In the convolution method, the local difference between the field and its
+reference value is first filtered so that only the large scales are
+retained, before being scaled by the same weight :math:`w`:
+
+.. math:: :label: eq:nudging_convolution
+
+   \Delta X = w \, \mathcal{L}\left[ X_{ref} - X \right]
+
+where :math:`\mathcal{L}` is a scale-selective filter that removes the
+short wavelengths. Rather than performing this filtering with a global
+spectral transform, :math:`\mathcal{L}` is applied as a convolution in
+physical space. The motivation for this choice, and the derivation of the
+convolution kernel on the sphere, are described in
+:ref:`nudging_science_convolution`.
+
+.. _nudging_science_weights:
+
+Nudging weights
+----------------
+
+In both of the above methods, the increment is scaled by a single weight
+:math:`w`, which combines a vertical tapering factor, a spin-up ramp, and a
+relaxation-timescale factor:
+
+.. math:: :label: eq:nudging_weight
+
+   w = w_{taper} \; r_{spinup} \; \min\left(1, \frac{\Delta t}{\tau}\right)
+
+where:
+
+* :math:`w_{taper}` is a height-dependent factor, computed as described in
+  :ref:`nudging_science_vertical_treatment`, that ramps the nudging on and
+  off over a chosen range of vertical levels;
+* :math:`r_{spinup}` is a factor that ramps linearly from 0 to 1 over a
+  configured spin-up period, described below;
+* :math:`\Delta t` is the model timestep and :math:`\tau` is a configured
+  relaxation timescale. If :math:`\tau` is set to zero, this factor is
+  instead set to 1, so that the field is relaxed fully towards the
+  reference state within a single timestep.
+
+Spin-up ramp
+~~~~~~~~~~~~
+
+The nudging weight can be ramped up gradually over a configurable spin-up
+period :math:`[t_{start}, t_{end}]`:, allowing forecasts over short time scales
+to be unperturbed.
+
+.. math:: :label: eq:nudging_spinup
+
+   r_{spinup} =
+   \begin{cases}
+     0 & t \le t_{start} \\
+     \dfrac{t - t_{start}}{t_{end} - t_{start}} & t_{start} < t < t_{end} \\
+     1 & t \ge t_{end}
+   \end{cases}
+
+where :math:`t` is the elapsed model time since the start of the run. If
+:math:`t_{start}` and :math:`t_{end}` are equal, the ramp is skipped and
+:math:`r_{spinup} = 1` once :math:`t` exceeds this value.

documentation/source/science_guide/nudging/index.rst

@@ -0,0 +1,39 @@
+.. -----------------------------------------------------------------------------
+    (c) Crown copyright Met Office. All rights reserved.
+    The file LICENCE, distributed with this code, contains details of the terms
+    under which the code may be used.
+
+    Some of the content of this file has been produced with the assistance of
+    Met Office GitHub Copilot Enterprise.
+   -----------------------------------------------------------------------------
+.. _nudging_science_index:
+
+Nudging
+================
+
+Nudging is a technique for relaxing a free-running simulation towards an
+external reference state. It can be targeted at a chosen
+range of horizontal scales and a chosen range of vertical levels. This
+allows the large-scale state of the simulation to be kept close to the
+reference state, while still permitting the model to generate its own
+small-scale variability. Two approaches are available for this nudging:
+selecting only certain horizontal scales (known as spectral nudging), or
+relaxing the model state directly towards the reference state at every
+horizontal scale (known as Newtonian relaxation).
+
+This section describes the scientific formulation of nudging in GungHo. It
+begins with an overview of the method and the relaxation increment, before
+describing how the large scales can be selected using a convolution, how this
+may be computed efficiently on a coarse mesh, and finally how the reference
+data is prepared in the vertical. For details of how to configure and
+enable nudging, see the :ref:`nudging_user_index` section of the User
+Guide.
+
+.. toctree::
+    :maxdepth: 1
+
+    overview
+    formulation
+    convolution
+    multiresolution
+    vertical_treatment