Processed diagnostics needed for PS49

applications/lfric_atm/metadata/field_def_diags.xml

@@ -971,6 +971,10 @@
   <field id="processed__total_snow_acc" name="total_snow_acc" long_name="surface_snowfall_amount" standard_name="snowfall_amount" unit="kg m-2" domain_ref="face"> $DT * processed__total_snow </field>
   <field id="processed__total_prec_acc" name="total_prec_acc" long_name="surface_precipitation_amount" standard_name="precipitation_amount" unit="kg m-2" domain_ref="face"> $DT * processed__total_prec </field>
   <field id="processed__virtual_theta" name="virtual_theta" long_name="virtual_potential_temperature" standard_name="virtual_potential_temperature" unit="K" grid_ref="full_level_face_grid"> theta * (1 + 0.61 * processed__qv) </field>
+  <field id="processed__pv_on_theta_levs" name="pv_on_theta_levs" long_name="ertel_potential_vorticity_at_theta_levels" standard_name="ertel_potential_vorticity" unit="m2 s-1 K kg-1" domain_ref="face" axis_ref="theta_levels"/>
+  <field id="processed__theta_on_pv2" name="theta_on_pv2" long_name="potential_temperature_at_potential_vorticity_2_surface" unit="K" domain_ref="face"/>
+  <field id="processed__u50m" name="u50m" long_name="eastward_wind_at_50m" standard_name="eastward_wind" unit="m s-1" domain_ref="face"/>
+  <field id="processed__v50m" name="v50m" long_name="northward_wind_at_50m" standard_name="northward_wind" unit="m s-1" domain_ref="face"/>
   <!-- Pressure level diagnostic group -->
   <field id="plev__heaviside" name="heaviside_func" long_name="heaviside_function_at_pressure_levels" unit="1" domain_ref="face" axis_ref="pressure_levels"/>
   <field id="plev__temp" name="plev_temp" long_name="temperature_at_pressure_levels" unit="K" domain_ref="face" axis_ref="pressure_levels"/>

interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90

@@ -0,0 +1,119 @@
+!-------------------------------------------------------------------------------
+! (c) Crown copyright 2026 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 Claude Code Enterprise.
+!> @brief Calculate diagnostics on height levels.
+!> @details Near surface wind components, obtained by linear interpolation of
+!!          the model level winds onto a fixed height above the ground.
+
+module height_lev_diags_alg_mod
+
+  use config_mod,                  only: config_type
+  use constants_mod,               only: r_def, i_def, l_def
+  use field_collection_mod,        only: field_collection_type
+  use field_mod,                   only: field_type
+  use fs_continuity_mod,           only: W3, Wtheta
+  use initialise_diagnostics_mod,  only: init_diag => init_diagnostic_field
+  use io_config_mod,               only: use_xios_io
+  use level_interp_mdi_kernel_mod, only: interp_order_linear
+  use sci_geometric_constants_mod, only: get_height_fv
+  use timing_mod,                  only: start_timing, stop_timing, tik, LPROF
+
+  implicit none
+
+  private
+
+  public :: height_lev_diags_alg
+
+  !> Height above the ground, in metres, at which the near surface wind
+  !> diagnostics are produced.
+  real(r_def), parameter :: near_surface_wind_height = 50.0_r_def
+
+contains
+
+  !> @brief   Calculate near surface wind diagnostics on a height surface.
+  !> @details The wind components already held on W3 are interpolated onto a
+  !!          single height above ground. Both diagnostics are skipped
+  !!          entirely unless requested. The lowest W3 level lies part way up
+  !!          the first layer, so on a vertical grid coarse enough to put that
+  !!          level above the target height the diagnostic is missing data
+  !!          rather than extrapolated.
+  !> @param[in] config          Application namelist configuration object
+  !> @param[in] derived_fields  Group of derived fields
+  subroutine height_lev_diags_alg(config, derived_fields)
+
+    use height_above_surface_kernel_mod,      &
+                              only: height_above_surface_kernel_type
+    use psykal_lite_phys_mod, only: invoke_level_interp_mdi_kernel_type
+
+    implicit none
+
+    ! Arguments
+    type(config_type),           intent(in) :: config
+    type(field_collection_type), intent(in) :: derived_fields
+
+    ! Model level fields
+    type(field_type), pointer :: u_in_w3
+    type(field_type), pointer :: v_in_w3
+    type(field_type), pointer :: height_w3
+    type(field_type), pointer :: height_wth
+
+    ! Diagnostic fields
+    type(field_type) :: u50m, v50m
+
+    ! Working fields
+    type(field_type) :: height_agl
+
+    logical(l_def) :: u50m_flag, v50m_flag
+
+    real(r_def) :: hlevs(1)
+
+    integer(tik) :: id
+
+    u50m_flag = init_diag(u50m, 'processed__u50m')
+    v50m_flag = init_diag(v50m, 'processed__v50m')
+
+    if (.not. (u50m_flag .or. v50m_flag)) return
+    if (.not. use_xios_io) return
+
+    if ( LPROF ) call start_timing( id, 'diags.height_levels' )
+
+    call derived_fields%get_field('u_in_w3', u_in_w3)
+    call derived_fields%get_field('v_in_w3', v_in_w3)
+    height_w3  => get_height_fv(config, u_in_w3%get_mesh(), W3)
+    height_wth => get_height_fv(config, u_in_w3%get_mesh(), Wtheta)
+
+    ! Heights are held above mean sea level, so remove the orography before
+    ! interpolating onto a fixed height above the ground.
+    call height_w3%copy_field_properties(height_agl)
+    call invoke( height_above_surface_kernel_type(   &
+                 height_agl, height_w3, height_wth) )
+
+    hlevs(1) = near_surface_wind_height
+
+    if (u50m_flag) then
+      call invoke_level_interp_mdi_kernel_type(u_in_w3,                 &
+                                               height_agl,              &
+                                               1_i_def, hlevs, u50m,    &
+                                               interp_order_linear)
+      call u50m%write_field()
+    end if
+
+    if (v50m_flag) then
+      call invoke_level_interp_mdi_kernel_type(v_in_w3,                 &
+                                               height_agl,              &
+                                               1_i_def, hlevs, v50m,    &
+                                               interp_order_linear)
+      call v50m%write_field()
+    end if
+
+    nullify(u_in_w3, v_in_w3, height_w3, height_wth)
+
+    if ( LPROF ) call stop_timing( id, 'diags.height_levels' )
+
+  end subroutine height_lev_diags_alg
+
+end module height_lev_diags_alg_mod

interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90

@@ -0,0 +1,122 @@
+!-------------------------------------------------------------------------------
+! (c) Crown copyright 2026 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 Claude Code Enterprise.
+!> @brief Potential vorticity diagnostics on isentropic and PV surfaces.
+!> @details Calculates Ertel potential vorticity on theta surfaces, and
+!!          potential temperature on the PV = +/-2 surface. Both are vertical
+!!          interpolations of fields the model already holds, and follow the
+!!          conventions of level_interp_mdi_kernel_mod: missing data where the
+!!          column does not span the target surface, and the uppermost
+!!          crossing where the coordinate is not monotonic.
+
+module pv_surface_diags_alg_mod
+
+  use constants_mod,              only: r_def, i_def, l_def
+  use field_mod,                  only: field_type
+  use io_config_mod,              only: use_xios_io
+  use initialise_diagnostics_mod, only: init_diag => init_diagnostic_field
+  use level_interp_mdi_kernel_mod, only: interp_order_linear, &
+                                         interp_order_cubic
+  use lfric_xios_diag_mod,        only: get_axis_dimension, get_axis_values
+  use timing_mod,                 only: start_timing, stop_timing, tik, LPROF
+
+  implicit none
+
+  private
+
+  public :: pv_surface_diags_alg
+
+  !> Potential vorticity of the surface used for the dynamical tropopause,
+  !> in PV units of m2 s-1 K kg-1. The absolute value of PV is used so that
+  !> the surface is PV = +/-2, giving a tropopause in both hemispheres.
+  real(r_def), parameter :: pv_trop_surface = 2.0e-6_r_def
+
+contains
+
+  !> @brief   Calculate potential vorticity diagnostics on theta and PV
+  !!          surfaces.
+  !> @details Called once potential vorticity is available on model levels. Both
+  !!          diagnostics are skipped entirely unless requested.
+  !> @param[in] pv          Potential vorticity on model levels (W3)
+  !> @param[in] theta_in_w3 Potential temperature mapped to W3
+  subroutine pv_surface_diags_alg(pv, theta_in_w3)
+
+    use psykal_lite_phys_mod, only: invoke_level_interp_mdi_kernel_type
+
+    implicit none
+
+    ! Arguments
+    type(field_type), intent(in) :: pv
+    type(field_type), intent(in) :: theta_in_w3
+
+    ! Diagnostic fields
+    type(field_type) :: pv_on_theta_levs, theta_on_pv2
+
+    ! Working fields
+    type(field_type) :: mod_pv
+
+    logical(l_def) :: pv_on_theta_levs_flag, theta_on_pv2_flag
+
+    integer(i_def) :: nthlev
+    real(r_def), allocatable :: thlevs(:)
+    real(r_def) :: pv_levs(1)
+
+    integer(tik) :: id
+
+    pv_on_theta_levs_flag = init_diag(pv_on_theta_levs,               &
+                                      'processed__pv_on_theta_levs')
+    theta_on_pv2_flag     = init_diag(theta_on_pv2,                   &
+                                      'processed__theta_on_pv2')
+
+    if (.not. (pv_on_theta_levs_flag .or. theta_on_pv2_flag)) return
+    if (.not. use_xios_io) return
+
+    if ( LPROF ) call start_timing( id, 'diags.pv_surfaces' )
+
+    ! Potential vorticity on theta surfaces. Theta increases with height
+    ! through most of the atmosphere, so the coordinate is monotonic and cubic
+    ! interpolation can be used.
+    if (pv_on_theta_levs_flag) then
+
+      nthlev = get_axis_dimension('theta_levels')
+      allocate(thlevs(nthlev))
+      thlevs = get_axis_values('theta_levels', nthlev)
+
+      call invoke_level_interp_mdi_kernel_type(pv, theta_in_w3,     &
+                                               nthlev, thlevs,      &
+                                               pv_on_theta_levs,    &
+                                               interp_order_cubic)
+      call pv_on_theta_levs%write_field()
+
+      deallocate(thlevs)
+
+    end if
+
+    ! Theta on the PV = +/-2 surface. |PV| is far from monotonic, so linear
+    ! interpolation to the uppermost crossing is used.
+    if (theta_on_pv2_flag) then
+
+      call pv%copy_field_properties(mod_pv)
+      call invoke( setval_X(mod_pv, pv),               &
+                   inc_X_powint_n(mod_pv, 2_i_def),    &
+                   inc_X_powreal_a(mod_pv, 0.5_r_def) )
+
+      pv_levs(1) = pv_trop_surface
+
+      call invoke_level_interp_mdi_kernel_type(theta_in_w3, mod_pv,  &
+                                               1_i_def, pv_levs,     &
+                                               theta_on_pv2,         &
+                                               interp_order_linear)
+      call theta_on_pv2%write_field()
+
+    end if
+
+    if ( LPROF ) call stop_timing( id, 'diags.pv_surfaces' )
+
+  end subroutine pv_surface_diags_alg
+
+end module pv_surface_diags_alg_mod

interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90

@@ -0,0 +1,97 @@
+!-------------------------------------------------------------------------------
+! (c) Crown copyright 2026 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 Claude Code Enterprise.
+!> @brief Height of the W3 levels above the surface.
+
+module height_above_surface_kernel_mod
+
+  use argument_mod,      only: arg_type,            &
+                               GH_FIELD,            &
+                               GH_READ, GH_WRITE,   &
+                               GH_REAL, CELL_COLUMN
+  use fs_continuity_mod, only: WTHETA, W3
+  use constants_mod,     only: r_def, i_def
+  use kernel_mod,        only: kernel_type
+
+  implicit none
+
+  private
+
+  !> Kernel metadata for PSyclone
+  type, public, extends(kernel_type) :: height_above_surface_kernel_type
+    private
+    type(arg_type) :: meta_args(3) = (/                 &
+         arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3),     & ! height_agl
+         arg_type(GH_FIELD, GH_REAL, GH_READ,  W3),     & ! height_w3
+         arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA)  & ! height_wth
+         /)
+    integer :: operates_on = CELL_COLUMN
+  contains
+    procedure, nopass :: height_above_surface_code
+  end type height_above_surface_kernel_type
+
+  public :: height_above_surface_code
+
+contains
+
+  !> @details Heights are held above mean sea level, and the lowest Wtheta
+  !!          level is the surface, so subtracting it from the W3 heights
+  !!          gives the height of each W3 level above ground. The lowest W3
+  !!          level sits part way up the first layer, so the result is
+  !!          strictly positive rather than zero at the bottom of the column.
+  !> @param[in]     nlayers     The number of layers
+  !> @param[in,out] height_agl  Height of the W3 levels above the surface
+  !> @param[in]     height_w3   Height of the W3 levels above mean sea level
+  !> @param[in]     height_wth  Height of the Wtheta levels above mean sea level
+  !> @param[in]     ndf_w3      Number of degrees of freedom per cell for W3
+  !> @param[in]     undf_w3     Number of total degrees of freedom for W3
+  !> @param[in]     map_w3      Dofmap for the cell at the base of the W3 column
+  !> @param[in]     ndf_wth     Number of degrees of freedom per cell for Wtheta
+  !> @param[in]     undf_wth    Number of total degrees of freedom for Wtheta
+  !> @param[in]     map_wth     Dofmap for the cell at the base of the Wtheta
+  !!                            column
+  subroutine height_above_surface_code(nlayers,    &
+                                       height_agl, &
+                                       height_w3,  &
+                                       height_wth, &
+                                       ndf_w3,     &
+                                       undf_w3,    &
+                                       map_w3,     &
+                                       ndf_wth,    &
+                                       undf_wth,   &
+                                       map_wth)
+
+    implicit none
+
+    ! Arguments added automatically in call to kernel
+    integer(kind=i_def), intent(in) :: nlayers
+    integer(kind=i_def), intent(in) :: ndf_w3
+    integer(kind=i_def), intent(in) :: undf_w3
+    integer(kind=i_def), intent(in), dimension(ndf_w3)  :: map_w3
+    integer(kind=i_def), intent(in) :: ndf_wth
+    integer(kind=i_def), intent(in) :: undf_wth
+    integer(kind=i_def), intent(in), dimension(ndf_wth) :: map_wth
+
+    ! Arguments passed explicitly from algorithm
+    real(kind=r_def), intent(inout), dimension(undf_w3)  :: height_agl
+    real(kind=r_def), intent(in),    dimension(undf_w3)  :: height_w3
+    real(kind=r_def), intent(in),    dimension(undf_wth) :: height_wth
+
+    ! Internal variables
+    integer(kind=i_def) :: k
+    real(kind=r_def)    :: surface
+
+    ! The lowest Wtheta level is the orography
+    surface = height_wth(map_wth(1))
+
+    do k = 0, nlayers - 1
+      height_agl(map_w3(1)+k) = height_w3(map_w3(1)+k) - surface
+    end do
+
+  end subroutine height_above_surface_code
+
+end module height_above_surface_kernel_mod