From fe835800de5ea76563fbf45860a23f945a38f7d6 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Mon, 27 Jul 2026 12:09:17 +0000 Subject: [PATCH 01/10] Add potential vorticity diagnostics on theta and PV surfaces Some of the content of this change has been produced with the assistance of Anthropic Claude Opus 5 (Claude Code). Adds two diagnostics: thlev__pv Ertel potential vorticity on theta surfaces dyn_trop__theta potential temperature on the PV = +/-2 surface Both are vertical interpolations of fields the model already holds, so no new science is introduced: gungho already computes Ertel PV on model levels (W3) via potential_vorticity_diagnostic_alg, and theta_in_w3 is already a derived field. level_interp_mdi_kernel_mod interpolates a field onto surfaces of constant value of a second field. It differs from the existing pres_interp_kernel_mod in three ways: * the bracketing level is the *highest* level at which the coordinate field is at or below the target, so for a non-monotonic coordinate such as |PV| the uppermost crossing is selected. This is what makes dyn_trop__theta a dynamical tropopause rather than a boundary layer artefact; * columns which do not span the target are set to missing data rather than extrapolated; * cubic interpolation degrades to linear at both ends of the column. Cubic interpolation is used for the theta surfaces, where the coordinate is monotonic, and linear for the PV surface, taken against |PV| = 2.0e-6 m2 s-1 K kg-1. The kernel takes an array of target levels, which PSyclone cannot pass (stfc/PSyclone#1312), so the PSy layer is hand written in psykal_lite_phys_mod alongside invoke_pres_interp_kernel_type. The theta_levels axis holds a single 320 K surface, following the TIGGE standard. Both fields are added to the oper_nwp_gl output stream. Co-Authored-By: Claude Opus 5 --- .../lfric_atm/metadata/field_def_diags.xml | 2 + .../algorithm/pv_surface_diags_alg_mod.x90 | 119 +++++++ .../kernel/level_interp_mdi_kernel_mod.F90 | 201 ++++++++++++ .../source/psy/psykal_lite_phys_mod.F90 | 79 +++++ .../level_interp_mdi_kernel_mod_test.pf | 299 ++++++++++++++++++ .../app/lfric_atm/file/axis_def_plev.xml | 5 +- .../file/file_def_diags_oper_nwp_gl.xml | 2 + .../diagnostics/diagnostics_calc_mod.F90 | 41 ++- .../driver/gungho_diagnostics_driver_mod.F90 | 5 +- 9 files changed, 740 insertions(+), 13 deletions(-) create mode 100644 interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 create mode 100644 interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 create mode 100644 interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 9b01795782..18089f65d6 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -992,6 +992,8 @@ + + plev__u_clim * plev__u_clim plev__u_clim * plev__v_clim plev__u_clim * plev__omega_clim diff --git a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 new file mode 100644 index 0000000000..d33edad670 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_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 +! Anthropic Claude Opus 5 (Claude Code). +!> @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) :: thlev_pv, theta_at_pv2 + + ! Working fields + type(field_type) :: mod_pv + + logical(l_def) :: thlev_pv_flag, theta_at_pv2_flag + + integer(i_def) :: nthlev + real(r_def), allocatable :: thlevs(:) + real(r_def) :: pv_levs(1) + + integer(tik) :: id + + thlev_pv_flag = init_diag(thlev_pv, 'thlev__pv') + theta_at_pv2_flag = init_diag(theta_at_pv2, 'dyn_trop__theta') + + if (.not. (thlev_pv_flag .or. theta_at_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 (thlev_pv_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, & + thlev_pv, & + interp_order_cubic) + call thlev_pv%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_at_pv2_flag) then + + call pv%copy_field_properties(mod_pv) + call invoke( X_times_Y(mod_pv, pv, pv), & + 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_at_pv2, & + interp_order_linear) + call theta_at_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 diff --git a/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 new file mode 100644 index 0000000000..c9b51b031c --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 @@ -0,0 +1,201 @@ +!------------------------------------------------------------------------------- +! (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 +! Anthropic Claude Opus 5 (Claude Code). +!> @brief Interpolate a field onto surfaces of constant value of another field. + +module level_interp_mdi_kernel_mod + + use argument_mod, only: arg_type, & + GH_FIELD, GH_SCALAR, & + GH_READ, GH_WRITE, & + GH_INTEGER, & + GH_REAL, CELL_COLUMN, & + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + use constants_mod, only: r_def, i_def, rmdi + use kernel_mod, only: kernel_type + + implicit none + + private + + !> Interpolation orders, matching the UM's interpor_mod + integer(kind=i_def), public, parameter :: interp_order_linear = 1_i_def + integer(kind=i_def), public, parameter :: interp_order_cubic = 3_i_def + + !> Kernel metadata for PSyclone + type, public, extends(kernel_type) :: level_interp_mdi_kernel_type + private + type(arg_type) :: meta_args(5) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_SCALAR,GH_INTEGER, GH_READ), & +! arg_type(GH_SCALAR_ARRAY,GH_REAL, GH_READ, 1), see PSyclone issue #1312 + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_SCALAR,GH_INTEGER, GH_READ) & + /) + integer :: operates_on = CELL_COLUMN + contains + procedure, nopass :: level_interp_mdi_code + end type level_interp_mdi_kernel_type + + public :: level_interp_mdi_code + +contains + + !> @details Interpolates a field onto surfaces of constant value of a second + !! "coordinate" field, for example potential vorticity onto theta + !! surfaces, or potential temperature onto a potential vorticity + !! surface. This is a port of the UM's vert_interp_mdi2 and shares + !! its conventions: + !! + !! * The bracketing level is the *highest* level at which the + !! coordinate field is less than or equal to the target value. For + !! a non-monotonic coordinate such as |PV| this selects the + !! uppermost crossing, which is what makes the theta on PV=+/-2 + !! diagnostic a dynamical tropopause rather than a boundary layer + !! artefact. + !! * Columns which do not span the target value are set to missing + !! data rather than extrapolated. This differs from + !! pres_interp_kernel_mod, which clamps to the end values. + !! * Cubic interpolation degrades to linear where the four point + !! stencil would run off either end of the column. + !> @param[in] nlayers The number of layers + !> @param[in] data_in Input data to interpolate + !> @param[in] coord_in Coordinate field at levels of the input data + !> @param[in] nlev_out Number of output surfaces + !> @param[in] target_levs Coordinate values of the output surfaces + !> @param[in,out] data_out Interpolated data + !> @param[in] interp_order Interpolation order, linear or cubic + !> @param[in] ndf_in Number of degrees of freedom per cell for in fields + !> @param[in] undf_in Number of total degrees of freedom for in fields + !> @param[in] map_in Dofmap for the cell at the base of the column for in fields + !> @param[in] ndf_out Number of degrees of freedom per cell for out fields + !> @param[in] undf_out Number of total degrees of freedom for out fields + !> @param[in] map_out Dofmap for the cell at the base of the column for out fields + subroutine level_interp_mdi_code(nlayers, & + data_in, & + coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order, & + ndf_in, undf_in, map_in, & + ndf_out, undf_out, map_out) + + implicit none + + ! Arguments added automatically in call to kernel + integer(kind=i_def), intent(in) :: nlayers, nlev_out + integer(kind=i_def), intent(in) :: ndf_out, undf_out + integer(kind=i_def), intent(in), dimension(ndf_out) :: map_out + integer(kind=i_def), intent(in) :: ndf_in, undf_in + integer(kind=i_def), intent(in), dimension(ndf_in) :: map_in + + ! Arguments passed explicitly from algorithm + real(kind=r_def), intent(in), dimension(undf_in) :: data_in + real(kind=r_def), intent(in), dimension(undf_in) :: coord_in + real(kind=r_def), intent(inout), dimension(undf_out) :: data_out + + ! Constants passed explicitly from algorithm + integer(kind=i_def), intent(in) :: interp_order + real(kind=r_def), intent(in), dimension(nlev_out) :: target_levs + + ! Internal variables + integer(kind=i_def) :: k, kl, level_below, top_df, base_in, base_out + real(kind=r_def) :: target_lev + real(kind=r_def) :: c_minus, c_here, c_plus, c_plus2 + + base_in = map_in(1) + base_out = map_out(1) + + ! For Wtheta ndf = 2, loop k = 0, nlayers + ! For W3 ndf = 1, loop k = 0, nlayers - 1 + top_df = nlayers + ndf_in - 2 + + do kl = 1, nlev_out + + target_lev = target_levs(kl) + + ! Highest level whose coordinate value is at or below the target, leaving + ! room for the level above. level_below = -1 flags missing data. + level_below = -1_i_def + do k = 0, top_df - 1 + if ( coord_in(base_in+k) <= target_lev ) then + level_below = k + end if + end do + + ! Target lies above the top of the column + if ( target_lev > coord_in(base_in+top_df) ) then + level_below = -1_i_def + end if + + if ( level_below == -1_i_def ) then + + ! Column does not span the target surface + data_out(base_out+kl-1) = rmdi + + else if ( interp_order == interp_order_linear .or. & + level_below == 0_i_def .or. & + level_below == top_df - 1_i_def ) then + + ! Linear interpolation. Also used at both ends of the column, where the + ! cubic stencil would be out of range. + data_out(base_out+kl-1) = & + ( (target_lev - coord_in(base_in+level_below)) & + * data_in(base_in+level_below+1) & + - (target_lev - coord_in(base_in+level_below+1)) & + * data_in(base_in+level_below) ) & + / ( coord_in(base_in+level_below+1) & + - coord_in(base_in+level_below) ) + + else + + ! Cubic interpolation over levels level_below-1 to level_below+2 + c_minus = coord_in(base_in+level_below-1) + c_here = coord_in(base_in+level_below) + c_plus = coord_in(base_in+level_below+1) + c_plus2 = coord_in(base_in+level_below+2) + + data_out(base_out+kl-1) = & + ( (target_lev - c_here) & + * (target_lev - c_plus) & + * (target_lev - c_plus2) ) & + / ( (c_minus - c_here) & + * (c_minus - c_plus) & + * (c_minus - c_plus2) ) & + * data_in(base_in+level_below-1) & + + ( (target_lev - c_minus) & + * (target_lev - c_plus) & + * (target_lev - c_plus2) ) & + / ( (c_here - c_minus) & + * (c_here - c_plus) & + * (c_here - c_plus2) ) & + * data_in(base_in+level_below) & + + ( (target_lev - c_minus) & + * (target_lev - c_here) & + * (target_lev - c_plus2) ) & + / ( (c_plus - c_minus) & + * (c_plus - c_here) & + * (c_plus - c_plus2) ) & + * data_in(base_in+level_below+1) & + + ( (target_lev - c_minus) & + * (target_lev - c_here) & + * (target_lev - c_plus) ) & + / ( (c_plus2 - c_minus) & + * (c_plus2 - c_here) & + * (c_plus2 - c_plus) ) & + * data_in(base_in+level_below+2) + + end if + + end do + + end subroutine level_interp_mdi_code + +end module level_interp_mdi_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 index a064523c8d..240c6226cb 100644 --- a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 +++ b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 @@ -3,6 +3,8 @@ ! 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 +! Anthropic Claude Opus 5 (Claude Code). !> @brief Provides an implementation of the Psy layer for physics !> @details Contains hand-rolled versions of the Psy layer that can be used for @@ -984,6 +986,83 @@ SUBROUTINE invoke_pres_interp_kernel_type(u_in_w3, exner_w3, nplev, plevs, plev_ ! END SUBROUTINE invoke_pres_interp_kernel_type + !--------------------------------------------------------------------- + !> Contains the PSy-layer to interpolate onto surfaces of constant value + !> of a second field. As above this requires passing an array + !> "target_levs" into the kernel which is currently unsupported by + !> PSyclone, see https://github.com/stfc/PSyclone/issues/1312 + SUBROUTINE invoke_level_interp_mdi_kernel_type(data_in, coord_in, nlev_out, target_levs, data_out, interp_order) + USE level_interp_mdi_kernel_mod, ONLY: level_interp_mdi_code + USE mesh_mod, ONLY: mesh_type + INTEGER(KIND=i_def), intent(in) :: nlev_out, interp_order + REAL(KIND=r_def), intent(in) :: target_levs(nlev_out) + TYPE(field_type), intent(in) :: data_in, coord_in, data_out + INTEGER(KIND=i_def) :: cell + INTEGER(KIND=i_def) :: loop0_start, loop0_stop + INTEGER(KIND=i_def) :: nlayers + REAL(KIND=r_def), pointer, dimension(:) :: data_out_data => null() + REAL(KIND=r_def), pointer, dimension(:) :: coord_in_data => null() + REAL(KIND=r_def), pointer, dimension(:) :: data_in_data => null() + TYPE(field_proxy_type) :: data_in_proxy, coord_in_proxy, data_out_proxy + INTEGER(KIND=i_def), pointer :: map_adspc1_data_in(:,:) => null(), map_adspc2_data_out(:,:) => null() + INTEGER(KIND=i_def) :: ndf_adspc1_data_in, undf_adspc1_data_in, ndf_adspc2_data_out, undf_adspc2_data_out + INTEGER(KIND=i_def) :: max_halo_depth_mesh + TYPE(mesh_type), pointer :: mesh => null() + ! + ! Initialise field and/or operator proxies + ! + data_in_proxy = data_in%get_proxy() + data_in_data => data_in_proxy%data + coord_in_proxy = coord_in%get_proxy() + coord_in_data => coord_in_proxy%data + data_out_proxy = data_out%get_proxy() + data_out_data => data_out_proxy%data + ! + ! Initialise number of layers + ! + nlayers = data_in_proxy%vspace%get_nlayers() + ! + ! Create a mesh object + ! + mesh => data_in_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + ! + ! Look-up dofmaps for each function space + ! + map_adspc1_data_in => data_in_proxy%vspace%get_whole_dofmap() + map_adspc2_data_out => data_out_proxy%vspace%get_whole_dofmap() + ! + ! Initialise number of DoFs for adspc1_data_in + ! + ndf_adspc1_data_in = data_in_proxy%vspace%get_ndf() + undf_adspc1_data_in = data_in_proxy%vspace%get_undf() + ! + ! Initialise number of DoFs for adspc2_data_out + ! + ndf_adspc2_data_out = data_out_proxy%vspace%get_ndf() + undf_adspc2_data_out = data_out_proxy%vspace%get_undf() + ! + ! Set-up all of the loop bounds + ! + loop0_start = 1 + loop0_stop = mesh%get_last_edge_cell() + ! + ! Call kernels and communication routines + ! + DO cell=loop0_start,loop0_stop + ! + CALL level_interp_mdi_code(nlayers, data_in_data, coord_in_data, nlev_out, target_levs, data_out_data, interp_order, & +&ndf_adspc1_data_in, undf_adspc1_data_in, map_adspc1_data_in(:,cell), ndf_adspc2_data_out, undf_adspc2_data_out, & +&map_adspc2_data_out(:,cell)) + END DO + ! + ! Set halos dirty/clean for fields modified in the above loop + ! + CALL data_out_proxy%set_dirty() + ! + ! + END SUBROUTINE invoke_level_interp_mdi_kernel_type + !--------------------------------------------------------------------- !> Contains the PSy-layer to build the pressure level diagnostics !> These require passing an array "plevs" into each kernel diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf new file mode 100644 index 0000000000..18c7291a94 --- /dev/null +++ b/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf @@ -0,0 +1,299 @@ +! ----------------------------------------------------------------------------- +! (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 +! Anthropic Claude Opus 5 (Claude Code). +!> @brief Test level_interp_mdi_kernel +!> +!> The kernel is a port of the UM's vert_interp_mdi2, so these tests pin down +!> the three behaviours that differ from pres_interp_kernel_mod: missing data +!> for columns which do not span the target surface, selection of the +!> uppermost crossing where the coordinate is not monotonic, and degradation +!> of cubic interpolation to linear at the ends of the column. + +module level_interp_mdi_kernel_mod_test + + use funit + use constants_mod, only: i_def, r_def, rmdi + use level_interp_mdi_kernel_mod, only: level_interp_mdi_code, & + interp_order_linear, & + interp_order_cubic + + implicit none + + private + public :: test_linear, test_cubic, test_missing_data, & + test_uppermost_crossing, test_cubic_edge_fallback + + @TestCase + type, extends(TestCase), public :: level_interp_mdi_test_type + private + real(kind=r_def), allocatable :: coord_in(:) + real(kind=r_def), allocatable :: data_in(:) + + contains + procedure setUp + procedure tearDown + procedure test_linear + procedure test_cubic + procedure test_missing_data + procedure test_uppermost_crossing + procedure test_cubic_edge_fallback + + end type level_interp_mdi_test_type + + ! A six level W3-like column. With ndf = 1 the kernel addresses levels + ! 0 to nlayers-1, so top_df = 5. + integer(kind=i_def), parameter :: nlayers = 6_i_def + integer(kind=i_def), parameter :: ndf_in = 1_i_def + integer(kind=i_def), parameter :: ndf_out = 1_i_def + integer(kind=i_def), parameter :: undf_in = nlayers + + integer(kind=i_def), parameter :: map_in(1) = (/ 1_i_def /) + integer(kind=i_def), parameter :: map_out(1) = (/ 1_i_def /) + + real(kind=r_def), parameter :: tolerance = 1.0e-12_r_def + +contains + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @before + subroutine setUp( this ) + implicit none + class(level_interp_mdi_test_type), intent(inout) :: this + + allocate( this%coord_in(undf_in) ) + allocate( this%data_in(undf_in) ) + + end subroutine setUp + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @after + subroutine tearDown( this ) + implicit none + class(level_interp_mdi_test_type), intent(inout) :: this + + deallocate( this%coord_in, this%data_in ) + + end subroutine tearDown + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Linear interpolation of a linear field is exact. + @test + subroutine test_linear( this ) + + implicit none + + class(level_interp_mdi_test_type), intent(inout) :: this + + integer(kind=i_def), parameter :: nlev_out = 2_i_def + + real(kind=r_def) :: data_out(nlev_out) + real(kind=r_def) :: expect(nlev_out) + real(kind=r_def) :: target_levs(nlev_out) + + ! Monotonically increasing coordinate, as theta is through most of the + ! atmosphere + this%coord_in = (/ 280.0_r_def, 290.0_r_def, 300.0_r_def, & + 310.0_r_def, 320.0_r_def, 330.0_r_def /) + this%data_in = 2.0_r_def*this%coord_in + + target_levs = (/ 295.0_r_def, 325.0_r_def /) + expect = 2.0_r_def*target_levs + + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_linear, & + ndf_in, undf_in, map_in, & + ndf_out, nlev_out, map_out ) + + @assertEqual(expect(:), data_out(:), tolerance) + + end subroutine test_linear + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Cubic interpolation of a cubic field is exact where the four point + !> stencil fits inside the column. The target of 3.5 gives level_below = 2, + !> so the stencil spans levels 1 to 4 and the cubic branch is taken. + @test + subroutine test_cubic( this ) + + implicit none + + class(level_interp_mdi_test_type), intent(inout) :: this + + integer(kind=i_def), parameter :: nlev_out = 1_i_def + + real(kind=r_def) :: data_out(nlev_out) + real(kind=r_def) :: target_levs(nlev_out) + real(kind=r_def) :: expect + integer(kind=i_def) :: k + + ! Unevenly spaced, so the test cannot be passed by a linear scheme + this%coord_in = (/ 0.0_r_def, 1.0_r_def, 3.0_r_def, & + 4.0_r_def, 7.0_r_def, 9.0_r_def /) + + do k = 1, nlayers + this%data_in(k) = cubic_fn( this%coord_in(k) ) + end do + + target_levs(1) = 3.5_r_def + expect = cubic_fn( target_levs(1) ) + + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_cubic, & + ndf_in, undf_in, map_in, & + ndf_out, nlev_out, map_out ) + + @assertEqual(expect, data_out(1), tolerance) + + end subroutine test_cubic + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Columns which do not span the target surface are set to missing data + !> rather than extrapolated, which is the key difference from + !> pres_interp_kernel_mod. + @test + subroutine test_missing_data( this ) + + implicit none + + class(level_interp_mdi_test_type), intent(inout) :: this + + integer(kind=i_def), parameter :: nlev_out = 3_i_def + + real(kind=r_def) :: data_out(nlev_out) + real(kind=r_def) :: target_levs(nlev_out) + + this%coord_in = (/ 280.0_r_def, 290.0_r_def, 300.0_r_def, & + 310.0_r_def, 320.0_r_def, 330.0_r_def /) + this%data_in = 2.0_r_def*this%coord_in + + ! Below the bottom of the column, above the top of the column, and inside + target_levs = (/ 270.0_r_def, 340.0_r_def, 305.0_r_def /) + + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_linear, & + ndf_in, undf_in, map_in, & + ndf_out, nlev_out, map_out ) + + @assertEqual(rmdi, data_out(1), tolerance) + @assertEqual(rmdi, data_out(2), tolerance) + @assertEqual(610.0_r_def, data_out(3), tolerance) + + end subroutine test_missing_data + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Where the coordinate is not monotonic, as for |PV|, the uppermost + !> crossing is selected. This is what makes theta on the PV = +/-2 surface + !> a dynamical tropopause rather than a boundary layer artefact. + @test + subroutine test_uppermost_crossing( this ) + + implicit none + + class(level_interp_mdi_test_type), intent(inout) :: this + + integer(kind=i_def), parameter :: nlev_out = 1_i_def + + real(kind=r_def) :: data_out(nlev_out) + real(kind=r_def) :: target_levs(nlev_out) + + ! |PV| in PV units, crossing 2.0 three times: up through the boundary + ! layer, back down, then up again at the tropopause + this%coord_in = (/ 1.0_r_def, 3.0_r_def, 1.0_r_def, & + 1.5_r_def, 2.5_r_def, 6.0_r_def /) + ! theta, increasing with height + this%data_in = (/ 280.0_r_def, 290.0_r_def, 300.0_r_def, & + 310.0_r_def, 320.0_r_def, 330.0_r_def /) + + target_levs(1) = 2.0_r_def + + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_linear, & + ndf_in, undf_in, map_in, & + ndf_out, nlev_out, map_out ) + + ! The uppermost crossing lies between levels 3 and 4, where the + ! coordinate runs 1.5 to 2.5 and theta runs 310 to 320. Taking the lowest + ! crossing instead would give 285.0 + @assertEqual(315.0_r_def, data_out(1), tolerance) + + end subroutine test_uppermost_crossing + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Cubic interpolation degrades to linear in the bottom and top intervals, + !> where the four point stencil would run off the end of the column. + @test + subroutine test_cubic_edge_fallback( this ) + + implicit none + + class(level_interp_mdi_test_type), intent(inout) :: this + + integer(kind=i_def), parameter :: nlev_out = 2_i_def + + real(kind=r_def) :: data_out(nlev_out) + real(kind=r_def) :: target_levs(nlev_out) + + this%coord_in = (/ 0.0_r_def, 1.0_r_def, 2.0_r_def, & + 3.0_r_def, 4.0_r_def, 5.0_r_def /) + ! Strongly curved, so a cubic fit would differ markedly from a linear one + this%data_in = (/ 0.0_r_def, 1.0_r_def, 8.0_r_def, & + 27.0_r_def, 64.0_r_def, 125.0_r_def /) + + ! Bottom interval (level_below = 0) and top interval (level_below = 4) + target_levs = (/ 0.5_r_def, 4.5_r_def /) + + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_cubic, & + ndf_in, undf_in, map_in, & + ndf_out, nlev_out, map_out ) + + ! Linear midpoints, not the cubic values + @assertEqual(0.5_r_def, data_out(1), tolerance) + @assertEqual(94.5_r_def, data_out(2), tolerance) + + end subroutine test_cubic_edge_fallback + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> A cubic in the coordinate, used to check that cubic interpolation is + !> exact rather than merely close. + function cubic_fn( c ) result( f ) + + implicit none + + real(kind=r_def), intent(in) :: c + real(kind=r_def) :: f + + f = c*c*c - 2.0_r_def*c*c + 3.0_r_def*c - 5.0_r_def + + end function cubic_fn + +end module level_interp_mdi_kernel_mod_test diff --git a/rose-stem/app/lfric_atm/file/axis_def_plev.xml b/rose-stem/app/lfric_atm/file/axis_def_plev.xml index b1f446983e..ea5ca39992 100644 --- a/rose-stem/app/lfric_atm/file/axis_def_plev.xml +++ b/rose-stem/app/lfric_atm/file/axis_def_plev.xml @@ -1,6 +1,9 @@ - + + + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index fb5ee25716..02d5377c86 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -65,6 +65,8 @@ + + diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 index 9e8704c78b..d7a13dce71 100755 --- a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 +++ b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 @@ -4,6 +4,8 @@ ! under which the code may be used. !----------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). !> @brief Module for computing and outputting derived diagnostics !! !! @details Computes various derived diagnostics that are written out @@ -29,6 +31,7 @@ module diagnostics_calc_mod use initialise_diagnostics_mod, only: diagnostic_to_be_sampled, & init_diag => init_diagnostic_field use field_mod, only: field_type + use field_collection_mod, only: field_collection_type use field_parent_mod, only: write_interface use fs_continuity_mod, only: W3 use model_clock_mod, only: model_clock_type @@ -244,17 +247,20 @@ end subroutine write_vorticity_diagnostic #ifdef UM_PHYSICS !------------------------------------------------------------------------------- !> @brief Potential vorticity diagnostic processing and output. -!> @details Optionally calculate and output both model level and pressure -!> level diagnostics. -!> @param[in] u_field The wind field -!> @param[in] theta The potential temperature field (K) -!> @param[in] rho The density field (kg/m3) -!> @param[in] exner The exner pressure (Pa) -!> @param[in] clock The model clock object +!> @details Optionally calculate and output model level, pressure level, +!> theta level and PV surface diagnostics. +!> @param[in] u_field The wind field +!> @param[in] theta The potential temperature field (K) +!> @param[in] rho The density field (kg/m3) +!> @param[in] exner The exner pressure (Pa) +!> @param[in] derived_fields Group of derived fields +!> @param[in] clock The model clock object !------------------------------------------------------------------------------- -subroutine write_pv_diagnostic(u_field, theta, rho, exner, clock) +subroutine write_pv_diagnostic(u_field, theta, rho, exner, & + derived_fields, clock) - use pres_lev_diags_alg_mod, only: pres_lev_field_alg + use pres_lev_diags_alg_mod, only: pres_lev_field_alg + use pv_surface_diags_alg_mod, only: pv_surface_diags_alg implicit none @@ -262,11 +268,13 @@ subroutine write_pv_diagnostic(u_field, theta, rho, exner, clock) type(field_type), intent(in) :: theta type(field_type), intent(in) :: rho type(field_type), intent(in) :: exner + type(field_collection_type), intent(in) :: derived_fields class(model_clock_type), intent(in) :: clock type(field_type) :: pv type(field_type) :: plev_pv - logical(l_def) :: pv_modlev_flag, plev_pv_flag + type(field_type), pointer :: theta_in_w3 => null() + logical(l_def) :: pv_modlev_flag, plev_pv_flag, pv_surface_flag logical(l_def), parameter :: xi3_axis = .false. logical(l_def), parameter :: add_W3_version = .false. @@ -275,12 +283,23 @@ subroutine write_pv_diagnostic(u_field, theta, rho, exner, clock) plev_pv_flag = init_diag(plev_pv, 'plev__pv') - if (pv_modlev_flag .or. plev_pv_flag) then + ! Theta surface and PV surface diagnostics, both of which need PV on + ! model levels + pv_surface_flag = diagnostic_to_be_sampled('thlev__pv') .or. & + diagnostic_to_be_sampled('dyn_trop__theta') + + if (pv_modlev_flag .or. plev_pv_flag .or. pv_surface_flag) then call potential_vorticity_diagnostic_alg(pv, u_field, theta, rho) if (plev_pv_flag) call pres_lev_field_alg(pv, exner, plev_pv, xi3_axis) + if (pv_surface_flag) then + call derived_fields%get_field('theta_in_w3', theta_in_w3) + call pv_surface_diags_alg(pv, theta_in_w3) + nullify(theta_in_w3) + end if + if (pv_modlev_flag) then call write_scalar_diagnostic('potential_vorticity', pv, clock, & pv%get_mesh(), add_W3_version) diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index 77917aaa9a..58b1f7224c 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -4,6 +4,8 @@ ! under which the code may be used. !----------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). !> @brief Outputs diagnostics from gungho/lfric_atm !> @details Calls the routine that generates diagnostic output for @@ -243,7 +245,8 @@ subroutine gungho_diagnostics_driver( modeldb, & end if call write_vorticity_diagnostic( u, exner, modeldb%clock ) #ifdef UM_PHYSICS - call write_pv_diagnostic( u, theta, rho, exner, modeldb%clock ) + call write_pv_diagnostic( u, theta, rho, exner, derived_fields, & + modeldb%clock ) #else call write_pv_diagnostic( u, theta, rho, modeldb%clock ) #endif From e89cf9aef6669bf64a982298dc95b7221022a664 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Mon, 27 Jul 2026 13:31:34 +0000 Subject: [PATCH 02/10] Rename PV surface diagnostics and extend the theta level set Some of the content of this change has been produced with the assistance of Anthropic Claude Opus 5 (Claude Code). Renames the two diagnostics into the processed__ namespace: thlev__pv -> processed__pv_on_theta_levs dyn_trop__theta -> processed__theta_on_pv2 The XIOS output names and the local field names in pv_surface_diags_alg_mod follow suit. The theta_levels axis now carries the 475, 425 and 315 K surfaces rather than a single 320 K surface. The axis values are listed in the order requested; the kernel treats each target level independently, so the ordering affects only the layout of the output. Co-Authored-By: Claude Opus 5 --- .../lfric_atm/metadata/field_def_diags.xml | 4 ++-- .../algorithm/pv_surface_diags_alg_mod.x90 | 24 ++++++++++--------- .../app/lfric_atm/file/axis_def_plev.xml | 5 ++-- .../file/file_def_diags_oper_nwp_gl.xml | 4 ++-- .../diagnostics/diagnostics_calc_mod.F90 | 5 ++-- 5 files changed, 22 insertions(+), 20 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 18089f65d6..d282959e3e 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -992,8 +992,8 @@ - - + + plev__u_clim * plev__u_clim plev__u_clim * plev__v_clim plev__u_clim * plev__omega_clim diff --git a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 index d33edad670..4b75a2e3eb 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 @@ -54,12 +54,12 @@ contains type(field_type), intent(in) :: theta_in_w3 ! Diagnostic fields - type(field_type) :: thlev_pv, theta_at_pv2 + type(field_type) :: pv_on_theta_levs, theta_on_pv2 ! Working fields type(field_type) :: mod_pv - logical(l_def) :: thlev_pv_flag, theta_at_pv2_flag + logical(l_def) :: pv_on_theta_levs_flag, theta_on_pv2_flag integer(i_def) :: nthlev real(r_def), allocatable :: thlevs(:) @@ -67,10 +67,12 @@ contains integer(tik) :: id - thlev_pv_flag = init_diag(thlev_pv, 'thlev__pv') - theta_at_pv2_flag = init_diag(theta_at_pv2, 'dyn_trop__theta') + 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. (thlev_pv_flag .or. theta_at_pv2_flag)) return + 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' ) @@ -78,7 +80,7 @@ contains ! 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 (thlev_pv_flag) then + if (pv_on_theta_levs_flag) then nthlev = get_axis_dimension('theta_levels') allocate(thlevs(nthlev)) @@ -86,9 +88,9 @@ contains call invoke_level_interp_mdi_kernel_type(pv, theta_in_w3, & nthlev, thlevs, & - thlev_pv, & + pv_on_theta_levs, & interp_order_cubic) - call thlev_pv%write_field() + call pv_on_theta_levs%write_field() deallocate(thlevs) @@ -96,7 +98,7 @@ contains ! Theta on the PV = +/-2 surface. |PV| is far from monotonic, so linear ! interpolation to the uppermost crossing is used. - if (theta_at_pv2_flag) then + if (theta_on_pv2_flag) then call pv%copy_field_properties(mod_pv) call invoke( X_times_Y(mod_pv, pv, pv), & @@ -106,9 +108,9 @@ contains call invoke_level_interp_mdi_kernel_type(theta_in_w3, mod_pv, & 1_i_def, pv_levs, & - theta_at_pv2, & + theta_on_pv2, & interp_order_linear) - call theta_at_pv2%write_field() + call theta_on_pv2%write_field() end if diff --git a/rose-stem/app/lfric_atm/file/axis_def_plev.xml b/rose-stem/app/lfric_atm/file/axis_def_plev.xml index ea5ca39992..02beb88a72 100644 --- a/rose-stem/app/lfric_atm/file/axis_def_plev.xml +++ b/rose-stem/app/lfric_atm/file/axis_def_plev.xml @@ -1,9 +1,8 @@ - - + + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index 02d5377c86..a27c6afa9d 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -65,8 +65,8 @@ - - + + diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 index d7a13dce71..bb764c4857 100755 --- a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 +++ b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 @@ -285,8 +285,9 @@ subroutine write_pv_diagnostic(u_field, theta, rho, exner, & ! Theta surface and PV surface diagnostics, both of which need PV on ! model levels - pv_surface_flag = diagnostic_to_be_sampled('thlev__pv') .or. & - diagnostic_to_be_sampled('dyn_trop__theta') + pv_surface_flag = & + diagnostic_to_be_sampled('processed__pv_on_theta_levs') .or. & + diagnostic_to_be_sampled('processed__theta_on_pv2') if (pv_modlev_flag .or. plev_pv_flag .or. pv_surface_flag) then From cac5e896a2e3be795ba8636ba44c400c4de09b1a Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Mon, 27 Jul 2026 15:02:53 +0100 Subject: [PATCH 03/10] fix xml files --- applications/lfric_atm/metadata/field_def_diags.xml | 4 ++-- rose-stem/app/lfric_atm/file/axis_def_plev.xml | 2 +- rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index d282959e3e..fa982b4d2a 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -971,6 +971,8 @@ $DT * processed__total_snow $DT * processed__total_prec theta * (1 + 0.61 * processed__qv) + + @@ -992,8 +994,6 @@ - - plev__u_clim * plev__u_clim plev__u_clim * plev__v_clim plev__u_clim * plev__omega_clim diff --git a/rose-stem/app/lfric_atm/file/axis_def_plev.xml b/rose-stem/app/lfric_atm/file/axis_def_plev.xml index 02beb88a72..b0ea4fb001 100644 --- a/rose-stem/app/lfric_atm/file/axis_def_plev.xml +++ b/rose-stem/app/lfric_atm/file/axis_def_plev.xml @@ -2,7 +2,7 @@ - + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index a27c6afa9d..a0a6ffb166 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -65,8 +65,6 @@ - - @@ -142,7 +140,9 @@ - + + + From 77e7ee14445fd29a1e975da87139e566fcb7d60a Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Mon, 27 Jul 2026 16:30:57 +0100 Subject: [PATCH 04/10] a couple more fixes --- .../source/algorithm/pv_surface_diags_alg_mod.x90 | 3 ++- rose-stem/app/lfric_atm/file/axis_def_plev.xml | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 index 4b75a2e3eb..6f96448719 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 @@ -101,7 +101,8 @@ contains if (theta_on_pv2_flag) then call pv%copy_field_properties(mod_pv) - call invoke( X_times_Y(mod_pv, pv, 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 diff --git a/rose-stem/app/lfric_atm/file/axis_def_plev.xml b/rose-stem/app/lfric_atm/file/axis_def_plev.xml index b0ea4fb001..6047e4e7d1 100644 --- a/rose-stem/app/lfric_atm/file/axis_def_plev.xml +++ b/rose-stem/app/lfric_atm/file/axis_def_plev.xml @@ -1,10 +1,10 @@ - - + + From 08f8247210555a09509e276ba55c27d7b95de3d7 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Tue, 28 Jul 2026 09:16:14 +0000 Subject: [PATCH 05/10] Add near surface wind diagnostics at 50 m above ground Adds processed__u50m and processed__v50m, the eastward and northward wind components interpolated onto a height surface 50 m above the ground. Both reuse the existing u_in_w3 and v_in_w3 derived fields and the level_interp_mdi kernel, so only the vertical coordinate is new: height_above_surface_kernel_mod subtracts the lowest Wtheta level, this being the orography, from the W3 heights above mean sea level. Interpolation is linear. Some of the content of this change has been produced with the assistance of Anthropic Claude Opus 5 (Claude Code). --- .../lfric_atm/metadata/field_def_diags.xml | 2 + .../algorithm/height_lev_diags_alg_mod.x90 | 119 +++++++++++ .../height_above_surface_kernel_mod.F90 | 97 +++++++++ .../height_above_surface_kernel_mod_test.pf | 185 ++++++++++++++++++ .../file/file_def_diags_oper_nwp_gl.xml | 2 + .../driver/gungho_diagnostics_driver_mod.F90 | 3 + 6 files changed, 408 insertions(+) create mode 100644 interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 create mode 100644 interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 create mode 100644 interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index fa982b4d2a..c041e5cdb8 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -973,6 +973,8 @@ theta * (1 + 0.61 * processed__qv) + + diff --git a/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 new file mode 100644 index 0000000000..598effb29a --- /dev/null +++ b/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 +! Anthropic Claude Opus 5 (Claude Code). +!> @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 => null() + type(field_type), pointer :: v_in_w3 => null() + type(field_type), pointer :: height_w3 => null() + type(field_type), pointer :: height_wth => null() + + ! 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 diff --git a/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 new file mode 100644 index 0000000000..f0190d362d --- /dev/null +++ b/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 +! Anthropic Claude Opus 5 (Claude Code). +!> @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 diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf new file mode 100644 index 0000000000..a6545e3940 --- /dev/null +++ b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf @@ -0,0 +1,185 @@ +! ----------------------------------------------------------------------------- +! (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 +! Anthropic Claude Opus 5 (Claude Code). +!> @brief Test height_above_surface_kernel +!> +!> The kernel turns W3 heights above mean sea level into heights above the +!> ground, using the lowest Wtheta level as the orography. The tests check +!> that the orography is removed, that the answer over high ground matches +!> the answer at sea level, and that only the W3 part of the column is +!> written. + +module height_above_surface_kernel_mod_test + + use funit + use constants_mod, only: i_def, r_def + use height_above_surface_kernel_mod, only: height_above_surface_code + + implicit none + + private + public :: test_over_orography, test_matches_sea_level, test_second_column + + @TestCase + type, extends(TestCase), public :: height_above_surface_test_type + private + real(kind=r_def), allocatable :: height_w3(:) + real(kind=r_def), allocatable :: height_wth(:) + real(kind=r_def), allocatable :: height_agl(:) + + contains + procedure setUp + procedure tearDown + procedure test_over_orography + procedure test_matches_sea_level + procedure test_second_column + + end type height_above_surface_test_type + + ! Six layers, with room for two columns side by side in the same array. + integer(kind=i_def), parameter :: nlayers = 6_i_def + integer(kind=i_def), parameter :: ndf_w3 = 1_i_def + integer(kind=i_def), parameter :: ndf_wth = 2_i_def + integer(kind=i_def), parameter :: undf_w3 = 2_i_def * nlayers + integer(kind=i_def), parameter :: undf_wth = 2_i_def * ( nlayers + 1_i_def ) + + ! Height of the Wtheta levels above the orography, so entry one is zero + real(kind=r_def), parameter :: eta_wth(nlayers+1) = & + (/ 0.0_r_def, 20.0_r_def, 60.0_r_def, 140.0_r_def, & + 300.0_r_def, 620.0_r_def, 1260.0_r_def /) + + ! W3 levels sit midway between the Wtheta levels, so the lowest is at 10 m + real(kind=r_def), parameter :: eta_w3(nlayers) = & + (/ 10.0_r_def, 40.0_r_def, 100.0_r_def, 220.0_r_def, & + 460.0_r_def, 940.0_r_def /) + + real(kind=r_def), parameter :: tolerance = 1.0e-12_r_def + +contains + + @before + subroutine setUp( this ) + + implicit none + + class(height_above_surface_test_type), intent(inout) :: this + + allocate( this%height_w3(undf_w3) ) + allocate( this%height_wth(undf_wth) ) + allocate( this%height_agl(undf_w3) ) + + this%height_w3(:) = 0.0_r_def + this%height_wth(:) = 0.0_r_def + this%height_agl(:) = -999.0_r_def + + end subroutine setUp + + @after + subroutine tearDown( this ) + + implicit none + + class(height_above_surface_test_type), intent(inout) :: this + + deallocate( this%height_w3 ) + deallocate( this%height_wth ) + deallocate( this%height_agl ) + + end subroutine tearDown + + !> A column over 800 m of orography. The orography is removed, and the + !> lowest level comes back at 10 m rather than zero. + @test + subroutine test_over_orography( this ) + + implicit none + + class(height_above_surface_test_type), intent(inout) :: this + + real(kind=r_def), parameter :: orography = 800.0_r_def + + integer(kind=i_def) :: k + + this%height_w3(1:nlayers) = eta_w3(:) + orography + this%height_wth(1:nlayers+1) = eta_wth(:) + orography + + call height_above_surface_code( nlayers, this%height_agl, & + this%height_w3, this%height_wth, & + ndf_w3, undf_w3, (/ 1_i_def /), & + ndf_wth, undf_wth, & + (/ 1_i_def, 2_i_def /) ) + + @assertEqual( 10.0_r_def, this%height_agl(1), tolerance ) + @assertEqual( 940.0_r_def, this%height_agl(nlayers), tolerance ) + + do k = 1, nlayers + @assertEqual( eta_w3(k), this%height_agl(k), tolerance ) + end do + + ! Only the W3 levels of this column have been written + @assertEqual( -999.0_r_def, this%height_agl(nlayers + 1_i_def), tolerance ) + + end subroutine test_over_orography + + !> The same terrain following column at sea level gives the same answer. + @test + subroutine test_matches_sea_level( this ) + + implicit none + + class(height_above_surface_test_type), intent(inout) :: this + + integer(kind=i_def) :: k + + this%height_w3(1:nlayers) = eta_w3(:) + this%height_wth(1:nlayers+1) = eta_wth(:) + + call height_above_surface_code( nlayers, this%height_agl, & + this%height_w3, this%height_wth, & + ndf_w3, undf_w3, (/ 1_i_def /), & + ndf_wth, undf_wth, & + (/ 1_i_def, 2_i_def /) ) + + do k = 1, nlayers + @assertEqual( eta_w3(k), this%height_agl(k), tolerance ) + end do + + end subroutine test_matches_sea_level + + !> Each column uses its own orography, taken from its own dofmap. + @test + subroutine test_second_column( this ) + + implicit none + + class(height_above_surface_test_type), intent(inout) :: this + + real(kind=r_def), parameter :: orography = 1500.0_r_def + + integer(kind=i_def) :: k + + this%height_w3(nlayers+1:undf_w3) = eta_w3(:) + orography + this%height_wth(nlayers+2:undf_wth) = eta_wth(:) + orography + + call height_above_surface_code( nlayers, this%height_agl, & + this%height_w3, this%height_wth, & + ndf_w3, undf_w3, & + (/ nlayers + 1_i_def /), & + ndf_wth, undf_wth, & + (/ nlayers + 2_i_def, & + nlayers + 3_i_def /) ) + + do k = 1, nlayers + @assertEqual( eta_w3(k), this%height_agl(k + nlayers), tolerance ) + end do + + ! The first column has not been touched + @assertEqual( -999.0_r_def, this%height_agl(1), tolerance ) + + end subroutine test_second_column + +end module height_above_surface_kernel_mod_test diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index a0a6ffb166..e9999189b9 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -16,6 +16,8 @@ + + diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index 58b1f7224c..4e2337a301 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -57,6 +57,7 @@ module gungho_diagnostics_driver_mod use driver_modeldb_mod, only : modeldb_type #ifdef UM_PHYSICS + use height_lev_diags_alg_mod, only : height_lev_diags_alg use pres_lev_diags_alg_mod, only : pres_lev_diags_alg use pmsl_alg_mod, only : pmsl_alg use rh_diag_alg_mod, only : rh_diag_alg @@ -337,6 +338,8 @@ subroutine gungho_diagnostics_driver( modeldb, & ! Pressure level diagnostics call pres_lev_diags_alg(modeldb%config, derived_fields, theta, exner, & mr, moist_dyn) + ! Height level diagnostics + call height_lev_diags_alg(modeldb%config, derived_fields) ! Wet bulb freezing level call freeze_lev_alg(modeldb%config,theta, mr, moist_dyn, exner_in_wth) #endif From 0d8c936a696a81c3d6850ae9d7e7297eda321eef Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Tue, 28 Jul 2026 12:27:11 +0100 Subject: [PATCH 06/10] remove null pointers --- .../source/algorithm/height_lev_diags_alg_mod.x90 | 8 ++++---- .../gungho/source/diagnostics/diagnostics_calc_mod.F90 | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 index 598effb29a..7f39cf08ec 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 @@ -56,10 +56,10 @@ contains type(field_collection_type), intent(in) :: derived_fields ! Model level fields - type(field_type), pointer :: u_in_w3 => null() - type(field_type), pointer :: v_in_w3 => null() - type(field_type), pointer :: height_w3 => null() - type(field_type), pointer :: height_wth => null() + 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 diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 index bb764c4857..864d971be3 100755 --- a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 +++ b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 @@ -273,7 +273,7 @@ subroutine write_pv_diagnostic(u_field, theta, rho, exner, & type(field_type) :: pv type(field_type) :: plev_pv - type(field_type), pointer :: theta_in_w3 => null() + type(field_type), pointer :: theta_in_w3 logical(l_def) :: pv_modlev_flag, plev_pv_flag, pv_surface_flag logical(l_def), parameter :: xi3_axis = .false. logical(l_def), parameter :: add_W3_version = .false. From b85aef505fc761ebbbfda5c531bccacf38ebd0f2 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Tue, 28 Jul 2026 14:52:35 +0100 Subject: [PATCH 07/10] fix to axis def --- rose-stem/app/lfric_atm/file/axis_def_plev.xml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/rose-stem/app/lfric_atm/file/axis_def_plev.xml b/rose-stem/app/lfric_atm/file/axis_def_plev.xml index 6047e4e7d1..ac8988daa9 100644 --- a/rose-stem/app/lfric_atm/file/axis_def_plev.xml +++ b/rose-stem/app/lfric_atm/file/axis_def_plev.xml @@ -1,12 +1,16 @@ + + + + + - - + From 1fd7b570b5ab1eb9ce230e39abf75d0fd4f388e1 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Fri, 14 Aug 2026 12:45:30 +0100 Subject: [PATCH 08/10] amersand alignment --- .../height_above_surface_kernel_mod.F90 | 2 +- .../height_above_surface_kernel_mod_test.pf | 6 +- .../level_interp_mdi_kernel_mod_test.pf | 70 +++++++++---------- .../diagnostics/diagnostics_calc_mod.F90 | 2 +- 4 files changed, 40 insertions(+), 40 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 index f0190d362d..950c4c4e6b 100644 --- a/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 @@ -24,7 +24,7 @@ module height_above_surface_kernel_mod !> Kernel metadata for PSyclone type, public, extends(kernel_type) :: height_above_surface_kernel_type private - type(arg_type) :: meta_args(3) = (/ & + 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 diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf index a6545e3940..026a0cf455 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf @@ -107,7 +107,7 @@ contains this%height_w3(1:nlayers) = eta_w3(:) + orography this%height_wth(1:nlayers+1) = eta_wth(:) + orography - call height_above_surface_code( nlayers, this%height_agl, & + call height_above_surface_code( nlayers, this%height_agl, & this%height_w3, this%height_wth, & ndf_w3, undf_w3, (/ 1_i_def /), & ndf_wth, undf_wth, & @@ -138,7 +138,7 @@ contains this%height_w3(1:nlayers) = eta_w3(:) this%height_wth(1:nlayers+1) = eta_wth(:) - call height_above_surface_code( nlayers, this%height_agl, & + call height_above_surface_code( nlayers, this%height_agl, & this%height_w3, this%height_wth, & ndf_w3, undf_w3, (/ 1_i_def /), & ndf_wth, undf_wth, & @@ -165,7 +165,7 @@ contains this%height_w3(nlayers+1:undf_w3) = eta_w3(:) + orography this%height_wth(nlayers+2:undf_wth) = eta_wth(:) + orography - call height_above_surface_code( nlayers, this%height_agl, & + call height_above_surface_code( nlayers, this%height_agl, & this%height_w3, this%height_wth, & ndf_w3, undf_w3, & (/ nlayers + 1_i_def /), & diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf index 18c7291a94..62465ebec0 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf @@ -103,13 +103,13 @@ contains target_levs = (/ 295.0_r_def, 325.0_r_def /) expect = 2.0_r_def*target_levs - call level_interp_mdi_code( nlayers, & - this%data_in, & - this%coord_in, & - nlev_out, & - target_levs, & - data_out, & - interp_order_linear, & + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_linear, & ndf_in, undf_in, map_in, & ndf_out, nlev_out, map_out ) @@ -146,13 +146,13 @@ contains target_levs(1) = 3.5_r_def expect = cubic_fn( target_levs(1) ) - call level_interp_mdi_code( nlayers, & - this%data_in, & - this%coord_in, & - nlev_out, & - target_levs, & - data_out, & - interp_order_cubic, & + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_cubic, & ndf_in, undf_in, map_in, & ndf_out, nlev_out, map_out ) @@ -183,13 +183,13 @@ contains ! Below the bottom of the column, above the top of the column, and inside target_levs = (/ 270.0_r_def, 340.0_r_def, 305.0_r_def /) - call level_interp_mdi_code( nlayers, & - this%data_in, & - this%coord_in, & - nlev_out, & - target_levs, & - data_out, & - interp_order_linear, & + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_linear, & ndf_in, undf_in, map_in, & ndf_out, nlev_out, map_out ) @@ -225,13 +225,13 @@ contains target_levs(1) = 2.0_r_def - call level_interp_mdi_code( nlayers, & - this%data_in, & - this%coord_in, & - nlev_out, & - target_levs, & - data_out, & - interp_order_linear, & + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_linear, & ndf_in, undf_in, map_in, & ndf_out, nlev_out, map_out ) @@ -266,13 +266,13 @@ contains ! Bottom interval (level_below = 0) and top interval (level_below = 4) target_levs = (/ 0.5_r_def, 4.5_r_def /) - call level_interp_mdi_code( nlayers, & - this%data_in, & - this%coord_in, & - nlev_out, & - target_levs, & - data_out, & - interp_order_cubic, & + call level_interp_mdi_code( nlayers, & + this%data_in, & + this%coord_in, & + nlev_out, & + target_levs, & + data_out, & + interp_order_cubic, & ndf_in, undf_in, map_in, & ndf_out, nlev_out, map_out ) diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 index 864d971be3..f4f0655ef8 100755 --- a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 +++ b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 @@ -259,7 +259,7 @@ end subroutine write_vorticity_diagnostic subroutine write_pv_diagnostic(u_field, theta, rho, exner, & derived_fields, clock) - use pres_lev_diags_alg_mod, only: pres_lev_field_alg + use pres_lev_diags_alg_mod, only: pres_lev_field_alg use pv_surface_diags_alg_mod, only: pv_surface_diags_alg implicit none From 22c6ac18f9b9da2a77b2adc1771db112ba3ef2be Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Tue, 18 Aug 2026 14:48:20 +0100 Subject: [PATCH 09/10] science review comments --- .../kernel/level_interp_mdi_kernel_mod.F90 | 15 ++--- .../height_above_surface_kernel_mod_test.pf | 64 +------------------ 2 files changed, 8 insertions(+), 71 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 index c9b51b031c..e10105f81d 100644 --- a/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 @@ -124,15 +124,12 @@ subroutine level_interp_mdi_code(nlayers, & ! Highest level whose coordinate value is at or below the target, leaving ! room for the level above. level_below = -1 flags missing data. level_below = -1_i_def - do k = 0, top_df - 1 - if ( coord_in(base_in+k) <= target_lev ) then - level_below = k - end if - end do - - ! Target lies above the top of the column - if ( target_lev > coord_in(base_in+top_df) ) then - level_below = -1_i_def + if (target_lev <= coord_in(base_in+top_df) ) then + do k = 0, top_df - 1 + if ( coord_in(base_in+k) <= target_lev ) then + level_below = k + end if + end do end if if ( level_below == -1_i_def ) then diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf index 026a0cf455..149542ca1c 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf @@ -10,8 +10,7 @@ !> The kernel turns W3 heights above mean sea level into heights above the !> ground, using the lowest Wtheta level as the orography. The tests check !> that the orography is removed, that the answer over high ground matches -!> the answer at sea level, and that only the W3 part of the column is -!> written. +!> the answer at sea level. module height_above_surface_kernel_mod_test @@ -22,7 +21,7 @@ module height_above_surface_kernel_mod_test implicit none private - public :: test_over_orography, test_matches_sea_level, test_second_column + public :: test_over_orography @TestCase type, extends(TestCase), public :: height_above_surface_test_type @@ -35,8 +34,6 @@ module height_above_surface_kernel_mod_test procedure setUp procedure tearDown procedure test_over_orography - procedure test_matches_sea_level - procedure test_second_column end type height_above_surface_test_type @@ -125,61 +122,4 @@ contains end subroutine test_over_orography - !> The same terrain following column at sea level gives the same answer. - @test - subroutine test_matches_sea_level( this ) - - implicit none - - class(height_above_surface_test_type), intent(inout) :: this - - integer(kind=i_def) :: k - - this%height_w3(1:nlayers) = eta_w3(:) - this%height_wth(1:nlayers+1) = eta_wth(:) - - call height_above_surface_code( nlayers, this%height_agl, & - this%height_w3, this%height_wth, & - ndf_w3, undf_w3, (/ 1_i_def /), & - ndf_wth, undf_wth, & - (/ 1_i_def, 2_i_def /) ) - - do k = 1, nlayers - @assertEqual( eta_w3(k), this%height_agl(k), tolerance ) - end do - - end subroutine test_matches_sea_level - - !> Each column uses its own orography, taken from its own dofmap. - @test - subroutine test_second_column( this ) - - implicit none - - class(height_above_surface_test_type), intent(inout) :: this - - real(kind=r_def), parameter :: orography = 1500.0_r_def - - integer(kind=i_def) :: k - - this%height_w3(nlayers+1:undf_w3) = eta_w3(:) + orography - this%height_wth(nlayers+2:undf_wth) = eta_wth(:) + orography - - call height_above_surface_code( nlayers, this%height_agl, & - this%height_w3, this%height_wth, & - ndf_w3, undf_w3, & - (/ nlayers + 1_i_def /), & - ndf_wth, undf_wth, & - (/ nlayers + 2_i_def, & - nlayers + 3_i_def /) ) - - do k = 1, nlayers - @assertEqual( eta_w3(k), this%height_agl(k + nlayers), tolerance ) - end do - - ! The first column has not been touched - @assertEqual( -999.0_r_def, this%height_agl(1), tolerance ) - - end subroutine test_second_column - end module height_above_surface_kernel_mod_test From a88c7f2b829f67c0f37ed7e57756acdc029d6397 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Tue, 15 Sep 2026 10:18:07 +0100 Subject: [PATCH 10/10] change claude copyright statement --- .../source/algorithm/height_lev_diags_alg_mod.x90 | 2 +- .../source/algorithm/pv_surface_diags_alg_mod.x90 | 2 +- .../source/kernel/height_above_surface_kernel_mod.F90 | 2 +- .../source/kernel/level_interp_mdi_kernel_mod.F90 | 2 +- .../source/psy/psykal_lite_phys_mod.F90 | 2 +- .../unit-test/kernel/height_above_surface_kernel_mod_test.pf | 2 +- .../unit-test/kernel/level_interp_mdi_kernel_mod_test.pf | 2 +- science/gungho/source/diagnostics/diagnostics_calc_mod.F90 | 2 +- science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 | 2 +- 9 files changed, 9 insertions(+), 9 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 index 7f39cf08ec..a90274ba49 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/height_lev_diags_alg_mod.x90 @@ -4,7 +4,7 @@ ! under which the code may be used. !------------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! 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. diff --git a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 index 6f96448719..c2801f3d9f 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pv_surface_diags_alg_mod.x90 @@ -4,7 +4,7 @@ ! under which the code may be used. !------------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! 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 diff --git a/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 index 950c4c4e6b..226547a0ab 100644 --- a/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/height_above_surface_kernel_mod.F90 @@ -4,7 +4,7 @@ ! under which the code may be used. !------------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Height of the W3 levels above the surface. module height_above_surface_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 index e10105f81d..61514aa17e 100644 --- a/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 @@ -4,7 +4,7 @@ ! under which the code may be used. !------------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Interpolate a field onto surfaces of constant value of another field. module level_interp_mdi_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 index 240c6226cb..effe4b0547 100644 --- a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 +++ b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 @@ -4,7 +4,7 @@ ! under which the code may be used. !---------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Provides an implementation of the Psy layer for physics !> @details Contains hand-rolled versions of the Psy layer that can be used for diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf index 149542ca1c..b46d8e4af3 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf @@ -4,7 +4,7 @@ ! under which the code may be used. ! ----------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Test height_above_surface_kernel !> !> The kernel turns W3 heights above mean sea level into heights above the diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf index 62465ebec0..7d0adc5968 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/level_interp_mdi_kernel_mod_test.pf @@ -4,7 +4,7 @@ ! under which the code may be used. ! ----------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Test level_interp_mdi_kernel !> !> The kernel is a port of the UM's vert_interp_mdi2, so these tests pin down diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 index f4f0655ef8..dbe4a46fb8 100755 --- a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 +++ b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 @@ -5,7 +5,7 @@ !----------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Module for computing and outputting derived diagnostics !! !! @details Computes various derived diagnostics that are written out diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index 4e2337a301..9cec4417de 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -5,7 +5,7 @@ !----------------------------------------------------------------------------- ! Some of the content of this file has been produced with the assistance of -! Anthropic Claude Opus 5 (Claude Code). +! Met Office Claude Code Enterprise. !> @brief Outputs diagnostics from gungho/lfric_atm !> @details Calls the routine that generates diagnostic output for