Addition of a new diagnostic for meridional atmospheric temperature

gradient
+ fix of the documentation

documentation/source/conf.py

@@ -79,9 +79,9 @@
 # -- Extension configuration -------------------------------------------------
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {'python': ('https://docs.python.org/3', None),
-                       'numpy': ('https://docs.scipy.org/doc/numpy/', None),
-                       'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
-                       'matplotlib': ('https://matplotlib.org', None),
+                       'numpy': ('https://numpy.org/doc/stable/', None),
+                       'scipy': ('https://docs.scipy.org/doc/scipy/', None),
+                       'matplotlib': ('https://matplotlib.org/stable/', None),
                        'sparse': ('https://sparse.pydata.org/en/stable/', None),
                        'sympy': ('https://docs.sympy.org/latest/', None),
                        'ipython': ('https://ipython.readthedocs.io/en/stable/', None),

documentation/source/files/technical/diagnostics.rst

@@ -26,6 +26,10 @@ We describe here the modules used to analyze the output of the model with diagno
 .. automodule:: qgs.diagnostics.multi
     :members:
 
+.. automodule:: qgs.diagnostics.differential
+    :members:
+
+
 .. automodule:: qgs.diagnostics.wind
     :members:
 

documentation/source/files/user_guide.rst

@@ -491,6 +491,7 @@ Note that it is also possible to use other ordinary differential equations integ
 
     * :class:`.MiddleAtmosphericTemperatureAnomalyDiagnostic`: Diagnostic giving the middle atmospheric temperature anomaly fields :math:`\delta T_{\rm a}`.
     * :class:`.MiddleAtmosphericTemperatureDiagnostic`: Diagnostic giving the middle atmospheric temperature fields :math:`T_{\rm a}`.
+    * :class:`.MiddleAtmosphericTemperatureMeridionalGradientDiagnostic`: Diagnostic giving the meridional gradient of the middle atmospheric temperature fields :math:`\partial_y T_{\rm a}`.
     * :class:`.OceanicLayerTemperatureAnomalyDiagnostic`: Diagnostic giving the oceanic layer temperature anomaly fields :math:`\delta T_{\rm o}`.
     * :class:`.OceanicLayerTemperatureDiagnostic`: Diagnostic giving the oceanic layer temperature fields :math:`T_{\rm o}`.
     * :class:`.GroundTemperatureAnomalyDiagnostic`: Diagnostic giving the ground layer temperature anomaly fields :math:`\delta T_{\rm g}`.

qgs/diagnostics/differential.py

@@ -1,6 +1,6 @@
 """
-    Differential base class
-    =======================
+    Differential diagnostic base class
+    ==================================
 
     Abstract base classes defining diagnostics on differnentiated grids.
 
@@ -48,7 +48,7 @@ def __init__(self, model_params, dimensional):
 
         FieldDiagnostic.__init__(self, model_params, dimensional)
 
-    def _configure_differential(self, basis, derivative, order, delta_x=None, delta_y=None):
+    def _configure_differential_grid(self, basis, derivative, order, delta_x=None, delta_y=None):
 
         self._compute_grid(delta_x, delta_y)
 

qgs/diagnostics/temperatures.py

@@ -7,14 +7,16 @@
     Description of the classes
     --------------------------
 
-    * :class:`AtmosphericTemperatureDiagnostic`: General base class for atmospheric temperature fields diagnostic.
+    * :class:`AtmosphericTemperatureDiagnostic`: General base class for atmospheric temperature fields diagnostics.
     * :class:`MiddleAtmosphericTemperatureDiagnostic`: Diagnostic giving the middle atmospheric anomaly fields :math:`T_{\\rm a}`.
     * :class:`MiddleAtmosphericTemperatureAnomalyDiagnostic`: Diagnostic giving the middle atmospheric temperature  anomaly fields :math:`\\delta T_{\\rm a}`.
     * :class:`OceanicTemperatureDiagnostic`: General base class for oceanic temperature fields diagnostic.
     * :class:`OceanicLayerTemperatureDiagnostic`: Diagnostic giving the oceanic layer temperature fields :math:`T_{\\rm o}`.
     * :class:`OceanicLayerTemperatureAnomalyDiagnostic`: Diagnostic giving the oceanic layer temperature anomaly fields :math:`\\delta T_{\\rm o}`.
     * :class:`GroundTemperatureDiagnostic`: Diagnostic giving the ground layer temperature fields :math:`T_{\\rm g}`.
     * :class:`GroundTemperatureAnomalyDiagnostic`: Diagnostic giving the ground layer temperature anomaly fields :math:`\\delta T_{\\rm g}`.
+    * :class:`AtmosphericTemperatureMeridionalGradientDiagnostic`: General base class for meridional gradient of atmospheric temperature fields diagnostics.
+    * :class:`MiddleAtmosphericTemperatureMeridionalGradientDiagnostic`: Diagnostic giving the meridional gradient of the middle atmospheric temperature fields :math:`\\partial_y T_{\\rm a}`.
 
 """
 import warnings
@@ -24,6 +26,7 @@
 from qgs.diagnostics.util import create_grid_basis
 
 from qgs.diagnostics.base import FieldDiagnostic
+from qgs.diagnostics.differential import DifferentialFieldDiagnostic
 
 
 class AtmosphericTemperatureDiagnostic(FieldDiagnostic):
@@ -549,6 +552,136 @@ def _get_diagnostic(self, dimensional):
         return self._diagnostic_data
 
 
+class AtmosphericTemperatureMeridionalGradientDiagnostic(DifferentialFieldDiagnostic):
+    """General base class for atmospheric temperature fields meridional gradient diagnostic.
+    Provide a spatial gridded representation of the fields.
+    This is an `abstract base class`_, it must be subclassed to create new diagnostics!
+
+    .. _abstract base class: https://docs.python.org/3/glossary.html#term-abstract-base-class
+
+    Parameters
+    ----------
+
+    model_params: QgParams
+        An instance of the model parameters.
+    delta_x: float, optional
+        Spatial step in the zonal direction `x` for the gridded representation of the field.
+        If not provided, take an optimal guess based on the provided model's parameters.
+    delta_y: float, optional
+        Spatial step in the meridional direction `y` for the gridded representation of the field.
+        If not provided, take an optimal guess based on the provided model's parameters.
+    dimensional: bool
+        Indicate if the output diagnostic must be dimensionalized or not.
+
+    Attributes
+    ----------
+
+    dimensional: bool
+        Indicate if the output diagnostic must be dimensionalized or not.
+    """
+
+    def __init__(self, model_params, delta_x=None, delta_y=None, dimensional=True):
+
+        DifferentialFieldDiagnostic.__init__(self, model_params, dimensional)
+
+        self._configure(delta_x=delta_x, delta_y=delta_y)
+        self._default_plot_kwargs['cmap'] = plt.get_cmap('coolwarm')
+
+    def _compute_grid(self, delta_x=None, delta_y=None):
+
+        if delta_x is None:
+            ams = self._model_params.ablocks
+            if ams is None:
+                warnings.warn("AtmosphericTemperatureDiagnostic: Unable to configure the grid automatically. Atmospheric wavenumbers information not " +
+                              "present in the model's parameters ! Please call the compute_grid method with the delta_x and delta_y parameters.")
+                return 1
+            xwn = [ams[i][0] for i in range(len(ams))]
+            mxwn = max(xwn)
+            n_point_x = 4 * mxwn + 2
+        else:
+            n_point_x = int(np.ceil((2 * np.pi / self._model_params.scale_params.n) / delta_x) + 1)
+
+        if delta_y is None:
+            ams = self._model_params.ablocks
+            if ams is None:
+                warnings.warn("AtmosphericTemperatureDiagnostic: Unable to configure the grid automatically. Atmospheric wavenumbers information not " +
+                              "present in the model's parameters ! Please call the compute_grid method with the delta_x and delta_y parameters.")
+                return 1
+            ywn = [ams[i][1] for i in range(len(ams))]
+            mywn = max(ywn)
+            n_point_y = 4 * mywn + 2
+        else:
+            n_point_y = int(np.ceil(np.pi / delta_y) + 1)
+
+        x = np.linspace(0., 2 * np.pi / self._model_params.scale_params.n, n_point_x)
+        y = np.linspace(0., np.pi, n_point_y)
+        self._X, self._Y = np.meshgrid(x, y)
+
+    def _configure(self, delta_x=None, delta_y=None):
+
+        basis = self._model_params.atmospheric_basis
+
+        self._configure_differential_grid(basis, "dy", 1, delta_x, delta_y)
+
+        if self._orography and self._X is not None and self._Y is not None:
+            if self._model_params.ground_params.orographic_basis == "atmospheric":
+                self._oro_basis = create_grid_basis(basis, self._X, self._Y, self._subs)
+            else:
+                self._oro_basis = create_grid_basis(self._model_params.ground_basis, self._X, self._Y, self._subs)
+        else:
+            self._oro_basis = None
+
+
+class MiddleAtmosphericTemperatureMeridionalGradientDiagnostic(AtmosphericTemperatureMeridionalGradientDiagnostic):
+    """Diagnostic giving the meridional gradient of the middle atmospheric temperature fields :math:`\\partial_y T_{\\rm a}` at 500hPa.
+    It is identified with the meridional gradient of the baroclinic streamfunction :math:`\\partial_y \\theta_{\\rm a}` of the system.
+    See also :ref:`files/model/oro_model:Mid-layer equations and the thermal wind relation` sections.
+
+    Parameters
+    ----------
+
+    model_params: QgParams
+        An instance of the model parameters.
+    delta_x: float, optional
+        Spatial step in the zonal direction `x` for the gridded representation of the field.
+        If not provided, take an optimal guess based on the provided model's parameters.
+    delta_y: float, optional
+        Spatial step in the meridional direction `y` for the gridded representation of the field.
+        If not provided, take an optimal guess based on the provided model's parameters.
+    dimensional: bool, optional
+        Indicate if the output diagnostic must be dimensionalized or not.
+        Default to `True`.
+
+    Attributes
+    ----------
+
+    dimensional: bool
+        Indicate if the output diagnostic must be dimensionalized or not.
+    """
+
+    def __init__(self, model_params, delta_x=None, delta_y=None, dimensional=True):
+
+        AtmosphericTemperatureMeridionalGradientDiagnostic.__init__(self, model_params, delta_x, delta_y, dimensional)
+        vr = self._model_params.variables_range
+        self._plot_title = r'Atmospheric 500hPa Temperature Meridional Gradient'
+        self._plot_units = r" (in " + self._model_params.get_variable_units(vr[0]) + r")"
+
+        self._color_bar_format = False
+
+    def _get_diagnostic(self, dimensional):
+
+        natm = self._model_params.nmod[0]
+        theta = np.swapaxes(self._data[natm:2 * natm, ...].T @ np.swapaxes(self._grid_basis, 0, 1), 0, 1)
+
+        if dimensional:
+            self._diagnostic_data = theta * self._model_params.temperature_scaling * 2
+            self._diagnostic_data_dimensional = True
+        else:
+            self._diagnostic_data = theta
+            self._diagnostic_data_dimensional = False
+        return self._diagnostic_data
+
+
 if __name__ == '__main__':
     from qgs.params.params import QgParams
     from qgs.params.params import QgParams
@@ -565,5 +698,8 @@ def _get_diagnostic(self, dimensional):
     time, traj = integrator.get_trajectories()
     integrator.terminate()
 
-    theta = MiddleAtmosphericTemperatureDiagnostic(pars)
+    theta = MiddleAtmosphericTemperatureAnomalyDiagnostic(pars)
     theta(time, traj)
+
+    dytheta = MiddleAtmosphericTemperatureMeridionalGradientDiagnostic(pars)
+    dytheta(time, traj)

qgs/diagnostics/wind.py

@@ -104,10 +104,10 @@ def _configure(self, delta_x=None, delta_y=None):
         basis = self._model_params.atmospheric_basis
 
         if self.type == "V":
-            self._configure_differential(basis, "dx", 1, delta_x, delta_y)
+            self._configure_differential_grid(basis, "dx", 1, delta_x, delta_y)
 
         elif self.type == "U":
-            self._configure_differential(basis, "dy", 1, delta_x, delta_y)
+            self._configure_differential_grid(basis, "dy", 1, delta_x, delta_y)
 
         elif self.type is None:
             warnings.warn("AtmosphericWindDiagnostic: Basis type note specified." +

qgs/params/params.py

@@ -835,6 +835,7 @@ class QgParams(Params):
     T4: bool
         Use or not the :math:`T^4` forcing for the evolution of the temperature field if the heat exchange is activated.
 
+
     .. _Gas constant: https://en.wikipedia.org/wiki/Gas_constant
     .. _dry air: https://en.wikipedia.org/wiki/Gas_constant#Specific_gas_constant
     .. _Stefan-Boltzmann constant: https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_constant