diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 9b01795782..c041e5cdb8 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -971,6 +971,10 @@ $DT * processed__total_snow $DT * processed__total_prec 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..a90274ba49 --- /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 +! 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 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..c2801f3d9f --- /dev/null +++ b/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 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..226547a0ab --- /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 +! 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 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..61514aa17e --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/level_interp_mdi_kernel_mod.F90 @@ -0,0 +1,198 @@ +!------------------------------------------------------------------------------- +! (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 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 + 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 + + ! 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..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 @@ -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 +! 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 @@ -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/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..b46d8e4af3 --- /dev/null +++ b/interfaces/physics_schemes_interface/unit-test/kernel/height_above_surface_kernel_mod_test.pf @@ -0,0 +1,125 @@ +! ----------------------------------------------------------------------------- +! (C) Crown copyright Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +! ----------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Met Office Claude Code Enterprise. +!> @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. + +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 + + @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 + + 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 + +end module height_above_surface_kernel_mod_test 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..7d0adc5968 --- /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 +! 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 +!> 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..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,10 +1,16 @@ - + + + + + + + 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..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 @@ + + @@ -140,7 +142,9 @@ - + + + diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 index 9e8704c78b..dbe4a46fb8 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 +! Met Office Claude Code Enterprise. !> @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 + 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,24 @@ 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('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 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..9cec4417de 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 +! Met Office Claude Code Enterprise. !> @brief Outputs diagnostics from gungho/lfric_atm !> @details Calls the routine that generates diagnostic output for @@ -55,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 @@ -243,7 +246,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 @@ -334,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