From 1390b6900e35b7ffbe460e3a9a4b75f840475e6f Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Tue, 5 May 2026 09:34:23 +0100 Subject: [PATCH 01/27] Rose stem test config file --- rose-stem/app/lfric_atm/opt/rose-app-regrav.conf | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 rose-stem/app/lfric_atm/opt/rose-app-regrav.conf diff --git a/rose-stem/app/lfric_atm/opt/rose-app-regrav.conf b/rose-stem/app/lfric_atm/opt/rose-app-regrav.conf new file mode 100644 index 000000000..fa2c9058e --- /dev/null +++ b/rose-stem/app/lfric_atm/opt/rose-app-regrav.conf @@ -0,0 +1,5 @@ +[namelist:formulation] +shallow=.false. + +[namelist:initialization] +regravitate='geopot' From 95b08efa7e589e33490fbb8738cc0de98ea4dafd Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Wed, 6 May 2026 14:12:23 +0100 Subject: [PATCH 02/27] Metadata and upgrade macro --- .../lfric-gungho/HEAD/rose-meta.conf | 22 +++++++++++++++++++ .../gungho/rose-meta/lfric-gungho/versions.py | 15 +++++-------- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf index 4380951b7..7e5d04fed 100644 --- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf +++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf @@ -1971,6 +1971,7 @@ help=When true, the gravitational field is constant, as would be the under the s =from planet centre. !kind=default sort-key=Panel-A04 +trigger=namelist:initialization=regravitate: this == .false.; type=logical [namelist:formulation=si_momentum_equation] @@ -3220,6 +3221,27 @@ help=Horizontal winds are read in on a W2h space from a start dump. If sort-key=02b type=logical +[namelist:initialization=regravitate] +compulsory=true +description=Method for switching gravity +!enumeration=true +help=Method for switching gravity from constant to height-varying. + =Used when the shallow switch is false and reading a start dump + =that has been created by a constant gravity model, such as the UM. + =All methods calculate an Exner field in hydrostatic balance with + =height-varying gravity. This is done by upward integr + =Isothermodynamic: Preserve absolute temperature and pressure on + = model levels by allowing the heights of levels + = to change. Absolute temperature is then interpolated + = back onto the original model levels. + =Geopotential remap: + = Physical height and geopotential height are the same for data in the + = start dump, due to gravity being constant. +!kind=default +sort-key=Panel-A04b +value-titles=Geopotential remap, Isothermal vn1, Isothermal vn2, Isothermal vn3 +values='geopot', 'isotherm1', 'isotherm2', 'isotherm3' + [namelist:initialization=sea_ice_source] compulsory=true description=Where to read the sea ice fields from diff --git a/science/gungho/rose-meta/lfric-gungho/versions.py b/science/gungho/rose-meta/lfric-gungho/versions.py index 01798ad2b..e10b1524d 100644 --- a/science/gungho/rose-meta/lfric-gungho/versions.py +++ b/science/gungho/rose-meta/lfric-gungho/versions.py @@ -18,16 +18,13 @@ def __repr__(self): __str__ = __repr__ -""" -Copy this template and complete to add your macro +class vn31_t218(MacroUpgrade): + """Upgrade macro for issue #218 by Chris Smith.""" -class vnXX_txxx(MacroUpgrade): - # Upgrade macro for by - - BEFORE_TAG = "vnX.X" - AFTER_TAG = "vnX.X_txxx" + BEFORE_TAG = "vn3.1_t135" + AFTER_TAG = "vn3.1_t218" def upgrade(self, config, meta_config=None): - # Add settings + self.add_setting(config, ["namelist:initialization", "regravitate"], "'isotherm1'") return config, self.reports -""" + From b2c25b334e58091e7b04a5b12475fedc80ddfe9d Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Wed, 6 May 2026 14:20:06 +0100 Subject: [PATCH 03/27] Regravitation code --- .../map_fd_to_prognostics_alg_mod-merge.x90 | 608 ++++++++++++++++++ .../hstat_cori_g_const_kernel_mod.F90 | 210 ++++++ .../regrav_geopot_kernel_mod.F90 | 226 +++++++ .../regrav_isotherm_kernel_mod.F90 | 153 +++++ 4 files changed, 1197 insertions(+) create mode 100644 science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 create mode 100644 science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 create mode 100644 science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 create mode 100644 science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 new file mode 100644 index 000000000..2e60f5a3b --- /dev/null +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 @@ -0,0 +1,608 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2019 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!> @brief Set prognostic fields from FD start dump +module map_fd_to_prognostics_alg_mod + + use constants_mod, only: r_def, i_def, l_def + use log_mod, only: log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_DEBUG, & + LOG_LEVEL_TRACE + use namelist_collection_mod, only: namelist_collection_type + use namelist_mod, only: namelist_type + + ! Derived Types + use field_mod, only: field_type + use integer_field_mod, only: integer_field_type + use field_collection_mod, only: field_collection_type + use mesh_mod, only: mesh_type + use operator_mod, only: operator_type + use function_space_mod, only: function_space_type + use function_space_collection_mod, only: function_space_collection + + use quadrature_xyoz_mod, only: quadrature_xyoz_type + use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type + use sci_fem_constants_mod, only: get_rmultiplicity_fv + use sci_geometric_constants_mod, only: get_height_fv, & + get_coordinates, & + get_panel_id, & + get_da_at_w2, & + get_face_selector_ew, & + get_face_selector_ns + use sci_mapping_constants_mod, only: get_u_lon_map, & + get_u_lat_map, & + get_u_up_map, & + get_u_lon_sample, & + get_u_lat_sample, & + get_u_up_sample, & + get_w3_to_w2_displacement + use dycore_constants_mod, only: get_vert_coriolis, & + get_geopotential + use fs_continuity_mod, only: W3, W2, W2broken, Wtheta + use sci_sort_ref_kernel_mod, only: sort_ref_kernel_type + use combine_w2_field_kernel_mod, only: combine_w2_field_kernel_type + use idealised_config_mod, only: perturb_init, & + perturb_magnitude + use random_perturb_field_alg_mod, only: random_perturb_field + use mr_indices_mod, only: imr_v, imr_cl, imr_ci, imr_r, & + imr_s, imr_g + use moist_dyn_mod, only: num_moist_factors, gas_law + use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg + + ! Configuration modules + use base_mesh_config_mod, only: geometry, & + geometry_spherical, & + topology, & + topology_fully_periodic + use formulation_config_mod, only: rotating, shallow + use finite_element_config_mod, only: coord_system, coord_system_native + use initialization_config_mod, only: model_eos_height, & + read_w2h_wind, & + regravitate, & + regravitate_geopot, & + regravitate_isotherm1, & + regravitate_isotherm2, & + regravitate_isotherm3 + use physics_config_mod, only: sample_physics_winds, & + sample_physics_winds_correction + use planet_config_mod, only: gravity, p_zero, kappa, rd, cp + + implicit none + + private + public :: map_fd_to_prognostics, hydrostatic_balance, set_wind + +contains + + !> @details Setting FE prognostic fields from FD fields + !> @param[in] configuration The application configuration object + !> @param[in,out] prognostic_fields Collection of prognostic fields + !> @param[in,out] mr Array of moisture mixing ratios + !> @param[in,out] fd_fields Collection of input Finite Difference + !> fields + subroutine map_fd_to_prognostics(configuration, prognostic_fields, mr, & + moist_dyn, fd_fields) + + implicit none + + ! FE Prognostic fields + type(namelist_collection_type), intent(in) :: configuration + type(field_collection_type), intent( inout ) :: prognostic_fields + type(field_type), intent( inout ) :: mr(:) + type(field_type), intent( inout ) :: moist_dyn(num_moist_factors) + + type( field_collection_type), intent( inout ) :: fd_fields + + ! Local dereferenced fields + type( field_type ), pointer :: theta + type( field_type ), pointer :: rho + type( field_type ), pointer :: u + type( field_type ), pointer :: exner + + ! FD fields + type( field_type ), pointer :: h_wind_in_w2h + type( field_type ), pointer :: ew_wind_in_w3 + type( field_type ), pointer :: ns_wind_in_w3 + type( field_type ), pointer :: dry_rho_in_w3 + type( field_type ), pointer :: upward_wind_in_wtheta + type( field_type ), pointer :: theta_in_wtheta + type( field_type ), pointer :: mv_in_wtheta + type( field_type ), pointer :: mcl_in_wtheta + type( field_type ), pointer :: mcf_in_wtheta + type( field_type ), pointer :: mr_in_wtheta + type( field_type ), pointer :: dA + integer(i_def) :: mesh_id + + type(integer_field_type), pointer :: face_selector_ew + type(integer_field_type), pointer :: face_selector_ns + + ! FD (source) fields + if ( read_w2h_wind ) then + call fd_fields%get_field('h_wind', h_wind_in_w2h) + else + call fd_fields%get_field('ew_wind_in_w3', ew_wind_in_w3) + call fd_fields%get_field('ns_wind_in_w3', ns_wind_in_w3) + end if + + call fd_fields%get_field('dry_rho_in_w3', dry_rho_in_w3) + call fd_fields%get_field('upward_wind_in_wtheta', upward_wind_in_wtheta) + call fd_fields%get_field('theta_in_wtheta', theta_in_wtheta) + call fd_fields%get_field('mv_in_wtheta', mv_in_wtheta) + call fd_fields%get_field('mcl_in_wtheta', mcl_in_wtheta) + call fd_fields%get_field('mcf_in_wtheta', mcf_in_wtheta) + call fd_fields%get_field('mr_in_wtheta', mr_in_wtheta) + + ! Prognostic (target) fields + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('theta', theta) + call prognostic_fields%get_field('exner', exner) + call prognostic_fields%get_field('rho', rho) + + ! Some variables should coming in from the restart file, but currently + ! they are not output by UM2LFRIC and are set to zero or a default value + ! This will need to be reviewed if this becomes a checkpoint-restart + ! routine and/or um2lfric is updated for GAL physics + call invoke( setval_x(theta, theta_in_wtheta), & + setval_x(rho, dry_rho_in_w3), & + setval_x(mr(imr_v), mv_in_wtheta), & + setval_x(mr(imr_cl), mcl_in_wtheta), & + setval_x(mr(imr_s), mcf_in_wtheta), & + setval_x(mr(imr_r), mr_in_wtheta), & + setval_c(mr(imr_ci), 0.0_r_def), & + setval_c(mr(imr_g), 0.0_r_def), & + setval_c(exner, 0.0_r_def), & + setval_c(upward_wind_in_wtheta, 0.0_r_def)) + + if ( read_w2h_wind )then + ! We read in the horizontal (normal to face) winds on the w2h points + mesh_id = u%get_mesh_id() + dA => get_da_at_w2(mesh_id) + face_selector_ew => get_face_selector_ew(mesh_id) + face_selector_ns => get_face_selector_ns(mesh_id) + + call invoke( combine_w2_field_kernel_type(u, h_wind_in_w2h, & + upward_wind_in_wtheta, & + face_selector_ew, & + face_selector_ns), & + inc_X_times_Y(u, dA) ) + else + ! We read in cell centred zonal/meridional winds + call set_wind( u, ew_wind_in_w3, ns_wind_in_w3, upward_wind_in_wtheta ) + end if + + ! Update factors for moist dynamics + call moist_dyn_factors_alg( moist_dyn, mr ) + ! Remove any static instability + call invoke( sort_ref_kernel_type(theta) ) + ! Initialize hydrostatically balanced exner field + call hydrostatic_balance( exner, rho, theta, u, moist_dyn, & + model_eos_height ) + + if ( perturb_init ) then + call log_event( "Gungho: Applying initial perturbation to theta", & + LOG_LEVEL_INFO ) + call random_perturb_field( theta, perturb_magnitude ) + end if + + ! Free up any prognostics not required + if ( read_w2h_wind ) then + call fd_fields%remove_field("h_wind") + else + call fd_fields%remove_field("ew_wind_in_w3") + call fd_fields%remove_field("ns_wind_in_w3") + end if + call fd_fields%remove_field("dry_rho_in_w3") + call fd_fields%remove_field("upward_wind_in_wtheta") + call fd_fields%remove_field("theta_in_wtheta") + call fd_fields%remove_field("mv_in_wtheta") + call fd_fields%remove_field("mcl_in_wtheta") + call fd_fields%remove_field("mcf_in_wtheta") + call fd_fields%remove_field("mr_in_wtheta") + + call log_event( "Gungho: Set prognostic fields from FD fields", LOG_LEVEL_INFO ) + + end subroutine map_fd_to_prognostics + + !> @param[in,out] u Wind (FE) + !> @param[in] u_lon Zonal Wind (FD) + !> @param[in] u_lat Meridional Wind (FD) + !> @param[in] u_up Upward Wind (FD) + subroutine set_wind(u, u_lon, u_lat, u_up) + + use sci_convert_phys_to_hdiv_kernel_mod, & + only: convert_phys_to_hdiv_kernel_type + use dg_matrix_vector_kernel_mod, only: dg_matrix_vector_kernel_type + use dg_inc_matrix_vector_kernel_mod, only: dg_inc_matrix_vector_kernel_type + use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type + use sci_mass_matrix_solver_alg_mod, only: mass_matrix_solver_alg + use matrix_vector_kernel_mod, only: matrix_vector_kernel_type + use sci_average_w2b_to_w2_kernel_mod, & + only: average_w2b_to_w2_kernel_type + use sci_w3_to_w2_correction_kernel_mod, & + only: w3_to_w2_correction_kernel_type + use timing_mod, only: start_timing, stop_timing, & + tik, LPROF + + implicit none + + type( field_type ), intent( inout ) :: u + type( field_type ), intent( in ) :: u_lat, u_lon, u_up + + type( field_type ) :: r_u + type( field_type ) :: u_broken + type( field_type ) :: u_lat_w2, u_lon_w2, u_up_w2 + type( field_type ) :: u_lat_large, u_lon_large + type( field_type ) :: u_correction + type( field_type ), pointer :: rmultiplicity_w2 + type( field_type ), pointer :: displacement + type( field_type ), pointer :: chi(:) + type( field_type ), pointer :: panel_id + type( operator_type ), pointer :: u_lon_map + type( operator_type ), pointer :: u_lat_map + type( operator_type ), pointer :: u_up_map + type( function_space_type ), pointer :: w2_fs + type( function_space_type ), pointer :: w2b_fs + + type( mesh_type ), pointer :: mesh + integer(tik) :: id + + if ( LPROF ) call start_timing( id, 'mappings.set_wind' ) + + mesh => u%get_mesh() + + if (sample_physics_winds) then + + w2b_fs => function_space_collection%get_fs(mesh, 0, 0, W2broken) + rmultiplicity_w2 => get_rmultiplicity_fv(W2, mesh%get_id()) + u_lon_map => get_u_lon_sample(mesh%get_id()) + u_lat_map => get_u_lat_sample(mesh%get_id()) + u_up_map => get_u_up_sample(mesh%get_id()) + + call u_broken%initialise( w2b_fs ) + + call invoke( name="sample_physics_winds_in_w2", & + setval_c(u, 0.0_r_def), & + dg_matrix_vector_kernel_type(u_broken, u_lon, u_lon_map), & + dg_inc_matrix_vector_kernel_type(u_broken, u_lat, u_lat_map), & + dg_inc_matrix_vector_kernel_type(u_broken, u_up, u_up_map), & + average_w2b_to_w2_kernel_type(u, u_broken, rmultiplicity_w2) ) + + if (sample_physics_winds_correction & + .and. geometry == geometry_spherical & + .and. topology == topology_fully_periodic) then + + w2_fs => u%get_function_space() + call u_lon_w2%initialise( w2_fs ) + call u_lat_w2%initialise( w2_fs ) + call u_up_w2%initialise( w2_fs ) + call u_correction%initialise( w2_fs ) + call u_lon_large%initialise( u_lon%get_function_space(), halo_depth=2 ) + call u_lat_large%initialise( u_lat%get_function_space(), halo_depth=2 ) + + displacement => get_w3_to_w2_displacement(mesh%get_id()) + chi => get_coordinates(mesh%get_id()) + panel_id => get_panel_id(mesh%get_id()) + + call invoke( setval_c(u_lon_w2, 0.0_r_def), & + setval_c(u_lat_w2, 0.0_r_def), & + setval_c(u_up_w2, 0.0_r_def), & + ! Copy u_lon and u_lat to fields with appropriate depth + setval_X(u_lon_large, u_lon), & + setval_X(u_lat_large, u_lat), & + ! Correct horizontal wind components (but not vertical) + w3_to_w2_correction_kernel_type(u_lon_w2, u_lon_large, 1, & + displacement, & + panel_id, 1), & + w3_to_w2_correction_kernel_type(u_lat_w2, u_lat_large, 1, & + displacement, & + panel_id, 1), & + ! Convert components to HDiv field + convert_phys_to_hdiv_kernel_type(u_correction, & + u_lon_w2, u_lat_w2, & + u_up_w2, chi, panel_id, & + geometry), & + ! Increment wind with correction field + inc_X_plus_Y(u, u_correction) ) + end if + + else + + call u%copy_field_properties(r_u) + + u_lon_map => get_u_lon_map( mesh%get_id() ) + u_lat_map => get_u_lat_map( mesh%get_id() ) + u_up_map => get_u_up_map( mesh%get_id() ) + + call invoke( name="galerkin_projection_from_lonlatr", & + setval_c( u, 0.0_r_def ), & + setval_c( r_u, 0.0_r_def ), & + matrix_vector_kernel_type(r_u, u_lon, u_lon_map), & + matrix_vector_kernel_type(r_u, u_lat, u_lat_map), & + matrix_vector_kernel_type(r_u, u_up, u_up_map), & + enforce_bc_kernel_type( r_u ) ) + call mass_matrix_solver_alg(u, r_u) + end if + + if ( LPROF ) call stop_timing( id, 'mappings.set_wind' ) + + end subroutine set_wind + + !> @param[in,out] exner Exner pressure + !> @param[in] rho Dry density + !> @param[in] theta Theta + !> @param[in] u Wind + !> @param[in] moist_dyn Moisture arrays for dynamics + !> @param[in] lbc_mask Optional LBC mask + subroutine hydrostatic_balance(exner, rho, theta, u, moist_dyn, & + eos_height, lbc_mask) + use dg_inc_matrix_vector_kernel_mod, & + only : dg_inc_matrix_vector_kernel_type + use sci_enforce_bc_kernel_mod, only : enforce_bc_kernel_type + use hydrostatic_coriolis_kernel_mod, & + only : hydrostatic_coriolis_kernel_type + use hstat_cori_g_const_kernel_mod, & + only : hstat_cori_g_const_kernel_type + use hydro_shallow_to_deep_kernel_mod, & + only : hydro_shallow_to_deep_kernel_type + use regrav_isotherm_kernel_mod, only : regrav_isotherm_kernel_type + use regrav_geopot_kernel_mod, only : regrav_geopot_kernel_type + use moist_dyn_mod, only : num_moist_factors + use matrix_vector_kernel_mod, only : matrix_vector_kernel_type + use sci_nodal_xyz_coordinates_kernel_mod, & + only : nodal_xyz_coordinates_kernel_type + use sci_mass_matrix_solver_alg_mod, & + only : mass_matrix_solver_alg + use sample_eos_pressure_kernel_mod, & + only : sample_eos_pressure_kernel_type + use sample_eos_rho_kernel_mod, only : sample_eos_rho_kernel_type + use sci_sample_w3_to_wtheta_kernel_mod, & + only : sample_w3_to_wtheta_kernel_type + use sci_convert_cart2sphere_vector_kernel_mod, & + only: convert_cart2sphere_vector_kernel_type + + implicit none + + type( field_type ), intent( inout ) :: exner + type( field_type ), intent( inout ) :: theta, rho, u + type( field_type ), intent( in ) :: moist_dyn(num_moist_factors) + integer(i_def), intent( in ) :: eos_height + type( field_type ), intent( in ), optional, target :: lbc_mask + + type(field_type ) :: c_vert_term, r_c + type(field_type ) :: exner_in_wth + type(field_type ) :: temperature + type(field_type), pointer :: geopotential + type(field_type), pointer :: chi(:) + type(field_type), pointer :: panel_id + type(field_type), pointer :: height_w3 + type(field_type), pointer :: height_wth + type(operator_type), pointer :: vert_coriolis + type(mesh_type), pointer :: mesh + type(field_type), pointer :: mask + type(field_type), target :: dummy_mask + integer(i_def) :: eos_level + + nullify(geopotential, chi, panel_id, height_w3, & + height_wth, vert_coriolis, mesh, mask ) + + mesh => exner%get_mesh() + height_w3 => get_height_fv( W3, mesh%get_id() ) + height_wth => get_height_fv( Wtheta, mesh%get_id() ) + geopotential => get_geopotential( mesh%get_id() ) + + ! If lbc_mask wasn't supplied, create a dummy mask + if (present(lbc_mask)) then + mask => lbc_mask + else + call dummy_mask%initialise( vector_space = exner%get_function_space() ) + call invoke( setval_c(dummy_mask, 1.0_r_def) ) + mask => dummy_mask + end if + + chi => get_coordinates( mesh%get_id() ) + panel_id => get_panel_id( mesh%get_id() ) + vert_coriolis => get_vert_coriolis( mesh%get_id() ) + + ! Create intermediate fields. + call theta%copy_field_properties(c_vert_term) + call invoke( setval_c(c_vert_term, 0.0_r_def) ) + + ! Compute contribution of Coriolis force to vertical balance + if ( associated(vert_coriolis) ) then + + call r_c%initialise( theta%get_function_space() ) + + ! Apply Coriolis mass matrix, and solve for vertical term in Wtheta + call invoke( setval_c(r_c, 0.0_r_def), & + dg_inc_matrix_vector_kernel_type(r_c, u, vert_coriolis), & + enforce_bc_kernel_type(r_c) ) + call mass_matrix_solver_alg(c_vert_term, r_c, bc_flag=.true.) + ! Enforce that result is zero at domain bottom and top + call invoke( enforce_bc_kernel_type(c_vert_term) ) + end if ! associated(vert_coriolis) + + ! Get the index of the level closest to the height, ignoring the top eta + ! level as eos_level is to be used for a W3 space rather than Wtheta + eos_level = get_nearest_level( mesh, eos_height ) + if ( eos_level == (mesh%get_nlayers() + 1) ) eos_level = eos_level - 1 + + if (shallow) then + ! Initialize the vertical profile by starting with the equation of state at + ! eos_level, and then integrating from this level to the top and the bottom + ! using hydrostatic balance with Coriolis terms. + call invoke( hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ) ) + else + ! Initialize the vertical profile of exner pressure by initializing + ! with the equation of state, and then integrating from eos_level to + ! the surface using hydrostatic balance with Coriolis terms and + ! absolute temperature and constant gravity, followed by integrating + ! from the surface to the top with the model geopotential, followed + ! by rederiving potential temperature (theta) and density (rho). + call theta%copy_field_properties(exner_in_wth) + call theta%copy_field_properties(temperature) + + select case( regravitate ) + case( regravitate_isotherm1 ) + call invoke( & + ! Compute exner from EoS + sample_eos_pressure_kernel_type( exner, rho, & + theta, moist_dyn(gas_law) ), & + + ! Compute exner in Wtheta (required to calculate T) + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + + ! Rebalance exner in vertical, using hydrostatic balance. + ! This integrates down using g (shallow) from eos_level and + ! then integrates up using phi (deep). + hydro_shallow_to_deep_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn, geopotential, & + height_w3, height_wth, & + mask, gravity, & + cp, eos_level ) ) + case( regravitate_isotherm2 ) + call invoke( & + ! Compute exner from EoS + sample_eos_pressure_kernel_type( exner, rho, & + theta, moist_dyn(gas_law) ), & + + ! Compute exner in Wtheta (required to calculate T) + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + ! Rebalance exner in vertical, using hydrostatic balance. + ! This integrates down using g (shallow) from eos_level and + ! then integrates up using phi (deep). + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask ) ) + case( regravitate_isotherm3 ) + call invoke( & + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask ) ) + case( regravitate_geopot ) + call invoke( & + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + regrav_geopot_kernel_type( temperature, theta, & + exner_in_wth, exner, & + c_vert_term, & + moist_dyn, & + height_w3, height_wth, & + mask ), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn, geopotential, & + height_w3, height_wth, & + mask ) ) + end select + call invoke( & + ! Recompute exner in Wtheta (required to calculate theta) + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + + ! Recompute potential temperature, theta = T /exner + X_divideby_Y( theta, temperature, exner_in_wth ), & + + ! Recompute density using the equation of state + sample_eos_rho_kernel_type(rho, exner, theta, & + moist_dyn(gas_law), & + kappa, rd, p_zero) & + ) + + end if + + nullify( mesh ) + + end subroutine hydrostatic_balance + + !> @brief Returns the index of the vertical level closest to the input height + !> @details Find the index of the vertical layer that is closest to a given + !! height that is specified as an integer percentage of the total height. + !> @param[in] mesh The mesh from which to obtain the eta vertical levels + !> @param[in] height The height at which to apply the equation of state, + !> specified as a percentage of the total height of the + !> domain + !> @return level Index of vertical level nearest to specified [%] height. + function get_nearest_level(mesh, height) result(level) + + type( mesh_type ), intent( in ) :: mesh + integer(i_def), intent( in ) :: height + integer(i_def) :: level + real(r_def), allocatable :: levels(:), abs_levels(:) + integer(i_def) :: eta_level_number + + eta_level_number = mesh%get_nlayers() + 1 + allocate( levels( eta_level_number ) ) + allocate( abs_levels( eta_level_number ) ) + + call mesh%get_eta( levels ) + + ! Find the index of the eta height level that is closest to the height + ! (given as a percentage) + abs_levels = abs( levels - ( 0.01_r_def * height ) ) + level = minloc( abs_levels, dim = 1 ) + + deallocate( levels, abs_levels ) + + return + + end function get_nearest_level + +end module map_fd_to_prognostics_alg_mod diff --git a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 new file mode 100644 index 000000000..f7809f661 --- /dev/null +++ b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 @@ -0,0 +1,210 @@ +!----------------------------------------------------------------------------- +! (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. +!----------------------------------------------------------------------------- + +module hstat_cori_g_const_kernel_mod + +!> @brief Computes exner from equation of state and vertical balance including +!! Coriolis term. +!> @details Calculate exner on the top level, using equation of state. Then +!! use finite differencing to calculate exner on successive levels +!! below. Unlike hydrostatic_eos_kernel which balances upwards, the +!! Coriolis term is also added to the vertical balance equation. + +use argument_mod, only : arg_type, func_type, & + GH_FIELD, GH_REAL, & + GH_SCALAR, GH_INTEGER, & + GH_READ, GH_WRITE, & + ANY_SPACE_9, ANY_SPACE_1, & + GH_BASIS, CELL_COLUMN, GH_EVALUATOR +use constants_mod, only : r_def, i_def +use fs_continuity_mod, only : Wtheta, W3 +use kernel_mod, only : kernel_type +use planet_config_mod, only : gravity + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +type, public, extends(kernel_type) :: hstat_cori_g_const_kernel_type + private + type(arg_type) :: meta_args(13) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & + /) + type(func_type) :: meta_funcs(2) = (/ & + func_type(W3, GH_BASIS), & + func_type(Wtheta, GH_BASIS) & + /) + integer :: operates_on = CELL_COLUMN + integer :: gh_shape = GH_EVALUATOR +contains + procedure, nopass :: hstat_cori_g_const_code +end type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- +public :: hstat_cori_g_const_code + +contains + +!! @param[in] nlayers Number of layers +!! @param[in,out] exner Exner pressure field +!! @param[in] rho Density field +!! @param[in] theta Potential temperature field +!! @param[in] coriolis_term Vertical component of the coriolis term +!! @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!! @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!! @param[in] height_w3 Height coordinate in w3 +!! @param[in] w3_mask LBC mask or Dummy mask for w3 space +!! @param[in] p_zero Reference surface pressure +!! @param[in] kappa Ratio of Rd and cp +!! @param[in] rd Gas constant for dry air +!! @param[in] cp Specific heat of dry air at constant pressure +!! @param[in] eos_index Vertical level at which the equation of state +!! is satisfied +!! @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!! @param[in] undf_w3 Total number of degrees of freedom for w3 +!! @param[in] map_w3 Dofmap for the cell at column base for w3 +!! @param[in] basis_w3 Basis functions evaluated at w3 nodes +!! @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!! @param[in] undf_wt Total number of degrees of freedom for wtheta +!! @param[in] map_wt Dofmap for the cell at column base for wt +!! @param[in] basis_wt Basis functions evaluated at wt nodes +subroutine hstat_cori_g_const_code( nlayers, & + exner, & + rho, & + theta, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + height_w3, & + w3_mask, & + p_zero, & + kappa, & + rd, & + cp, & + eos_index, & + ndf_w3, & + undf_w3, & + map_w3, & + basis_w3, & + ndf_wt, & + undf_wt, & + map_wt, & + basis_wt ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + integer(kind=i_def), intent(in) :: eos_index + + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + real(kind=r_def), dimension(undf_w3), intent(in) :: rho, & + height_w3, & + w3_mask + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot + real(kind=r_def), dimension(undf_wt), intent(in) :: theta + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), dimension(1,ndf_w3,ndf_w3), intent(in) :: basis_w3 + real(kind=r_def), dimension(1,ndf_wt,ndf_w3), intent(in) :: basis_wt + real(kind=r_def), intent(in) :: p_zero + real(kind=r_def), intent(in) :: kappa + real(kind=r_def), intent(in) :: rd + real(kind=r_def), intent(in) :: cp + + ! Internal variables + integer(kind=i_def) :: k, df, dft, df3, eos_index_m1 + real(kind=r_def), dimension(ndf_w3) :: rho_e + real(kind=r_def), dimension(ndf_wt) :: theta_moist_e + real(kind=r_def) :: rho_cell, theta_moist, dz + + ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) + if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then + return + end if + + ! Compute exner from eqn of state at vertical-level eos_index + ! Levels in the kernel are numbered 0 to nlayers-1, rather than 1 to nlayers + + eos_index_m1 = eos_index - 1 + + k = eos_index_m1 + + do df3 = 1, ndf_w3 + rho_e( df3 ) = rho( map_w3(df3) + k) + end do + + do dft = 1, ndf_wt + theta_moist_e( dft ) = theta( map_wt(dft) + k) * moist_dyn_gas( map_wt(dft) + k ) + end do + + do df = 1, ndf_w3 + + rho_cell = 0.0_r_def + do df3 = 1, ndf_w3 + rho_cell = rho_cell + rho_e( df3 )*basis_w3( 1,df3,df ) + end do + + theta_moist = 0.0_r_def + do dft = 1, ndf_wt + theta_moist = theta_moist + theta_moist_e( dft )*basis_wt( 1,dft,df ) + end do + + exner( map_w3(df) + k ) = ( rd*rho_cell*theta_moist/p_zero ) & + **( kappa/( 1.0_r_def-kappa ) ) + end do + + ! Exner on other levels from hydrostatic balance + ! Compute exner at level k-1 from exner at level k plus other terms. + ! Add the coriolis term using the bottom face (B) from the cell at level k-1 + do k = eos_index_m1, 1, -1 + + dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) + theta_moist = moist_dyn_gas( map_wt(1) + k ) * theta( map_wt(1) + k ) / & + moist_dyn_tot( map_wt(1) + k ) + exner( map_w3(1) + k - 1 ) = exner( map_w3 (1) + k ) & + + ( gravity - coriolis_term( map_wt(1) + k ) ) * dz & + / ( cp * theta_moist ) + + end do + + do k = eos_index_m1, nlayers-2 + + dz = height_w3( map_w3(1) + k + 1 ) - height_w3( map_w3(1) + k ) + theta_moist = moist_dyn_gas( map_wt(1) + k + 1 ) * theta( map_wt(1) + k + 1) / & + moist_dyn_tot( map_wt(1) + k + 1) + exner( map_w3(1) + k + 1 ) = exner( map_w3 (1) + k ) & + - ( gravity - coriolis_term( map_wt(1) + k + 1 ) ) * dz & + / ( cp * theta_moist ) + + end do + +end subroutine hstat_cori_g_const_code + +end module hstat_cori_g_const_kernel_mod diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 new file mode 100644 index 000000000..4467e8b81 --- /dev/null +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -0,0 +1,226 @@ +!----------------------------------------------------------------------------- +! (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. +!----------------------------------------------------------------------------- + +!> @brief +!> @details +module regrav_geopot_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_READWRITE, & + CELL_COLUMN +use constants_mod, only: r_def, i_def +use extrusion_config_mod, only: planet_radius +use fs_continuity_mod, only: Wtheta, W3 +use kernel_mod, only: kernel_type +use planet_config_mod, only: gravity, cp + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +type, public, extends(kernel_type) :: regrav_geopot_kernel_type + private + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: regrav_geopot_code +end type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- +public :: regrav_geopot_code + +contains + +!> @param[in] nlayers Number of layers +!> @param[in,out] theta Potential temperature field +!> @param[in,out] exner Exner pressure field +!> @param[in] exner_in_wth The Exner pressure in Wtheta +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] moist_dyn_fac Water factor +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +subroutine regrav_geopot_code( & + nlayers, & + temperature, & + theta, & + exner_in_wth, & + exner, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + moist_dyn_fac, & + height_w3, & + height_wth, & + w3_mask, & + ndf_wt, & + undf_wt, & + map_wt, & + ndf_w3, & + undf_w3, & + map_w3 ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot, & + moist_dyn_fac + real(kind=r_def), dimension(undf_wt), intent(in) :: exner_in_wth, & + height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + + ! Internal variables + integer(kind=i_def) :: k + real(kind=r_def) :: temp_virt + + integer(kind=i_def) :: itn, nitns + real(kind=r_def) :: ht_wt(0:nlayers), ht_w3(nlayers) + real(kind=r_def) :: g_wt + real(kind=r_def) :: th(0:nlayers), ex(nlayers) + real(kind=r_def) :: exner_surf + real(kind=r_def) :: weight1 + + ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) + ! setting exner to 1 + if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then + do k = 0, nlayers-1 + exner(map_w3(1) + k ) = 1.0_r_def + enddo + return + end if + + do k = 0, nlayers + ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & + ( planet_radius - height_wth(map_wt(1)+k) ) + end do + + temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & + moist_dyn_tot( map_wt(1) ) + weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) + exner_surf = exner( map_w3(1) ) * ( cp * temp_virt ) / & + ( cp * temp_virt - ( gravity - coriolis_term( map_wt(1) ) ) * weight1 ) + theta( map_wt(1) ) = temperature( map_wt(1) ) / exner_surf + + +! Reuse th as temporary storage for T + do k = 0, nlayers + th(k) = temperature( map_wt(1) + k ) + end do + + ! Map temperature from newly computed grid back onto current model grid + call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & + temperature(map_wt(1)) ) + + ! Once T has been interpolated, exner needs recalculating for hydrostatic + ! balance on the current model grid. +! temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & +! moist_dyn_tot( map_wt(1) ) +! weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) +! g_wt = gravity * ( planet_radius / ( planet_radius + height_wth(map_wt(1)) ) )**2 & +! - coriolis_term( map_wt(1) ) +! +! exner( map_w3(1) ) = exner_surf * & +! ( cp * temp_virt - g_wt * weight1 ) / & +! ( cp * temp_virt ) + +end subroutine regrav_geopot_code + +subroutine interp( nin, zin, fin, nout, zout, fout ) + +implicit none + +integer(kind=i_def), intent(in) :: nin, nout +real(kind=r_def), intent(in) :: zin(nin), fin(nin), zout(nout) +real(kind=r_def), intent(inout) :: fout(nout) + +integer(kind=i_def) :: kout, kin +real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max +real(kind=r_def) :: eta_out, eta_1, eta_2 + +if ( zin(1) > zin(nin) ) then + zmax = zin(1) + zmin = zin(nin) + f_at_max = fin(1) + f_at_min = fin(nin) +else + zmax = zin(nin) + zmin = zin(1) + f_at_max = fin(nin) + f_at_min = fin(1) +end if + + +do kout = 1, nout + if ( zout(kout) <= zmin ) then + fout(kout) = f_at_min + else if ( zout(kout) >= zmax ) then + fout(kout) = f_at_max + else + do kin = 1, nin - 1 + if ( ( zout(kout) - zin(kin) ) * & + ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then + eta_out = sqrt(max(zout(kout), 0.0_r_def)) + eta_1 = sqrt(max(zin(kin), 0.0_r_def)) + eta_2 = sqrt(max(zin(kin+1), 0.0_r_def)) + eta_out = zout(kout) + eta_1 = zin(kin) + eta_2 = zin(kin+1) + xi = ( eta_out - eta_1 ) / ( eta_2 - eta_1 ) + fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) + fout(kout) = exp(fout(kout)) + end if + end do + end if +end do + +end subroutine interp + +! Does the input UM data, T(z), pi(z), have the correct z information? +! If T and pi have been interpolated onto the LFRic grid by UM2LFRic, why +! is there a problem with the lower boundary? I'm assuming this problem is +! that the orographic height of UM and LFRic are different. but if they are, +! then the UM and LFRic grids are different and the T(z), pi(z) are not +! valid at LFRic grid heights. Need to look at precisely what data is +! coming from um2lfric. + +end module regrav_geopot_kernel_mod diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 new file mode 100644 index 000000000..4128f3e1e --- /dev/null +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -0,0 +1,153 @@ +!----------------------------------------------------------------------------- +! (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. +!----------------------------------------------------------------------------- + +!> @brief Defines exner pressure from vertical balance with a deep gravity, +!> but using input data based on a shallow gravity. +!> @details Using the input exner pressure at level eos_index, integrate down +!> to the surface using hydrostatic balance, plus Coriolis terms, +!> using a discretisation using absolute temperature and constant +!> gravity g (shallow definition). Then integrate from the surface to +!> the top but using the model geopotential (may be defined as shallow +!> or deep depending on configuration). +module regrav_isotherm_kernel_mod + +use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_SCALAR, GH_INTEGER, & + GH_READ, GH_READWRITE, & + CELL_COLUMN +use constants_mod, only : r_def, i_def +use fs_continuity_mod, only : Wtheta, W3 +use kernel_mod, only : kernel_type + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +type, public, extends(kernel_type) :: regrav_isotherm_kernel_type + private + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: regrav_isotherm_code +end type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- +public :: regrav_isotherm_code + +contains + +!> @param[in] nlayers Number of layers +!> @param[in,out] exner Exner pressure field +!> @param[in] temperature Absolute temperature field +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] moist_dyn_fac Water factor +!> @param[in] phi Geopotential field +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +subroutine regrav_isotherm_code( & + nlayers, & + exner, & + temperature, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + moist_dyn_fac, & + phi, & + height_w3, & + height_wth, & + w3_mask, & + ndf_w3, & + undf_w3, & + map_w3, & + ndf_wt, & + undf_wt, & + map_wt ) + + use planet_config_mod, only: gravity, cp + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask, & + phi + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot, & + moist_dyn_fac + real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & + height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + + ! Internal variables + integer(kind=i_def) :: k + real(kind=r_def) :: temp_virtual + real(kind=r_def) :: dz, weight1, weight2 + + ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) + ! setting exner to 1 + if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then + do k = 0, nlayers-1 + exner(map_w3(1) + k ) = 1.0_r_def + enddo + return + end if + + ! Integrate from surface to top using model geopotential phi + do k = 1, nlayers-1 + + temp_virtual = moist_dyn_gas( map_wt(1) + k ) * temperature( map_wt(1) + k ) / & + moist_dyn_tot( map_wt(1) + k ) + + dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) + weight1 = ( height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) ) / dz + weight2 = ( height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) ) /dz + + ! Pi_k = Pi_k-1 * ( cp T_k-1/2 - ( (phi_k - phi_k-1) - F_k-1/2 * dz) * (1-w) ) / + ! ( cp T_k-1/2 + ( (phi_k - phi_k-1) - F_k-1/2 * dz) * w ) + exner( map_w3(1) + k ) = exner( map_w3 (1) + k-1 ) * & + ( cp * temp_virtual - ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & + - coriolis_term( map_wt(1) + k ) * dz ) * weight1 ) / & + ( cp * temp_virtual + ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & + - coriolis_term( map_wt(1) + k ) * dz ) * weight2 ) + + end do + +end subroutine regrav_isotherm_code + +end module regrav_isotherm_kernel_mod From eed41cc4681c416890d42ff1f830fd244f1fabd4 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 09:05:49 +0100 Subject: [PATCH 04/27] Geopotential height diagnostics --- .../algorithm/pres_lev_diags_alg_mod.x90 | 13 +++++++------ .../source/kernel/geo_on_pres_kernel_mod.F90 | 18 ++++++++++++------ 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index 0907509b0..b3032d320 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -12,6 +12,7 @@ module pres_lev_diags_alg_mod use io_config_mod, only: write_diag, use_xios_io use timing_mod, only: start_timing, stop_timing, tik, LPROF use constants_mod, only: r_def, i_def, l_def, str_def + use dycore_constants_mod, only: get_geopotential use initialise_diagnostics_mod, only: init_diag => init_diagnostic_field, & diag_samp => diagnostic_to_be_sampled use lfric_xios_diag_mod, only: get_axis_dimension, get_axis_values @@ -65,8 +66,8 @@ contains type( field_type ), pointer :: v_in_w3 type( field_type ), pointer :: w_in_wth type( field_type ), pointer :: wetrho_in_wth - type( field_type ), pointer :: height_w3 type( field_type ), pointer :: height_wth + type( field_type ), pointer :: geopotential ! Define fields type( field_type ) :: plev_heaviside, plev_temp, plev_u, plev_v, & @@ -261,11 +262,11 @@ contains plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', activate=plev_geopot_clim_flag) if ((plev_geopot_flag .or. plev_geopot_clim_flag) .and. use_xios_io) then - height_w3 => get_height_fv(W3, exner_w3%get_mesh_id()) - call invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, & - height_wth, exner_wth, & - nplev, plevs, plev_geopot, & - p_zero, kappa, cp, gravity, & + geopotential => get_geopotential( exner_w3%get_mesh_id() ) + call invoke_geo_on_pres_kernel_type(geopotential, exner_w3, theta_wth, & + height_wth, exner_wth, & + nplev, plevs, plev_geopot, & + p_zero, kappa, cp, gravity, & ex_power) if (plev_geopot_flag) call plev_geopot%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 index da1bf2dcf..8b10927d2 100644 --- a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 @@ -17,6 +17,7 @@ module geo_on_pres_kernel_mod use fs_continuity_mod, only: WTHETA use constants_mod, only: r_def, i_def use kernel_mod, only: kernel_type + use planet_constants_mod, only: planet_radius implicit none @@ -116,6 +117,7 @@ subroutine geo_on_pres_code(nlayers, & ! Internal variables integer(kind=i_def) :: k, level_above, top_df, kp, level_extrap real(kind=r_def) :: desired_ex, h_ref_lev + real(kind=r_def) :: geo_ht_extrap, geo_ht_lowest real(kind=r_def), parameter:: extrap_height = 2000.0_r_def do kp = 1, nplev @@ -138,9 +140,9 @@ subroutine geo_on_pres_code(nlayers, & if (level_above == -1_i_def) then ! Desired level is above model top, extrapolate up - data_out(map_out(1)+kp-1) = data_in(map_in(1)+top_df) - & - (desired_ex - ex_at_data(map_in(1)+top_df)) & - * cp * theta(map_wth(1)+nlayers) / gravity + data_out(map_out(1)+kp-1) = ( data_in(map_in(1)+top_df) - & + (desired_ex - ex_at_data(map_in(1)+top_df)) & + * cp * theta(map_wth(1)+nlayers) ) / gravity else if (level_above == 0_i_def) then ! Desired level is below surface, extrapolate down do k = 1, nlayers @@ -150,10 +152,14 @@ subroutine geo_on_pres_code(nlayers, & exit end if end do - h_ref_lev = (height_wth(map_wth(1)+level_extrap) - data_in(map_in(1))) & + geo_ht_extrap = height_wth(map_wth(1)+level_extrap) / & + ( 1.0_r_def + height_wth(map_wth(1)+level_extrap) / & + planet_radius ) + geo_ht_lowest = data_in(map_in(1)) / gravity + h_ref_lev = ( geo_ht_extrap - geo_ht_lowest ) & / (1.0_r_def - (exner_wth(map_wth(1)+level_extrap) / & ex_at_data(map_in(1)))**ex_power ) - data_out(map_out(1)+kp-1) = data_in(map_in(1)+level_above) & + data_out(map_out(1)+kp-1) = geo_ht_lowest & + h_ref_lev * (1.0_r_def - & (desired_ex/ex_at_data(map_in(1)))**ex_power) else @@ -165,7 +171,7 @@ subroutine geo_on_pres_code(nlayers, & ex_at_data(map_in(1)+level_above)) * & data_in(map_in(1)+level_above-1) ) / & (ex_at_data(map_in(1)+level_above) - & - ex_at_data(map_in(1)+level_above-1)) + ex_at_data(map_in(1)+level_above-1)) / gravity end if end do From dbc822fd3c0fc5433da901ce6cda75f64ade402a Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 09:09:39 +0100 Subject: [PATCH 05/27] Remove spurious merge --- .../map_fd_to_prognostics_alg_mod-merge.x90 | 608 ------------------ .../physics/map_fd_to_prognostics_alg_mod.x90 | 104 ++- 2 files changed, 93 insertions(+), 619 deletions(-) delete mode 100644 science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 deleted file mode 100644 index 2e60f5a3b..000000000 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod-merge.x90 +++ /dev/null @@ -1,608 +0,0 @@ -!----------------------------------------------------------------------------- -! (C) Crown copyright 2019 Met Office. All rights reserved. -! The file LICENCE, distributed with this code, contains details of the terms -! under which the code may be used. -!----------------------------------------------------------------------------- - -!> @brief Set prognostic fields from FD start dump -module map_fd_to_prognostics_alg_mod - - use constants_mod, only: r_def, i_def, l_def - use log_mod, only: log_event, & - log_scratch_space, & - LOG_LEVEL_INFO, & - LOG_LEVEL_DEBUG, & - LOG_LEVEL_TRACE - use namelist_collection_mod, only: namelist_collection_type - use namelist_mod, only: namelist_type - - ! Derived Types - use field_mod, only: field_type - use integer_field_mod, only: integer_field_type - use field_collection_mod, only: field_collection_type - use mesh_mod, only: mesh_type - use operator_mod, only: operator_type - use function_space_mod, only: function_space_type - use function_space_collection_mod, only: function_space_collection - - use quadrature_xyoz_mod, only: quadrature_xyoz_type - use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type - use sci_fem_constants_mod, only: get_rmultiplicity_fv - use sci_geometric_constants_mod, only: get_height_fv, & - get_coordinates, & - get_panel_id, & - get_da_at_w2, & - get_face_selector_ew, & - get_face_selector_ns - use sci_mapping_constants_mod, only: get_u_lon_map, & - get_u_lat_map, & - get_u_up_map, & - get_u_lon_sample, & - get_u_lat_sample, & - get_u_up_sample, & - get_w3_to_w2_displacement - use dycore_constants_mod, only: get_vert_coriolis, & - get_geopotential - use fs_continuity_mod, only: W3, W2, W2broken, Wtheta - use sci_sort_ref_kernel_mod, only: sort_ref_kernel_type - use combine_w2_field_kernel_mod, only: combine_w2_field_kernel_type - use idealised_config_mod, only: perturb_init, & - perturb_magnitude - use random_perturb_field_alg_mod, only: random_perturb_field - use mr_indices_mod, only: imr_v, imr_cl, imr_ci, imr_r, & - imr_s, imr_g - use moist_dyn_mod, only: num_moist_factors, gas_law - use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg - - ! Configuration modules - use base_mesh_config_mod, only: geometry, & - geometry_spherical, & - topology, & - topology_fully_periodic - use formulation_config_mod, only: rotating, shallow - use finite_element_config_mod, only: coord_system, coord_system_native - use initialization_config_mod, only: model_eos_height, & - read_w2h_wind, & - regravitate, & - regravitate_geopot, & - regravitate_isotherm1, & - regravitate_isotherm2, & - regravitate_isotherm3 - use physics_config_mod, only: sample_physics_winds, & - sample_physics_winds_correction - use planet_config_mod, only: gravity, p_zero, kappa, rd, cp - - implicit none - - private - public :: map_fd_to_prognostics, hydrostatic_balance, set_wind - -contains - - !> @details Setting FE prognostic fields from FD fields - !> @param[in] configuration The application configuration object - !> @param[in,out] prognostic_fields Collection of prognostic fields - !> @param[in,out] mr Array of moisture mixing ratios - !> @param[in,out] fd_fields Collection of input Finite Difference - !> fields - subroutine map_fd_to_prognostics(configuration, prognostic_fields, mr, & - moist_dyn, fd_fields) - - implicit none - - ! FE Prognostic fields - type(namelist_collection_type), intent(in) :: configuration - type(field_collection_type), intent( inout ) :: prognostic_fields - type(field_type), intent( inout ) :: mr(:) - type(field_type), intent( inout ) :: moist_dyn(num_moist_factors) - - type( field_collection_type), intent( inout ) :: fd_fields - - ! Local dereferenced fields - type( field_type ), pointer :: theta - type( field_type ), pointer :: rho - type( field_type ), pointer :: u - type( field_type ), pointer :: exner - - ! FD fields - type( field_type ), pointer :: h_wind_in_w2h - type( field_type ), pointer :: ew_wind_in_w3 - type( field_type ), pointer :: ns_wind_in_w3 - type( field_type ), pointer :: dry_rho_in_w3 - type( field_type ), pointer :: upward_wind_in_wtheta - type( field_type ), pointer :: theta_in_wtheta - type( field_type ), pointer :: mv_in_wtheta - type( field_type ), pointer :: mcl_in_wtheta - type( field_type ), pointer :: mcf_in_wtheta - type( field_type ), pointer :: mr_in_wtheta - type( field_type ), pointer :: dA - integer(i_def) :: mesh_id - - type(integer_field_type), pointer :: face_selector_ew - type(integer_field_type), pointer :: face_selector_ns - - ! FD (source) fields - if ( read_w2h_wind ) then - call fd_fields%get_field('h_wind', h_wind_in_w2h) - else - call fd_fields%get_field('ew_wind_in_w3', ew_wind_in_w3) - call fd_fields%get_field('ns_wind_in_w3', ns_wind_in_w3) - end if - - call fd_fields%get_field('dry_rho_in_w3', dry_rho_in_w3) - call fd_fields%get_field('upward_wind_in_wtheta', upward_wind_in_wtheta) - call fd_fields%get_field('theta_in_wtheta', theta_in_wtheta) - call fd_fields%get_field('mv_in_wtheta', mv_in_wtheta) - call fd_fields%get_field('mcl_in_wtheta', mcl_in_wtheta) - call fd_fields%get_field('mcf_in_wtheta', mcf_in_wtheta) - call fd_fields%get_field('mr_in_wtheta', mr_in_wtheta) - - ! Prognostic (target) fields - call prognostic_fields%get_field('u', u) - call prognostic_fields%get_field('theta', theta) - call prognostic_fields%get_field('exner', exner) - call prognostic_fields%get_field('rho', rho) - - ! Some variables should coming in from the restart file, but currently - ! they are not output by UM2LFRIC and are set to zero or a default value - ! This will need to be reviewed if this becomes a checkpoint-restart - ! routine and/or um2lfric is updated for GAL physics - call invoke( setval_x(theta, theta_in_wtheta), & - setval_x(rho, dry_rho_in_w3), & - setval_x(mr(imr_v), mv_in_wtheta), & - setval_x(mr(imr_cl), mcl_in_wtheta), & - setval_x(mr(imr_s), mcf_in_wtheta), & - setval_x(mr(imr_r), mr_in_wtheta), & - setval_c(mr(imr_ci), 0.0_r_def), & - setval_c(mr(imr_g), 0.0_r_def), & - setval_c(exner, 0.0_r_def), & - setval_c(upward_wind_in_wtheta, 0.0_r_def)) - - if ( read_w2h_wind )then - ! We read in the horizontal (normal to face) winds on the w2h points - mesh_id = u%get_mesh_id() - dA => get_da_at_w2(mesh_id) - face_selector_ew => get_face_selector_ew(mesh_id) - face_selector_ns => get_face_selector_ns(mesh_id) - - call invoke( combine_w2_field_kernel_type(u, h_wind_in_w2h, & - upward_wind_in_wtheta, & - face_selector_ew, & - face_selector_ns), & - inc_X_times_Y(u, dA) ) - else - ! We read in cell centred zonal/meridional winds - call set_wind( u, ew_wind_in_w3, ns_wind_in_w3, upward_wind_in_wtheta ) - end if - - ! Update factors for moist dynamics - call moist_dyn_factors_alg( moist_dyn, mr ) - ! Remove any static instability - call invoke( sort_ref_kernel_type(theta) ) - ! Initialize hydrostatically balanced exner field - call hydrostatic_balance( exner, rho, theta, u, moist_dyn, & - model_eos_height ) - - if ( perturb_init ) then - call log_event( "Gungho: Applying initial perturbation to theta", & - LOG_LEVEL_INFO ) - call random_perturb_field( theta, perturb_magnitude ) - end if - - ! Free up any prognostics not required - if ( read_w2h_wind ) then - call fd_fields%remove_field("h_wind") - else - call fd_fields%remove_field("ew_wind_in_w3") - call fd_fields%remove_field("ns_wind_in_w3") - end if - call fd_fields%remove_field("dry_rho_in_w3") - call fd_fields%remove_field("upward_wind_in_wtheta") - call fd_fields%remove_field("theta_in_wtheta") - call fd_fields%remove_field("mv_in_wtheta") - call fd_fields%remove_field("mcl_in_wtheta") - call fd_fields%remove_field("mcf_in_wtheta") - call fd_fields%remove_field("mr_in_wtheta") - - call log_event( "Gungho: Set prognostic fields from FD fields", LOG_LEVEL_INFO ) - - end subroutine map_fd_to_prognostics - - !> @param[in,out] u Wind (FE) - !> @param[in] u_lon Zonal Wind (FD) - !> @param[in] u_lat Meridional Wind (FD) - !> @param[in] u_up Upward Wind (FD) - subroutine set_wind(u, u_lon, u_lat, u_up) - - use sci_convert_phys_to_hdiv_kernel_mod, & - only: convert_phys_to_hdiv_kernel_type - use dg_matrix_vector_kernel_mod, only: dg_matrix_vector_kernel_type - use dg_inc_matrix_vector_kernel_mod, only: dg_inc_matrix_vector_kernel_type - use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type - use sci_mass_matrix_solver_alg_mod, only: mass_matrix_solver_alg - use matrix_vector_kernel_mod, only: matrix_vector_kernel_type - use sci_average_w2b_to_w2_kernel_mod, & - only: average_w2b_to_w2_kernel_type - use sci_w3_to_w2_correction_kernel_mod, & - only: w3_to_w2_correction_kernel_type - use timing_mod, only: start_timing, stop_timing, & - tik, LPROF - - implicit none - - type( field_type ), intent( inout ) :: u - type( field_type ), intent( in ) :: u_lat, u_lon, u_up - - type( field_type ) :: r_u - type( field_type ) :: u_broken - type( field_type ) :: u_lat_w2, u_lon_w2, u_up_w2 - type( field_type ) :: u_lat_large, u_lon_large - type( field_type ) :: u_correction - type( field_type ), pointer :: rmultiplicity_w2 - type( field_type ), pointer :: displacement - type( field_type ), pointer :: chi(:) - type( field_type ), pointer :: panel_id - type( operator_type ), pointer :: u_lon_map - type( operator_type ), pointer :: u_lat_map - type( operator_type ), pointer :: u_up_map - type( function_space_type ), pointer :: w2_fs - type( function_space_type ), pointer :: w2b_fs - - type( mesh_type ), pointer :: mesh - integer(tik) :: id - - if ( LPROF ) call start_timing( id, 'mappings.set_wind' ) - - mesh => u%get_mesh() - - if (sample_physics_winds) then - - w2b_fs => function_space_collection%get_fs(mesh, 0, 0, W2broken) - rmultiplicity_w2 => get_rmultiplicity_fv(W2, mesh%get_id()) - u_lon_map => get_u_lon_sample(mesh%get_id()) - u_lat_map => get_u_lat_sample(mesh%get_id()) - u_up_map => get_u_up_sample(mesh%get_id()) - - call u_broken%initialise( w2b_fs ) - - call invoke( name="sample_physics_winds_in_w2", & - setval_c(u, 0.0_r_def), & - dg_matrix_vector_kernel_type(u_broken, u_lon, u_lon_map), & - dg_inc_matrix_vector_kernel_type(u_broken, u_lat, u_lat_map), & - dg_inc_matrix_vector_kernel_type(u_broken, u_up, u_up_map), & - average_w2b_to_w2_kernel_type(u, u_broken, rmultiplicity_w2) ) - - if (sample_physics_winds_correction & - .and. geometry == geometry_spherical & - .and. topology == topology_fully_periodic) then - - w2_fs => u%get_function_space() - call u_lon_w2%initialise( w2_fs ) - call u_lat_w2%initialise( w2_fs ) - call u_up_w2%initialise( w2_fs ) - call u_correction%initialise( w2_fs ) - call u_lon_large%initialise( u_lon%get_function_space(), halo_depth=2 ) - call u_lat_large%initialise( u_lat%get_function_space(), halo_depth=2 ) - - displacement => get_w3_to_w2_displacement(mesh%get_id()) - chi => get_coordinates(mesh%get_id()) - panel_id => get_panel_id(mesh%get_id()) - - call invoke( setval_c(u_lon_w2, 0.0_r_def), & - setval_c(u_lat_w2, 0.0_r_def), & - setval_c(u_up_w2, 0.0_r_def), & - ! Copy u_lon and u_lat to fields with appropriate depth - setval_X(u_lon_large, u_lon), & - setval_X(u_lat_large, u_lat), & - ! Correct horizontal wind components (but not vertical) - w3_to_w2_correction_kernel_type(u_lon_w2, u_lon_large, 1, & - displacement, & - panel_id, 1), & - w3_to_w2_correction_kernel_type(u_lat_w2, u_lat_large, 1, & - displacement, & - panel_id, 1), & - ! Convert components to HDiv field - convert_phys_to_hdiv_kernel_type(u_correction, & - u_lon_w2, u_lat_w2, & - u_up_w2, chi, panel_id, & - geometry), & - ! Increment wind with correction field - inc_X_plus_Y(u, u_correction) ) - end if - - else - - call u%copy_field_properties(r_u) - - u_lon_map => get_u_lon_map( mesh%get_id() ) - u_lat_map => get_u_lat_map( mesh%get_id() ) - u_up_map => get_u_up_map( mesh%get_id() ) - - call invoke( name="galerkin_projection_from_lonlatr", & - setval_c( u, 0.0_r_def ), & - setval_c( r_u, 0.0_r_def ), & - matrix_vector_kernel_type(r_u, u_lon, u_lon_map), & - matrix_vector_kernel_type(r_u, u_lat, u_lat_map), & - matrix_vector_kernel_type(r_u, u_up, u_up_map), & - enforce_bc_kernel_type( r_u ) ) - call mass_matrix_solver_alg(u, r_u) - end if - - if ( LPROF ) call stop_timing( id, 'mappings.set_wind' ) - - end subroutine set_wind - - !> @param[in,out] exner Exner pressure - !> @param[in] rho Dry density - !> @param[in] theta Theta - !> @param[in] u Wind - !> @param[in] moist_dyn Moisture arrays for dynamics - !> @param[in] lbc_mask Optional LBC mask - subroutine hydrostatic_balance(exner, rho, theta, u, moist_dyn, & - eos_height, lbc_mask) - use dg_inc_matrix_vector_kernel_mod, & - only : dg_inc_matrix_vector_kernel_type - use sci_enforce_bc_kernel_mod, only : enforce_bc_kernel_type - use hydrostatic_coriolis_kernel_mod, & - only : hydrostatic_coriolis_kernel_type - use hstat_cori_g_const_kernel_mod, & - only : hstat_cori_g_const_kernel_type - use hydro_shallow_to_deep_kernel_mod, & - only : hydro_shallow_to_deep_kernel_type - use regrav_isotherm_kernel_mod, only : regrav_isotherm_kernel_type - use regrav_geopot_kernel_mod, only : regrav_geopot_kernel_type - use moist_dyn_mod, only : num_moist_factors - use matrix_vector_kernel_mod, only : matrix_vector_kernel_type - use sci_nodal_xyz_coordinates_kernel_mod, & - only : nodal_xyz_coordinates_kernel_type - use sci_mass_matrix_solver_alg_mod, & - only : mass_matrix_solver_alg - use sample_eos_pressure_kernel_mod, & - only : sample_eos_pressure_kernel_type - use sample_eos_rho_kernel_mod, only : sample_eos_rho_kernel_type - use sci_sample_w3_to_wtheta_kernel_mod, & - only : sample_w3_to_wtheta_kernel_type - use sci_convert_cart2sphere_vector_kernel_mod, & - only: convert_cart2sphere_vector_kernel_type - - implicit none - - type( field_type ), intent( inout ) :: exner - type( field_type ), intent( inout ) :: theta, rho, u - type( field_type ), intent( in ) :: moist_dyn(num_moist_factors) - integer(i_def), intent( in ) :: eos_height - type( field_type ), intent( in ), optional, target :: lbc_mask - - type(field_type ) :: c_vert_term, r_c - type(field_type ) :: exner_in_wth - type(field_type ) :: temperature - type(field_type), pointer :: geopotential - type(field_type), pointer :: chi(:) - type(field_type), pointer :: panel_id - type(field_type), pointer :: height_w3 - type(field_type), pointer :: height_wth - type(operator_type), pointer :: vert_coriolis - type(mesh_type), pointer :: mesh - type(field_type), pointer :: mask - type(field_type), target :: dummy_mask - integer(i_def) :: eos_level - - nullify(geopotential, chi, panel_id, height_w3, & - height_wth, vert_coriolis, mesh, mask ) - - mesh => exner%get_mesh() - height_w3 => get_height_fv( W3, mesh%get_id() ) - height_wth => get_height_fv( Wtheta, mesh%get_id() ) - geopotential => get_geopotential( mesh%get_id() ) - - ! If lbc_mask wasn't supplied, create a dummy mask - if (present(lbc_mask)) then - mask => lbc_mask - else - call dummy_mask%initialise( vector_space = exner%get_function_space() ) - call invoke( setval_c(dummy_mask, 1.0_r_def) ) - mask => dummy_mask - end if - - chi => get_coordinates( mesh%get_id() ) - panel_id => get_panel_id( mesh%get_id() ) - vert_coriolis => get_vert_coriolis( mesh%get_id() ) - - ! Create intermediate fields. - call theta%copy_field_properties(c_vert_term) - call invoke( setval_c(c_vert_term, 0.0_r_def) ) - - ! Compute contribution of Coriolis force to vertical balance - if ( associated(vert_coriolis) ) then - - call r_c%initialise( theta%get_function_space() ) - - ! Apply Coriolis mass matrix, and solve for vertical term in Wtheta - call invoke( setval_c(r_c, 0.0_r_def), & - dg_inc_matrix_vector_kernel_type(r_c, u, vert_coriolis), & - enforce_bc_kernel_type(r_c) ) - call mass_matrix_solver_alg(c_vert_term, r_c, bc_flag=.true.) - ! Enforce that result is zero at domain bottom and top - call invoke( enforce_bc_kernel_type(c_vert_term) ) - end if ! associated(vert_coriolis) - - ! Get the index of the level closest to the height, ignoring the top eta - ! level as eos_level is to be used for a W3 space rather than Wtheta - eos_level = get_nearest_level( mesh, eos_height ) - if ( eos_level == (mesh%get_nlayers() + 1) ) eos_level = eos_level - 1 - - if (shallow) then - ! Initialize the vertical profile by starting with the equation of state at - ! eos_level, and then integrating from this level to the top and the bottom - ! using hydrostatic balance with Coriolis terms. - call invoke( hydrostatic_coriolis_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn, geopotential, & - height_w3, mask, & - p_zero, kappa, & - rd, cp, & - eos_level ) ) - else - ! Initialize the vertical profile of exner pressure by initializing - ! with the equation of state, and then integrating from eos_level to - ! the surface using hydrostatic balance with Coriolis terms and - ! absolute temperature and constant gravity, followed by integrating - ! from the surface to the top with the model geopotential, followed - ! by rederiving potential temperature (theta) and density (rho). - call theta%copy_field_properties(exner_in_wth) - call theta%copy_field_properties(temperature) - - select case( regravitate ) - case( regravitate_isotherm1 ) - call invoke( & - ! Compute exner from EoS - sample_eos_pressure_kernel_type( exner, rho, & - theta, moist_dyn(gas_law) ), & - - ! Compute exner in Wtheta (required to calculate T) - sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & - height_wth, height_w3 ), & - - ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - - ! Rebalance exner in vertical, using hydrostatic balance. - ! This integrates down using g (shallow) from eos_level and - ! then integrates up using phi (deep). - hydro_shallow_to_deep_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn, geopotential, & - height_w3, height_wth, & - mask, gravity, & - cp, eos_level ) ) - case( regravitate_isotherm2 ) - call invoke( & - ! Compute exner from EoS - sample_eos_pressure_kernel_type( exner, rho, & - theta, moist_dyn(gas_law) ), & - - ! Compute exner in Wtheta (required to calculate T) - sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & - height_wth, height_w3 ), & - - ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - ! Rebalance exner in vertical, using hydrostatic balance. - ! This integrates down using g (shallow) from eos_level and - ! then integrates up using phi (deep). - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - height_w3, mask, & - p_zero, kappa, & - rd, cp, & - eos_level ), & - regrav_isotherm_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - geopotential, & - height_w3, height_wth, & - mask ) ) - case( regravitate_isotherm3 ) - call invoke( & - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - height_w3, mask, & - p_zero, kappa, & - rd, cp, & - eos_level ), & - sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & - height_wth, height_w3 ), & - - ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - regrav_isotherm_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - geopotential, & - height_w3, height_wth, & - mask ) ) - case( regravitate_geopot ) - call invoke( & - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn, & - height_w3, mask, & - p_zero, kappa, & - rd, cp, & - eos_level ), & - sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & - height_wth, height_w3 ), & - ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - regrav_geopot_kernel_type( temperature, theta, & - exner_in_wth, exner, & - c_vert_term, & - moist_dyn, & - height_w3, height_wth, & - mask ), & - regrav_isotherm_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn, geopotential, & - height_w3, height_wth, & - mask ) ) - end select - call invoke( & - ! Recompute exner in Wtheta (required to calculate theta) - sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & - height_wth, height_w3 ), & - - ! Recompute potential temperature, theta = T /exner - X_divideby_Y( theta, temperature, exner_in_wth ), & - - ! Recompute density using the equation of state - sample_eos_rho_kernel_type(rho, exner, theta, & - moist_dyn(gas_law), & - kappa, rd, p_zero) & - ) - - end if - - nullify( mesh ) - - end subroutine hydrostatic_balance - - !> @brief Returns the index of the vertical level closest to the input height - !> @details Find the index of the vertical layer that is closest to a given - !! height that is specified as an integer percentage of the total height. - !> @param[in] mesh The mesh from which to obtain the eta vertical levels - !> @param[in] height The height at which to apply the equation of state, - !> specified as a percentage of the total height of the - !> domain - !> @return level Index of vertical level nearest to specified [%] height. - function get_nearest_level(mesh, height) result(level) - - type( mesh_type ), intent( in ) :: mesh - integer(i_def), intent( in ) :: height - integer(i_def) :: level - real(r_def), allocatable :: levels(:), abs_levels(:) - integer(i_def) :: eta_level_number - - eta_level_number = mesh%get_nlayers() + 1 - allocate( levels( eta_level_number ) ) - allocate( abs_levels( eta_level_number ) ) - - call mesh%get_eta( levels ) - - ! Find the index of the eta height level that is closest to the height - ! (given as a percentage) - abs_levels = abs( levels - ( 0.01_r_def * height ) ) - level = minloc( abs_levels, dim = 1 ) - - deallocate( levels, abs_levels ) - - return - - end function get_nearest_level - -end module map_fd_to_prognostics_alg_mod diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 8f81cc672..2e60f5a3b 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -61,6 +61,13 @@ module map_fd_to_prognostics_alg_mod topology_fully_periodic use formulation_config_mod, only: rotating, shallow use finite_element_config_mod, only: coord_system, coord_system_native + use initialization_config_mod, only: model_eos_height, & + read_w2h_wind, & + regravitate, & + regravitate_geopot, & + regravitate_isotherm1, & + regravitate_isotherm2, & + regravitate_isotherm3 use physics_config_mod, only: sample_physics_winds, & sample_physics_winds_correction use planet_config_mod, only: gravity, p_zero, kappa, rd, cp @@ -114,15 +121,6 @@ contains type(integer_field_type), pointer :: face_selector_ew type(integer_field_type), pointer :: face_selector_ns - integer(i_def) :: model_eos_height - logical(l_def) :: read_w2h_wind - type(namelist_type), pointer :: initialization_nml - - initialization_nml => configuration%get_namelist('initialization') - call initialization_nml%get_value( 'model_eos_height', model_eos_height ) - call initialization_nml%get_value( 'read_w2h_wind', read_w2h_wind ) - nullify( initialization_nml ) - ! FD (source) fields if ( read_w2h_wind ) then call fd_fields%get_field('h_wind', h_wind_in_w2h) @@ -347,8 +345,12 @@ contains use sci_enforce_bc_kernel_mod, only : enforce_bc_kernel_type use hydrostatic_coriolis_kernel_mod, & only : hydrostatic_coriolis_kernel_type + use hstat_cori_g_const_kernel_mod, & + only : hstat_cori_g_const_kernel_type use hydro_shallow_to_deep_kernel_mod, & only : hydro_shallow_to_deep_kernel_type + use regrav_isotherm_kernel_mod, only : regrav_isotherm_kernel_type + use regrav_geopot_kernel_mod, only : regrav_geopot_kernel_type use moist_dyn_mod, only : num_moist_factors use matrix_vector_kernel_mod, only : matrix_vector_kernel_type use sci_nodal_xyz_coordinates_kernel_mod, & @@ -449,7 +451,10 @@ contains ! by rederiving potential temperature (theta) and density (rho). call theta%copy_field_properties(exner_in_wth) call theta%copy_field_properties(temperature) - call invoke( & + + select case( regravitate ) + case( regravitate_isotherm1 ) + call invoke( & ! Compute exner from EoS sample_eos_pressure_kernel_type( exner, rho, & theta, moist_dyn(gas_law) ), & @@ -469,8 +474,85 @@ contains moist_dyn, geopotential, & height_w3, height_wth, & mask, gravity, & - cp, eos_level ), & + cp, eos_level ) ) + case( regravitate_isotherm2 ) + call invoke( & + ! Compute exner from EoS + sample_eos_pressure_kernel_type( exner, rho, & + theta, moist_dyn(gas_law) ), & + + ! Compute exner in Wtheta (required to calculate T) + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + ! Rebalance exner in vertical, using hydrostatic balance. + ! This integrates down using g (shallow) from eos_level and + ! then integrates up using phi (deep). + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask ) ) + case( regravitate_isotherm3 ) + call invoke( & + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask ) ) + case( regravitate_geopot ) + call invoke( & + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + regrav_geopot_kernel_type( temperature, theta, & + exner_in_wth, exner, & + c_vert_term, & + moist_dyn, & + height_w3, height_wth, & + mask ), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn, geopotential, & + height_w3, height_wth, & + mask ) ) + end select + call invoke( & ! Recompute exner in Wtheta (required to calculate theta) sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & From 0413dc46c3ffcb83d28ed18e4c77cad1affa371c Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 09:23:08 +0100 Subject: [PATCH 06/27] Regravitate rose stem test --- rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc | 4 ++-- rose-stem/site/meto/groups/groups_lfric_atm.cylc | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc index a49293d58..9e0fab571 100644 --- a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc +++ b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc @@ -50,10 +50,10 @@ "panel_decomp": "auto_nonuniform", }) %} -{% elif task_ns.conf_name == "nwp_gal9_deep-C48_MG" %} +{% elif task_ns.conf_name == "nwp_gal9_regrav-C48_MG" %} {% do task_dict.update({ - "opt_confs": ["um_dump","deep","test_diags","physics_segmentation"], + "opt_confs": ["um_dump","regrav","test_diags","physics_segmentation"], "resolution": "C48_MG", "DT": 1800, "tsteps": 36, diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index ea319d080..7e73477b1 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -19,6 +19,7 @@ "lfric_atm_nwp_azspice_extra": [ "lfric_atm_nwp_gal9-C48_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_debug-C48_MG_azspice_gnu_full-debug-32bit", + "lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit", "lfric_atm_nwp_gal9_coarse_aero-C48_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_da-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_eda-C12_azspice_gnu_fast-debug-32bit", From a4f131af5b2dc2b35ba13d7fee263c2dd70cbf7e Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 10:31:03 +0100 Subject: [PATCH 07/27] Simplify argument lists --- .../physics/map_fd_to_prognostics_alg_mod.x90 | 51 +++++++++---------- .../hstat_cori_g_const_kernel_mod.F90 | 8 +-- .../regrav_isotherm_kernel_mod.F90 | 24 ++++----- 3 files changed, 37 insertions(+), 46 deletions(-) diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 2e60f5a3b..346250fed 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -495,8 +495,6 @@ contains moist_dyn(gas_law), & moist_dyn(total_mass), & height_w3, mask, & - p_zero, kappa, & - rd, cp, & eos_level ), & regrav_isotherm_kernel_type( exner, temperature, & c_vert_term, & @@ -512,44 +510,43 @@ contains moist_dyn(gas_law), & moist_dyn(total_mass), & height_w3, mask, & - p_zero, kappa, & - rd, cp, & eos_level ), & sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - regrav_isotherm_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn(gas_law), & + X_times_Y( temperature, exner_in_wth, theta), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & moist_dyn(total_mass), & - geopotential, & - height_w3, height_wth, & + geopotential, & + height_w3, height_wth, & mask ) ) case( regravitate_geopot ) call invoke( & - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn, & - height_w3, mask, & - p_zero, kappa, & - rd, cp, & - eos_level ), & + hstat_cori_g_const_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + height_w3, mask, & + eos_level ), & sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & + X_times_Y( temperature, exner_in_wth, theta), & regrav_geopot_kernel_type( temperature, theta, & - exner_in_wth, exner, & - c_vert_term, & - moist_dyn, & - height_w3, height_wth, & - mask ), & - regrav_isotherm_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn, geopotential, & - height_w3, height_wth, & + exner_in_wth, exner, & + c_vert_term, & + moist_dyn, & + height_w3, height_wth, & + mask ), & + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & mask ) ) end select call invoke( & diff --git a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 index f7809f661..f732a2126 100644 --- a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 @@ -22,7 +22,7 @@ module hstat_cori_g_const_kernel_mod use constants_mod, only : r_def, i_def use fs_continuity_mod, only : Wtheta, W3 use kernel_mod, only : kernel_type -use planet_config_mod, only : gravity +use planet_config_mod, only : gravity, p_zero, kappa, rd, cp implicit none @@ -33,7 +33,7 @@ module hstat_cori_g_const_kernel_mod !------------------------------------------------------------------------------- type, public, extends(kernel_type) :: hstat_cori_g_const_kernel_type private - type(arg_type) :: meta_args(13) = (/ & + type(arg_type) :: meta_args(9) = (/ & arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3), & arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & @@ -42,10 +42,6 @@ module hstat_cori_g_const_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & /) type(func_type) :: meta_funcs(2) = (/ & diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 index 4128f3e1e..73321c439 100644 --- a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -32,15 +32,16 @@ module regrav_isotherm_kernel_mod !------------------------------------------------------------------------------- type, public, extends(kernel_type) :: regrav_isotherm_kernel_type private - type(arg_type) :: meta_args(8) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & /) integer :: operates_on = CELL_COLUMN contains @@ -60,7 +61,6 @@ module regrav_isotherm_kernel_mod !> @param[in] coriolis_term Vertical component of the coriolis term !> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon !> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!> @param[in] moist_dyn_fac Water factor !> @param[in] phi Geopotential field !> @param[in] height_w3 Height coordinate in w3 !> @param[in] height_wth Height coordinate in wth @@ -78,7 +78,6 @@ subroutine regrav_isotherm_code( & coriolis_term, & moist_dyn_gas, & moist_dyn_tot, & - moist_dyn_fac, & phi, & height_w3, & height_wth, & @@ -108,8 +107,7 @@ subroutine regrav_isotherm_code( & w3_mask, & phi real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac + moist_dyn_tot real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & height_wth real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term From 133d1f9d11bffb8de304d2c9d1d7f9c30c698adb Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 11:12:42 +0100 Subject: [PATCH 08/27] Properly delete items from arg list --- .../initialisation/hstat_cori_g_const_kernel_mod.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 index f732a2126..da7d29f19 100644 --- a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 @@ -70,10 +70,6 @@ module hstat_cori_g_const_kernel_mod !! @param[in] moist_dyn_tot Total mass factor 1 + sum m_x !! @param[in] height_w3 Height coordinate in w3 !! @param[in] w3_mask LBC mask or Dummy mask for w3 space -!! @param[in] p_zero Reference surface pressure -!! @param[in] kappa Ratio of Rd and cp -!! @param[in] rd Gas constant for dry air -!! @param[in] cp Specific heat of dry air at constant pressure !! @param[in] eos_index Vertical level at which the equation of state !! is satisfied !! @param[in] ndf_w3 Number of degrees of freedom per cell for w3 @@ -93,10 +89,6 @@ subroutine hstat_cori_g_const_code( nlayers, & moist_dyn_tot, & height_w3, & w3_mask, & - p_zero, & - kappa, & - rd, & - cp, & eos_index, & ndf_w3, & undf_w3, & @@ -129,10 +121,6 @@ subroutine hstat_cori_g_const_code( nlayers, & real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term real(kind=r_def), dimension(1,ndf_w3,ndf_w3), intent(in) :: basis_w3 real(kind=r_def), dimension(1,ndf_wt,ndf_w3), intent(in) :: basis_wt - real(kind=r_def), intent(in) :: p_zero - real(kind=r_def), intent(in) :: kappa - real(kind=r_def), intent(in) :: rd - real(kind=r_def), intent(in) :: cp ! Internal variables integer(kind=i_def) :: k, df, dft, df3, eos_index_m1 From 4f35a7eacbf61c41c6c724949878978a6092292b Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 11:40:50 +0100 Subject: [PATCH 09/27] Missing declaration --- .../algorithm/physics/map_fd_to_prognostics_alg_mod.x90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 346250fed..7861f2882 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -51,7 +51,9 @@ module map_fd_to_prognostics_alg_mod use random_perturb_field_alg_mod, only: random_perturb_field use mr_indices_mod, only: imr_v, imr_cl, imr_ci, imr_r, & imr_s, imr_g - use moist_dyn_mod, only: num_moist_factors, gas_law + use moist_dyn_mod, only: num_moist_factors, & + gas_law, & + total_mass use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg ! Configuration modules From ddc704153f172ac1a6f9d1c7916308ecc68b13b1 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 15:36:38 +0100 Subject: [PATCH 10/27] KGO and filename correction --- rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc | 2 +- ...p_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) create mode 100644 rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt diff --git a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc index 9e0fab571..02837ff0b 100644 --- a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc +++ b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc @@ -59,7 +59,7 @@ "tsteps": 36, "mpi_parts": 24, "kgo_checks": ["checksum"], - "plot_str": "plot_map.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR", + "plot_str": "plot_map.py $NODAL_DATA_DIR/lfric_diagnostics.nc $PLOT_DIR", }) %} {% elif task_ns.conf_name == "nwp_gal9_nomg-C48" %} diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt new file mode 100644 index 000000000..652a56c0d --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 411AFEB6C87F3787 +Inner product checksum theta = 42718B765093CD28 +Inner product checksum u = 4552C9D6D84CAF9C +Inner product checksum mr1 = 40399BB2B2F333E2 +Inner product checksum mr2 = 3F39BBA239439B78 +Inner product checksum mr3 = 3EF535A3F74F3252 +Inner product checksum mr4 = 3F2BE0BCF15A0966 +Inner product checksum mr5 = 0 +Inner product checksum mr6 = 0 From 85d473c43b0dcf857aa62991cbb66cea6bcb228f Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 7 May 2026 15:45:04 +0100 Subject: [PATCH 11/27] Upgrade macro correction --- science/gungho/rose-meta/lfric-gungho/versions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/science/gungho/rose-meta/lfric-gungho/versions.py b/science/gungho/rose-meta/lfric-gungho/versions.py index e10b1524d..6667b0aea 100644 --- a/science/gungho/rose-meta/lfric-gungho/versions.py +++ b/science/gungho/rose-meta/lfric-gungho/versions.py @@ -21,7 +21,7 @@ def __repr__(self): class vn31_t218(MacroUpgrade): """Upgrade macro for issue #218 by Chris Smith.""" - BEFORE_TAG = "vn3.1_t135" + BEFORE_TAG = "vn3.1" AFTER_TAG = "vn3.1_t218" def upgrade(self, config, meta_config=None): From e9855a3c97b2631337369c4fdc499d537598b616 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Fri, 8 May 2026 10:27:49 +0100 Subject: [PATCH 12/27] Remove unused variables --- .../lfric-gungho/HEAD/rose-meta.conf | 4 +- .../physics/map_fd_to_prognostics_alg_mod.x90 | 5 +- .../regrav_geopot_kernel_mod.F90 | 147 ++++++++---------- .../regrav_isotherm_kernel_mod.F90 | 78 +++++----- 4 files changed, 108 insertions(+), 126 deletions(-) diff --git a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf index 7e5d04fed..066513127 100644 --- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf +++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf @@ -1971,7 +1971,7 @@ help=When true, the gravitational field is constant, as would be the under the s =from planet centre. !kind=default sort-key=Panel-A04 -trigger=namelist:initialization=regravitate: this == .false.; +trigger=namelist:initialization=regravitate: this == ".false."; type=logical [namelist:formulation=si_momentum_equation] @@ -3236,7 +3236,7 @@ help=Method for switching gravity from constant to height-varying. = back onto the original model levels. =Geopotential remap: = Physical height and geopotential height are the same for data in the - = start dump, due to gravity being constant. + = start dump, due to gravity being constant. !kind=default sort-key=Panel-A04b value-titles=Geopotential remap, Isothermal vn1, Isothermal vn2, Isothermal vn3 diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 7861f2882..8049e99e9 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -538,9 +538,10 @@ contains ! Compute absolute temperature T = theta * exner X_times_Y( temperature, exner_in_wth, theta), & regrav_geopot_kernel_type( temperature, theta, & - exner_in_wth, exner, & + exner, & c_vert_term, & - moist_dyn, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & height_w3, height_wth, & mask ), & regrav_isotherm_kernel_type( exner, temperature, & diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 index 4467e8b81..d3ec0a2a5 100644 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -27,16 +27,16 @@ module regrav_geopot_kernel_mod !------------------------------------------------------------------------------- type, public, extends(kernel_type) :: regrav_geopot_kernel_type private - type(arg_type) :: meta_args(9) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & /) integer :: operates_on = CELL_COLUMN contains @@ -50,73 +50,65 @@ module regrav_geopot_kernel_mod contains -!> @param[in] nlayers Number of layers -!> @param[in,out] theta Potential temperature field -!> @param[in,out] exner Exner pressure field -!> @param[in] exner_in_wth The Exner pressure in Wtheta -!> @param[in] coriolis_term Vertical component of the coriolis term -!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!> @param[in] moist_dyn_fac Water factor -!> @param[in] height_w3 Height coordinate in w3 -!> @param[in] height_wth Height coordinate in wth -!> @param[in] w3_mask LBC mask or Dummy mask for w3 space -!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!> @param[in] undf_wt Total number of degrees of freedom for wtheta -!> @param[in] map_wt Dofmap for the cell at column base for wt -!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!> @param[in] undf_w3 Total number of degrees of freedom for w3 -!> @param[in] map_w3 Dofmap for the cell at column base for w3 -subroutine regrav_geopot_code( & - nlayers, & - temperature, & - theta, & - exner_in_wth, & - exner, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac, & - height_w3, & - height_wth, & - w3_mask, & - ndf_wt, & - undf_wt, & - map_wt, & - ndf_w3, & - undf_w3, & - map_w3 ) +!> @param[in] nlayers Number of layers +!> @param[in,out] temperature Absolute temperature field +!> @param[in,out] theta Potential temperature field +!> @param[in,out] exner Exner pressure field +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +subroutine regrav_geopot_code( nlayers, & + temperature, & + theta, & + exner, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + height_w3, & + height_wth, & + w3_mask, & + ndf_wt, & + undf_wt, & + map_wt, & + ndf_w3, & + undf_w3, & + map_w3 ) implicit none ! Arguments - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - - real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac - real(kind=r_def), dimension(undf_wt), intent(in) :: exner_in_wth, & - height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot + real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term ! Internal variables - integer(kind=i_def) :: k - real(kind=r_def) :: temp_virt - - integer(kind=i_def) :: itn, nitns - real(kind=r_def) :: ht_wt(0:nlayers), ht_w3(nlayers) - real(kind=r_def) :: g_wt - real(kind=r_def) :: th(0:nlayers), ex(nlayers) + integer(kind=i_def) :: k + real(kind=r_def) :: temp_virt + real(kind=r_def) :: ht_wt(0:nlayers) + real(kind=r_def) :: th(0:nlayers) real(kind=r_def) :: exner_surf real(kind=r_def) :: weight1 @@ -129,6 +121,7 @@ subroutine regrav_geopot_code( & return end if + ! Geopotential height do k = 0, nlayers ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & ( planet_radius - height_wth(map_wt(1)+k) ) @@ -151,18 +144,6 @@ subroutine regrav_geopot_code( & call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & temperature(map_wt(1)) ) - ! Once T has been interpolated, exner needs recalculating for hydrostatic - ! balance on the current model grid. -! temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & -! moist_dyn_tot( map_wt(1) ) -! weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) -! g_wt = gravity * ( planet_radius / ( planet_radius + height_wth(map_wt(1)) ) )**2 & -! - coriolis_term( map_wt(1) ) -! -! exner( map_w3(1) ) = exner_surf * & -! ( cp * temp_virt - g_wt * weight1 ) / & -! ( cp * temp_virt ) - end subroutine regrav_geopot_code subroutine interp( nin, zin, fin, nout, zout, fout ) diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 index 73321c439..7e5f3f9e4 100644 --- a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -71,51 +71,51 @@ module regrav_isotherm_kernel_mod !> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta !> @param[in] undf_wt Total number of degrees of freedom for wtheta !> @param[in] map_wt Dofmap for the cell at column base for wt -subroutine regrav_isotherm_code( & - nlayers, & - exner, & - temperature, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - phi, & - height_w3, & - height_wth, & - w3_mask, & - ndf_w3, & - undf_w3, & - map_w3, & - ndf_wt, & - undf_wt, & - map_wt ) - - use planet_config_mod, only: gravity, cp +subroutine regrav_isotherm_code( nlayers, & + exner, & + temperature, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + phi, & + height_w3, & + height_wth, & + w3_mask, & + ndf_w3, & + undf_w3, & + map_w3, & + ndf_wt, & + undf_wt, & + map_wt ) + + use planet_config_mod, only: gravity, cp implicit none ! Arguments - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask, & - phi - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot - real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & - height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask, & + phi + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot + real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & + height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term ! Internal variables - integer(kind=i_def) :: k - real(kind=r_def) :: temp_virtual - real(kind=r_def) :: dz, weight1, weight2 + integer(kind=i_def) :: k + real(kind=r_def) :: temp_virtual + real(kind=r_def) :: dz, weight1, weight2 ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) ! setting exner to 1 From 313680b9f3351bb66d8daf9e517e6594a94096b7 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Fri, 8 May 2026 10:34:25 +0100 Subject: [PATCH 13/27] Descriptions --- .../source/kernel/initialisation/regrav_geopot_kernel_mod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 index d3ec0a2a5..de4f4b702 100644 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -4,8 +4,8 @@ ! under which the code may be used. !----------------------------------------------------------------------------- -!> @brief -!> @details +!> @brief Need to say something here +!> @details More to be said here module regrav_geopot_kernel_mod use argument_mod, only: arg_type, & From 36c21ee059d383f6d357f32881f82e9142710403 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Mon, 11 May 2026 11:40:21 +0100 Subject: [PATCH 14/27] Things caught by developer tests --- .../algorithm/pres_lev_diags_alg_mod.x90 | 5 +- .../source/kernel/geo_on_pres_kernel_mod.F90 | 6 +- .../source/psy/psykal_lite_phys_mod.F90 | 7 +- .../regrav_geopot_kernel_mod.F90 | 207 --------------- .../regrav_isotherm_kernel_mod.F90 | 237 +++++++++++------- 5 files changed, 158 insertions(+), 304 deletions(-) delete mode 100644 science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index b3032d320..4f4c1fe92 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -21,7 +21,8 @@ module pres_lev_diags_alg_mod use fs_continuity_mod, only: W3, WTheta use mr_indices_mod, only: nummr, imr_v use moist_dyn_mod, only: num_moist_factors, total_mass - use planet_config_mod, only: p_zero, kappa, gravity, cp + use planet_config_mod, only: p_zero, kappa, gravity, cp, & + planet_radius use planet_constants_mod, only: ex_power implicit none @@ -267,7 +268,7 @@ contains height_wth, exner_wth, & nplev, plevs, plev_geopot, & p_zero, kappa, cp, gravity, & - ex_power) + planet_radius, ex_power) if (plev_geopot_flag) call plev_geopot%write_field() if (plev_geopot_clim_flag) then diff --git a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 index 8b10927d2..c045fa94f 100644 --- a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 @@ -17,7 +17,6 @@ module geo_on_pres_kernel_mod use fs_continuity_mod, only: WTHETA use constants_mod, only: r_def, i_def use kernel_mod, only: kernel_type - use planet_constants_mod, only: planet_radius implicit none @@ -66,6 +65,7 @@ module geo_on_pres_kernel_mod !> @param[in] kappa Rd / cp !> @param[in] cp Specific heat at constant pressure !> @param[in] gravity Acceleration due to gravity + !> @param[in] planet_radius Radius of the planet !> @param[in] ex_power (cp * lapse_rate ) / g !> @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 @@ -86,6 +86,7 @@ subroutine geo_on_pres_code(nlayers, & kappa, & cp, & gravity, & + planet_radius, & ex_power, & ndf_in, undf_in, map_in, & ndf_wth, undf_wth, map_wth, & @@ -111,7 +112,8 @@ subroutine geo_on_pres_code(nlayers, & real(kind=r_def), intent(inout), dimension(undf_out) :: data_out ! Constants passed explicitly from algorithm - real(kind=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power + real(kind=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power, & + planet_radius real(kind=r_def), intent(in), dimension(nplev) :: plevs ! Internal variables 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 a064523c8..2f19b600b 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 @@ -992,10 +992,11 @@ END SUBROUTINE invoke_pres_interp_kernel_type !> Hence this module could be removed once the PSyclone ticket is !> completed SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height_wth, exner_wth, nplev, plevs, plev_geopot, & -&p_zero, kappa, cp, gravity, ex_power) +&p_zero, kappa, cp, gravity, planet_radius, ex_power) USE geo_on_pres_kernel_mod, ONLY: geo_on_pres_code USE mesh_mod, ONLY: mesh_type - REAL(KIND=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power + REAL(KIND=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power, & + planet_radius INTEGER(KIND=i_def), intent(in) :: nplev REAL(KIND=r_def), intent(in) :: plevs(nplev) TYPE(field_type), intent(in) :: height_w3, exner_w3, theta_wth, height_wth, exner_wth, plev_geopot @@ -1071,7 +1072,7 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height DO cell=loop0_start,loop0_stop ! CALL geo_on_pres_code(nlayers, height_w3_data, exner_w3_data, theta_wth_data, height_wth_data, exner_wth_data, nplev, & -&plevs, plev_geopot_data, p_zero, kappa, cp, gravity, ex_power, ndf_adspc1_height_w3, undf_adspc1_height_w3, & +&plevs, plev_geopot_data, p_zero, kappa, cp, gravity, planet_radius, ex_power, ndf_adspc1_height_w3, undf_adspc1_height_w3, & &map_adspc1_height_w3(:,cell), ndf_wtheta, undf_wtheta, map_wtheta(:,cell), ndf_adspc2_plev_geopot, undf_adspc2_plev_geopot, & &map_adspc2_plev_geopot(:,cell)) END DO diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 deleted file mode 100644 index de4f4b702..000000000 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ /dev/null @@ -1,207 +0,0 @@ -!----------------------------------------------------------------------------- -! (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. -!----------------------------------------------------------------------------- - -!> @brief Need to say something here -!> @details More to be said here -module regrav_geopot_kernel_mod - -use argument_mod, only: arg_type, & - GH_FIELD, GH_REAL, & - GH_READ, GH_READWRITE, & - CELL_COLUMN -use constants_mod, only: r_def, i_def -use extrusion_config_mod, only: planet_radius -use fs_continuity_mod, only: Wtheta, W3 -use kernel_mod, only: kernel_type -use planet_config_mod, only: gravity, cp - -implicit none - -private - -!------------------------------------------------------------------------------- -! Public types -!------------------------------------------------------------------------------- -type, public, extends(kernel_type) :: regrav_geopot_kernel_type - private - type(arg_type) :: meta_args(9) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & - /) - integer :: operates_on = CELL_COLUMN -contains - procedure, nopass :: regrav_geopot_code -end type - -!------------------------------------------------------------------------------- -! Contained functions/subroutines -!------------------------------------------------------------------------------- -public :: regrav_geopot_code - -contains - -!> @param[in] nlayers Number of layers -!> @param[in,out] temperature Absolute temperature field -!> @param[in,out] theta Potential temperature field -!> @param[in,out] exner Exner pressure field -!> @param[in] coriolis_term Vertical component of the coriolis term -!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!> @param[in] height_w3 Height coordinate in w3 -!> @param[in] height_wth Height coordinate in wth -!> @param[in] w3_mask LBC mask or Dummy mask for w3 space -!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!> @param[in] undf_wt Total number of degrees of freedom for wtheta -!> @param[in] map_wt Dofmap for the cell at column base for wt -!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!> @param[in] undf_w3 Total number of degrees of freedom for w3 -!> @param[in] map_w3 Dofmap for the cell at column base for w3 -subroutine regrav_geopot_code( nlayers, & - temperature, & - theta, & - exner, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - height_w3, & - height_wth, & - w3_mask, & - ndf_wt, & - undf_wt, & - map_wt, & - ndf_w3, & - undf_w3, & - map_w3 ) - - implicit none - - ! Arguments - real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot - real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term - - ! Internal variables - integer(kind=i_def) :: k - real(kind=r_def) :: temp_virt - real(kind=r_def) :: ht_wt(0:nlayers) - real(kind=r_def) :: th(0:nlayers) - real(kind=r_def) :: exner_surf - real(kind=r_def) :: weight1 - - ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) - ! setting exner to 1 - if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then - do k = 0, nlayers-1 - exner(map_w3(1) + k ) = 1.0_r_def - enddo - return - end if - - ! Geopotential height - do k = 0, nlayers - ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & - ( planet_radius - height_wth(map_wt(1)+k) ) - end do - - temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & - moist_dyn_tot( map_wt(1) ) - weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) - exner_surf = exner( map_w3(1) ) * ( cp * temp_virt ) / & - ( cp * temp_virt - ( gravity - coriolis_term( map_wt(1) ) ) * weight1 ) - theta( map_wt(1) ) = temperature( map_wt(1) ) / exner_surf - - -! Reuse th as temporary storage for T - do k = 0, nlayers - th(k) = temperature( map_wt(1) + k ) - end do - - ! Map temperature from newly computed grid back onto current model grid - call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & - temperature(map_wt(1)) ) - -end subroutine regrav_geopot_code - -subroutine interp( nin, zin, fin, nout, zout, fout ) - -implicit none - -integer(kind=i_def), intent(in) :: nin, nout -real(kind=r_def), intent(in) :: zin(nin), fin(nin), zout(nout) -real(kind=r_def), intent(inout) :: fout(nout) - -integer(kind=i_def) :: kout, kin -real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max -real(kind=r_def) :: eta_out, eta_1, eta_2 - -if ( zin(1) > zin(nin) ) then - zmax = zin(1) - zmin = zin(nin) - f_at_max = fin(1) - f_at_min = fin(nin) -else - zmax = zin(nin) - zmin = zin(1) - f_at_max = fin(nin) - f_at_min = fin(1) -end if - - -do kout = 1, nout - if ( zout(kout) <= zmin ) then - fout(kout) = f_at_min - else if ( zout(kout) >= zmax ) then - fout(kout) = f_at_max - else - do kin = 1, nin - 1 - if ( ( zout(kout) - zin(kin) ) * & - ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then - eta_out = sqrt(max(zout(kout), 0.0_r_def)) - eta_1 = sqrt(max(zin(kin), 0.0_r_def)) - eta_2 = sqrt(max(zin(kin+1), 0.0_r_def)) - eta_out = zout(kout) - eta_1 = zin(kin) - eta_2 = zin(kin+1) - xi = ( eta_out - eta_1 ) / ( eta_2 - eta_1 ) - fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) - fout(kout) = exp(fout(kout)) - end if - end do - end if -end do - -end subroutine interp - -! Does the input UM data, T(z), pi(z), have the correct z information? -! If T and pi have been interpolated onto the LFRic grid by UM2LFRic, why -! is there a problem with the lower boundary? I'm assuming this problem is -! that the orographic height of UM and LFRic are different. but if they are, -! then the UM and LFRic grids are different and the T(z), pi(z) are not -! valid at LFRic grid heights. Need to look at precisely what data is -! coming from um2lfric. - -end module regrav_geopot_kernel_mod diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 index 7e5f3f9e4..00ed128fe 100644 --- a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -4,24 +4,19 @@ ! under which the code may be used. !----------------------------------------------------------------------------- -!> @brief Defines exner pressure from vertical balance with a deep gravity, -!> but using input data based on a shallow gravity. -!> @details Using the input exner pressure at level eos_index, integrate down -!> to the surface using hydrostatic balance, plus Coriolis terms, -!> using a discretisation using absolute temperature and constant -!> gravity g (shallow definition). Then integrate from the surface to -!> the top but using the model geopotential (may be defined as shallow -!> or deep depending on configuration). -module regrav_isotherm_kernel_mod - -use argument_mod, only : arg_type, & - GH_FIELD, GH_REAL, & - GH_SCALAR, GH_INTEGER, & - GH_READ, GH_READWRITE, & - CELL_COLUMN -use constants_mod, only : r_def, i_def -use fs_continuity_mod, only : Wtheta, W3 -use kernel_mod, only : kernel_type +!> @brief Need to say something here +!> @details More to be said here +module regrav_geopot_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_READWRITE, & + CELL_COLUMN +use constants_mod, only: r_def, i_def +use extrusion_config_mod, only: planet_radius +use fs_continuity_mod, only: Wtheta, W3 +use kernel_mod, only: kernel_type +use planet_config_mod, only: gravity, cp implicit none @@ -30,71 +25,67 @@ module regrav_isotherm_kernel_mod !------------------------------------------------------------------------------- ! Public types !------------------------------------------------------------------------------- -type, public, extends(kernel_type) :: regrav_isotherm_kernel_type +type, public, extends(kernel_type) :: regrav_geopot_kernel_type private - type(arg_type) :: meta_args(9) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & /) integer :: operates_on = CELL_COLUMN contains - procedure, nopass :: regrav_isotherm_code + procedure, nopass :: regrav_geopot_code end type !------------------------------------------------------------------------------- ! Contained functions/subroutines !------------------------------------------------------------------------------- -public :: regrav_isotherm_code +public :: regrav_geopot_code contains -!> @param[in] nlayers Number of layers -!> @param[in,out] exner Exner pressure field -!> @param[in] temperature Absolute temperature field -!> @param[in] coriolis_term Vertical component of the coriolis term -!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!> @param[in] phi Geopotential field -!> @param[in] height_w3 Height coordinate in w3 -!> @param[in] height_wth Height coordinate in wth -!> @param[in] w3_mask LBC mask or Dummy mask for w3 space -!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!> @param[in] undf_w3 Total number of degrees of freedom for w3 -!> @param[in] map_w3 Dofmap for the cell at column base for w3 -!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!> @param[in] undf_wt Total number of degrees of freedom for wtheta -!> @param[in] map_wt Dofmap for the cell at column base for wt -subroutine regrav_isotherm_code( nlayers, & - exner, & - temperature, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - phi, & - height_w3, & - height_wth, & - w3_mask, & - ndf_w3, & - undf_w3, & - map_w3, & - ndf_wt, & - undf_wt, & - map_wt ) - - use planet_config_mod, only: gravity, cp +!> @param[in] nlayers Number of layers +!> @param[in,out] temperature Absolute temperature field +!> @param[in,out] theta Potential temperature field +!> @param[in,out] exner Exner pressure field +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +subroutine regrav_geopot_code( nlayers, & + temperature, & + theta, & + exner, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + height_w3, & + height_wth, & + w3_mask, & + ndf_wt, & + undf_wt, & + map_wt, & + ndf_w3, & + undf_w3, & + map_w3 ) implicit none ! Arguments - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - integer(kind=i_def), intent(in) :: nlayers, & ndf_w3, & undf_w3, & @@ -103,19 +94,24 @@ subroutine regrav_isotherm_code( nlayers, & integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask, & - phi + w3_mask real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & moist_dyn_tot - real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & - height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term ! Internal variables integer(kind=i_def) :: k - real(kind=r_def) :: temp_virtual - real(kind=r_def) :: dz, weight1, weight2 + real(kind=r_def) :: temp_virt + real(kind=r_def) :: ht_wt(0:nlayers) + real(kind=r_def) :: th(0:nlayers) + real(kind=r_def) :: exner_surf + real(kind=r_def) :: weight1 ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) ! setting exner to 1 @@ -126,26 +122,87 @@ subroutine regrav_isotherm_code( nlayers, & return end if - ! Integrate from surface to top using model geopotential phi - do k = 1, nlayers-1 - - temp_virtual = moist_dyn_gas( map_wt(1) + k ) * temperature( map_wt(1) + k ) / & - moist_dyn_tot( map_wt(1) + k ) + ! Geopotential height + do k = 0, nlayers + ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & + ( planet_radius - height_wth(map_wt(1)+k) ) + end do - dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) - weight1 = ( height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) ) / dz - weight2 = ( height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) ) /dz + temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & + moist_dyn_tot( map_wt(1) ) + weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) + exner_surf = exner( map_w3(1) ) * ( cp * temp_virt ) / & + ( cp * temp_virt - ( gravity - coriolis_term( map_wt(1) ) ) * weight1 ) + theta( map_wt(1) ) = temperature( map_wt(1) ) / exner_surf - ! Pi_k = Pi_k-1 * ( cp T_k-1/2 - ( (phi_k - phi_k-1) - F_k-1/2 * dz) * (1-w) ) / - ! ( cp T_k-1/2 + ( (phi_k - phi_k-1) - F_k-1/2 * dz) * w ) - exner( map_w3(1) + k ) = exner( map_w3 (1) + k-1 ) * & - ( cp * temp_virtual - ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & - - coriolis_term( map_wt(1) + k ) * dz ) * weight1 ) / & - ( cp * temp_virtual + ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & - - coriolis_term( map_wt(1) + k ) * dz ) * weight2 ) +! Reuse th as temporary storage for T + do k = 0, nlayers + th(k) = temperature( map_wt(1) + k ) end do -end subroutine regrav_isotherm_code + ! Map temperature from newly computed grid back onto current model grid + call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & + temperature(map_wt(1)) ) + +end subroutine regrav_geopot_code + +subroutine interp( nin, zin, fin, nout, zout, fout ) + +implicit none + +integer(kind=i_def), intent(in) :: nin, nout +real(kind=r_def), intent(in) :: zin(nin), fin(nin), zout(nout) +real(kind=r_def), intent(inout) :: fout(nout) + +integer(kind=i_def) :: kout, kin +real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max +real(kind=r_def) :: eta_out, eta_1, eta_2 + +if ( zin(1) > zin(nin) ) then + zmax = zin(1) + zmin = zin(nin) + f_at_max = fin(1) + f_at_min = fin(nin) +else + zmax = zin(nin) + zmin = zin(1) + f_at_max = fin(nin) + f_at_min = fin(1) +end if + + +do kout = 1, nout + if ( zout(kout) <= zmin ) then + fout(kout) = f_at_min + else if ( zout(kout) >= zmax ) then + fout(kout) = f_at_max + else + do kin = 1, nin - 1 + if ( ( zout(kout) - zin(kin) ) * & + ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then + eta_out = sqrt(max(zout(kout), 0.0_r_def)) + eta_1 = sqrt(max(zin(kin), 0.0_r_def)) + eta_2 = sqrt(max(zin(kin+1), 0.0_r_def)) + eta_out = zout(kout) + eta_1 = zin(kin) + eta_2 = zin(kin+1) + xi = ( eta_out - eta_1 ) / ( eta_2 - eta_1 ) + fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) + fout(kout) = exp(fout(kout)) + end if + end do + end if +end do + +end subroutine interp + +! Does the input UM data, T(z), pi(z), have the correct z information? +! If T and pi have been interpolated onto the LFRic grid by UM2LFRic, why +! is there a problem with the lower boundary? I'm assuming this problem is +! that the orographic height of UM and LFRic are different. but if they are, +! then the UM and LFRic grids are different and the T(z), pi(z) are not +! valid at LFRic grid heights. Need to look at precisely what data is +! coming from um2lfric. -end module regrav_isotherm_kernel_mod +end module regrav_geopot_kernel_mod From 60cf96e9f3deaa621a334ac258735655251225cd Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Mon, 11 May 2026 12:10:20 +0100 Subject: [PATCH 15/27] Undo messing up of files --- .../regrav_geopot_kernel_mod.F90 | 208 +++++++++++++++ .../regrav_isotherm_kernel_mod.F90 | 245 +++++++----------- 2 files changed, 302 insertions(+), 151 deletions(-) create mode 100644 science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 new file mode 100644 index 000000000..00ed128fe --- /dev/null +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -0,0 +1,208 @@ +!----------------------------------------------------------------------------- +! (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. +!----------------------------------------------------------------------------- + +!> @brief Need to say something here +!> @details More to be said here +module regrav_geopot_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_READWRITE, & + CELL_COLUMN +use constants_mod, only: r_def, i_def +use extrusion_config_mod, only: planet_radius +use fs_continuity_mod, only: Wtheta, W3 +use kernel_mod, only: kernel_type +use planet_config_mod, only: gravity, cp + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +type, public, extends(kernel_type) :: regrav_geopot_kernel_type + private + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: regrav_geopot_code +end type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- +public :: regrav_geopot_code + +contains + +!> @param[in] nlayers Number of layers +!> @param[in,out] temperature Absolute temperature field +!> @param[in,out] theta Potential temperature field +!> @param[in,out] exner Exner pressure field +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +subroutine regrav_geopot_code( nlayers, & + temperature, & + theta, & + exner, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + height_w3, & + height_wth, & + w3_mask, & + ndf_wt, & + undf_wt, & + map_wt, & + ndf_w3, & + undf_w3, & + map_w3 ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + + + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot + real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + + ! Internal variables + integer(kind=i_def) :: k + real(kind=r_def) :: temp_virt + real(kind=r_def) :: ht_wt(0:nlayers) + real(kind=r_def) :: th(0:nlayers) + real(kind=r_def) :: exner_surf + real(kind=r_def) :: weight1 + + ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) + ! setting exner to 1 + if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then + do k = 0, nlayers-1 + exner(map_w3(1) + k ) = 1.0_r_def + enddo + return + end if + + ! Geopotential height + do k = 0, nlayers + ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & + ( planet_radius - height_wth(map_wt(1)+k) ) + end do + + temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & + moist_dyn_tot( map_wt(1) ) + weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) + exner_surf = exner( map_w3(1) ) * ( cp * temp_virt ) / & + ( cp * temp_virt - ( gravity - coriolis_term( map_wt(1) ) ) * weight1 ) + theta( map_wt(1) ) = temperature( map_wt(1) ) / exner_surf + + +! Reuse th as temporary storage for T + do k = 0, nlayers + th(k) = temperature( map_wt(1) + k ) + end do + + ! Map temperature from newly computed grid back onto current model grid + call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & + temperature(map_wt(1)) ) + +end subroutine regrav_geopot_code + +subroutine interp( nin, zin, fin, nout, zout, fout ) + +implicit none + +integer(kind=i_def), intent(in) :: nin, nout +real(kind=r_def), intent(in) :: zin(nin), fin(nin), zout(nout) +real(kind=r_def), intent(inout) :: fout(nout) + +integer(kind=i_def) :: kout, kin +real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max +real(kind=r_def) :: eta_out, eta_1, eta_2 + +if ( zin(1) > zin(nin) ) then + zmax = zin(1) + zmin = zin(nin) + f_at_max = fin(1) + f_at_min = fin(nin) +else + zmax = zin(nin) + zmin = zin(1) + f_at_max = fin(nin) + f_at_min = fin(1) +end if + + +do kout = 1, nout + if ( zout(kout) <= zmin ) then + fout(kout) = f_at_min + else if ( zout(kout) >= zmax ) then + fout(kout) = f_at_max + else + do kin = 1, nin - 1 + if ( ( zout(kout) - zin(kin) ) * & + ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then + eta_out = sqrt(max(zout(kout), 0.0_r_def)) + eta_1 = sqrt(max(zin(kin), 0.0_r_def)) + eta_2 = sqrt(max(zin(kin+1), 0.0_r_def)) + eta_out = zout(kout) + eta_1 = zin(kin) + eta_2 = zin(kin+1) + xi = ( eta_out - eta_1 ) / ( eta_2 - eta_1 ) + fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) + fout(kout) = exp(fout(kout)) + end if + end do + end if +end do + +end subroutine interp + +! Does the input UM data, T(z), pi(z), have the correct z information? +! If T and pi have been interpolated onto the LFRic grid by UM2LFRic, why +! is there a problem with the lower boundary? I'm assuming this problem is +! that the orographic height of UM and LFRic are different. but if they are, +! then the UM and LFRic grids are different and the T(z), pi(z) are not +! valid at LFRic grid heights. Need to look at precisely what data is +! coming from um2lfric. + +end module regrav_geopot_kernel_mod diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 index 00ed128fe..7ba70c7aa 100644 --- a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -4,19 +4,24 @@ ! under which the code may be used. !----------------------------------------------------------------------------- -!> @brief Need to say something here -!> @details More to be said here -module regrav_geopot_kernel_mod - -use argument_mod, only: arg_type, & - GH_FIELD, GH_REAL, & - GH_READ, GH_READWRITE, & - CELL_COLUMN -use constants_mod, only: r_def, i_def -use extrusion_config_mod, only: planet_radius -use fs_continuity_mod, only: Wtheta, W3 -use kernel_mod, only: kernel_type -use planet_config_mod, only: gravity, cp +!> @brief Defines exner pressure from vertical balance with a deep gravity, +!> but using input data based on a shallow gravity. +!> @details Using the input exner pressure at level eos_index, integrate down +!> to the surface using hydrostatic balance, plus Coriolis terms, +!> using a discretisation using absolute temperature and constant +!> gravity g (shallow definition). Then integrate from the surface to +!> the top but using the model geopotential (may be defined as shallow +!> or deep depending on configuration). +module regrav_isotherm_kernel_mod + +use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_SCALAR, GH_INTEGER, & + GH_READ, GH_READWRITE, & + CELL_COLUMN +use constants_mod, only : r_def, i_def +use fs_continuity_mod, only : Wtheta, W3 +use kernel_mod, only : kernel_type implicit none @@ -25,67 +30,70 @@ module regrav_geopot_kernel_mod !------------------------------------------------------------------------------- ! Public types !------------------------------------------------------------------------------- -type, public, extends(kernel_type) :: regrav_geopot_kernel_type +type, public, extends(kernel_type) :: regrav_isotherm_kernel_type private - type(arg_type) :: meta_args(9) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & /) integer :: operates_on = CELL_COLUMN contains - procedure, nopass :: regrav_geopot_code + procedure, nopass :: regrav_isotherm_code end type !------------------------------------------------------------------------------- ! Contained functions/subroutines !------------------------------------------------------------------------------- -public :: regrav_geopot_code +public :: regrav_isotherm_code contains -!> @param[in] nlayers Number of layers -!> @param[in,out] temperature Absolute temperature field -!> @param[in,out] theta Potential temperature field -!> @param[in,out] exner Exner pressure field -!> @param[in] coriolis_term Vertical component of the coriolis term -!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!> @param[in] height_w3 Height coordinate in w3 -!> @param[in] height_wth Height coordinate in wth -!> @param[in] w3_mask LBC mask or Dummy mask for w3 space -!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!> @param[in] undf_wt Total number of degrees of freedom for wtheta -!> @param[in] map_wt Dofmap for the cell at column base for wt -!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!> @param[in] undf_w3 Total number of degrees of freedom for w3 -!> @param[in] map_w3 Dofmap for the cell at column base for w3 -subroutine regrav_geopot_code( nlayers, & - temperature, & - theta, & - exner, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - height_w3, & - height_wth, & - w3_mask, & - ndf_wt, & - undf_wt, & - map_wt, & - ndf_w3, & - undf_w3, & - map_w3 ) +!> @param[in] nlayers Number of layers +!> @param[in,out] exner Exner pressure field +!> @param[in] temperature Absolute temperature field +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] phi Geopotential field +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +subroutine regrav_isotherm_code( nlayers, & + exner, & + temperature, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + phi, & + height_w3, & + height_wth, & + w3_mask, & + ndf_w3, & + undf_w3, & + map_w3, & + ndf_wt, & + undf_wt, & + map_wt ) + + use planet_config_mod, only: gravity, cp implicit none ! Arguments + integer(kind=i_def), intent(in) :: nlayers, & ndf_w3, & undf_w3, & @@ -94,24 +102,20 @@ subroutine regrav_geopot_code( nlayers, & integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - - - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot - real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask, & + phi + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot + real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & + height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term ! Internal variables integer(kind=i_def) :: k - real(kind=r_def) :: temp_virt - real(kind=r_def) :: ht_wt(0:nlayers) - real(kind=r_def) :: th(0:nlayers) - real(kind=r_def) :: exner_surf - real(kind=r_def) :: weight1 + real(kind=r_def) :: temp_virtual + real(kind=r_def) :: dz, weight1, weight2 ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) ! setting exner to 1 @@ -122,87 +126,26 @@ subroutine regrav_geopot_code( nlayers, & return end if - ! Geopotential height - do k = 0, nlayers - ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & - ( planet_radius - height_wth(map_wt(1)+k) ) - end do + ! Integrate from surface to top using model geopotential phi + do k = 1, nlayers-1 - temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & - moist_dyn_tot( map_wt(1) ) - weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) - exner_surf = exner( map_w3(1) ) * ( cp * temp_virt ) / & - ( cp * temp_virt - ( gravity - coriolis_term( map_wt(1) ) ) * weight1 ) - theta( map_wt(1) ) = temperature( map_wt(1) ) / exner_surf - - -! Reuse th as temporary storage for T - do k = 0, nlayers - th(k) = temperature( map_wt(1) + k ) - end do + temp_virtual = moist_dyn_gas( map_wt(1) + k ) * temperature( map_wt(1) + k ) / & + moist_dyn_tot( map_wt(1) + k ) - ! Map temperature from newly computed grid back onto current model grid - call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & - temperature(map_wt(1)) ) + dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) + weight1 = ( height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) ) / dz + weight2 = ( height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) ) /dz -end subroutine regrav_geopot_code + ! Pi_k = Pi_k-1 * ( cp T_k-1/2 - ( (phi_k - phi_k-1) - F_k-1/2 * dz) * (1-w) ) / + ! ( cp T_k-1/2 + ( (phi_k - phi_k-1) - F_k-1/2 * dz) * w ) + exner( map_w3(1) + k ) = exner( map_w3 (1) + k-1 ) * & + ( cp * temp_virtual - ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & + - coriolis_term( map_wt(1) + k ) * dz ) * weight1 ) / & + ( cp * temp_virtual + ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & + - coriolis_term( map_wt(1) + k ) * dz ) * weight2 ) -subroutine interp( nin, zin, fin, nout, zout, fout ) - -implicit none - -integer(kind=i_def), intent(in) :: nin, nout -real(kind=r_def), intent(in) :: zin(nin), fin(nin), zout(nout) -real(kind=r_def), intent(inout) :: fout(nout) - -integer(kind=i_def) :: kout, kin -real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max -real(kind=r_def) :: eta_out, eta_1, eta_2 - -if ( zin(1) > zin(nin) ) then - zmax = zin(1) - zmin = zin(nin) - f_at_max = fin(1) - f_at_min = fin(nin) -else - zmax = zin(nin) - zmin = zin(1) - f_at_max = fin(nin) - f_at_min = fin(1) -end if - - -do kout = 1, nout - if ( zout(kout) <= zmin ) then - fout(kout) = f_at_min - else if ( zout(kout) >= zmax ) then - fout(kout) = f_at_max - else - do kin = 1, nin - 1 - if ( ( zout(kout) - zin(kin) ) * & - ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then - eta_out = sqrt(max(zout(kout), 0.0_r_def)) - eta_1 = sqrt(max(zin(kin), 0.0_r_def)) - eta_2 = sqrt(max(zin(kin+1), 0.0_r_def)) - eta_out = zout(kout) - eta_1 = zin(kin) - eta_2 = zin(kin+1) - xi = ( eta_out - eta_1 ) / ( eta_2 - eta_1 ) - fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) - fout(kout) = exp(fout(kout)) - end if - end do - end if -end do - -end subroutine interp + end do -! Does the input UM data, T(z), pi(z), have the correct z information? -! If T and pi have been interpolated onto the LFRic grid by UM2LFRic, why -! is there a problem with the lower boundary? I'm assuming this problem is -! that the orographic height of UM and LFRic are different. but if they are, -! then the UM and LFRic grids are different and the T(z), pi(z) are not -! valid at LFRic grid heights. Need to look at precisely what data is -! coming from um2lfric. +end subroutine regrav_isotherm_code -end module regrav_geopot_kernel_mod +end module regrav_isotherm_kernel_mod From cd41d2b740c6e711273c7f76743c2f764b4241a1 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Wed, 13 May 2026 08:47:33 +0100 Subject: [PATCH 16/27] Unit test for pressure diag --- .../source/algorithm/pres_lev_diags_alg_mod.x90 | 4 ++-- .../unit-test/kernel/geo_on_pres_kernel_mod_test.pf | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index 4f4c1fe92..9c5879692 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -21,8 +21,8 @@ module pres_lev_diags_alg_mod use fs_continuity_mod, only: W3, WTheta use mr_indices_mod, only: nummr, imr_v use moist_dyn_mod, only: num_moist_factors, total_mass - use planet_config_mod, only: p_zero, kappa, gravity, cp, & - planet_radius + use planet_config_mod, only: p_zero, kappa, gravity, cp + use extrusion_config_mod, only: planet_radius use planet_constants_mod, only: ex_power implicit none diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf index 69d9d0086..52f594101 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf @@ -89,6 +89,7 @@ contains real(r_def), parameter :: ex_power = 1.0_r_def real(r_def), parameter :: cp = 1000.0_r_def real(r_def), parameter :: gravity = 10.0_r_def + real(r_def), parameter :: planet_radius = 6.4e6_r_def ! Input arrays this%map_in = 1_i_def @@ -129,6 +130,7 @@ contains kappa, & cp, & gravity, & + planet_radius, & ex_power, & ndf_in, & undf_in, & From c7a82275c8cc8b9dabe5af45eb9e856ba9fa2456 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 14 May 2026 09:31:14 +0100 Subject: [PATCH 17/27] Corrections to geopotential height diagnostic --- .../algorithm/pres_lev_diags_alg_mod.x90 | 4 +- .../source/kernel/geo_on_pres_kernel_mod.F90 | 60 +++++++++++-------- .../source/psy/psykal_lite_phys_mod.F90 | 47 ++++++++------- .../kernel/geo_on_pres_kernel_mod_test.pf | 20 ++++--- 4 files changed, 76 insertions(+), 55 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index 9c5879692..0fb3c8b46 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -23,6 +23,7 @@ module pres_lev_diags_alg_mod use moist_dyn_mod, only: num_moist_factors, total_mass use planet_config_mod, only: p_zero, kappa, gravity, cp use extrusion_config_mod, only: planet_radius + use formulation_config_mod, only: shallow use planet_constants_mod, only: ex_power implicit none @@ -268,7 +269,8 @@ contains height_wth, exner_wth, & nplev, plevs, plev_geopot, & p_zero, kappa, cp, gravity, & - planet_radius, ex_power) + planet_radius, ex_power, & + shallow) if (plev_geopot_flag) call plev_geopot%write_field() if (plev_geopot_clim_flag) then diff --git a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 index c045fa94f..8aa2ae93a 100644 --- a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 @@ -7,16 +7,16 @@ module geo_on_pres_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 fs_continuity_mod, only: WTHETA - use constants_mod, only: r_def, i_def - use kernel_mod, only: kernel_type + use argument_mod, only: arg_type, & + GH_FIELD, GH_SCALAR, & + GH_READ, GH_WRITE, & + GH_INTEGER, GH_LOGICAL, & + GH_REAL, CELL_COLUMN, & + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + use fs_continuity_mod, only: WTHETA + use constants_mod, only: r_def, i_def, l_def + use kernel_mod, only: kernel_type implicit none @@ -25,20 +25,21 @@ module geo_on_pres_kernel_mod !> Kernel metadata for PSyclone type, public, extends(kernel_type) :: geo_on_pres_kernel_type private - type(arg_type) :: meta_args(12) = (/ & - 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_FIELD, GH_REAL, GH_READ, WTHETA), & - arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & - arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & - arg_type(GH_SCALAR,GH_INTEGER, GH_READ), & + type(arg_type) :: meta_args(13) = (/ & + 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_FIELD, GH_REAL, GH_READ, WTHETA), & + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & + 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_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ) & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_LOGICAL, GH_READ) & /) integer :: operates_on = CELL_COLUMN contains @@ -67,6 +68,7 @@ module geo_on_pres_kernel_mod !> @param[in] gravity Acceleration due to gravity !> @param[in] planet_radius Radius of the planet !> @param[in] ex_power (cp * lapse_rate ) / g + !> @param[in] shallow If true, gravity is constant !> @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 @@ -88,6 +90,7 @@ subroutine geo_on_pres_code(nlayers, & gravity, & planet_radius, & ex_power, & + shallow, & ndf_in, undf_in, map_in, & ndf_wth, undf_wth, map_wth, & ndf_out, undf_out, map_out) @@ -115,6 +118,7 @@ subroutine geo_on_pres_code(nlayers, & real(kind=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power, & planet_radius real(kind=r_def), intent(in), dimension(nplev) :: plevs + logical(kind=l_def), intent(in) :: shallow ! Internal variables integer(kind=i_def) :: k, level_above, top_df, kp, level_extrap @@ -154,9 +158,13 @@ subroutine geo_on_pres_code(nlayers, & exit end if end do - geo_ht_extrap = height_wth(map_wth(1)+level_extrap) / & - ( 1.0_r_def + height_wth(map_wth(1)+level_extrap) / & - planet_radius ) + if ( shallow ) then + geo_ht_extrap = height_wth(map_wth(1)+level_extrap) + else + geo_ht_extrap = height_wth(map_wth(1)+level_extrap) / & + ( 1.0_r_def + height_wth(map_wth(1)+level_extrap) / & + planet_radius ) + end if geo_ht_lowest = data_in(map_in(1)) / gravity h_ref_lev = ( geo_ht_extrap - geo_ht_lowest ) & / (1.0_r_def - (exner_wth(map_wth(1)+level_extrap) / & 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 2f19b600b..b66bae444 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 @@ -10,7 +10,7 @@ module psykal_lite_phys_mod - use constants_mod, only : i_def, r_def + use constants_mod, only : i_def, r_def, l_def use field_mod, only : field_type, field_proxy_type use integer_field_mod, only : integer_field_type, integer_field_proxy_type use mesh_mod, only : mesh_type @@ -991,15 +991,18 @@ END SUBROUTINE invoke_pres_interp_kernel_type !> see https://github.com/stfc/PSyclone/issues/1312 !> Hence this module could be removed once the PSyclone ticket is !> completed - SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height_wth, exner_wth, nplev, plevs, plev_geopot, & -&p_zero, kappa, cp, gravity, planet_radius, ex_power) + SUBROUTINE invoke_geo_on_pres_kernel_type(geopot_w3, exner_w3, & + theta_wth, height_wth, exner_wth, & + nplev, plevs, plev_geopot, & + p_zero, kappa, cp, gravity, planet_radius, ex_power, shallow) USE geo_on_pres_kernel_mod, ONLY: geo_on_pres_code USE mesh_mod, ONLY: mesh_type REAL(KIND=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power, & planet_radius INTEGER(KIND=i_def), intent(in) :: nplev + LOGICAL(KIND=l_def), intent(in) :: shallow REAL(KIND=r_def), intent(in) :: plevs(nplev) - TYPE(field_type), intent(in) :: height_w3, exner_w3, theta_wth, height_wth, exner_wth, plev_geopot + TYPE(field_type), intent(in) :: geopot_w3, exner_w3, theta_wth, height_wth, exner_wth, plev_geopot INTEGER(KIND=i_def) cell INTEGER(KIND=i_def) loop0_start, loop0_stop INTEGER(KIND=i_def) nlayers @@ -1008,19 +1011,19 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height REAL(KIND=r_def), pointer, dimension(:) :: height_wth_data => null() REAL(KIND=r_def), pointer, dimension(:) :: theta_wth_data => null() REAL(KIND=r_def), pointer, dimension(:) :: exner_w3_data => null() - REAL(KIND=r_def), pointer, dimension(:) :: height_w3_data => null() - TYPE(field_proxy_type) height_w3_proxy, exner_w3_proxy, theta_wth_proxy, height_wth_proxy, exner_wth_proxy, plev_geopot_proxy - INTEGER(KIND=i_def), pointer :: map_adspc1_height_w3(:,:) => null(), map_adspc2_plev_geopot(:,:) => null(), & + REAL(KIND=r_def), pointer, dimension(:) :: geopot_w3_data => null() + TYPE(field_proxy_type) geopot_w3_proxy, exner_w3_proxy, theta_wth_proxy, height_wth_proxy, exner_wth_proxy, plev_geopot_proxy + INTEGER(KIND=i_def), pointer :: map_adspc1_geopot_w3(:,:) => null(), map_adspc2_plev_geopot(:,:) => null(), & &map_wtheta(:,:) => null() - INTEGER(KIND=i_def) ndf_adspc1_height_w3, undf_adspc1_height_w3, ndf_wtheta, undf_wtheta, ndf_adspc2_plev_geopot, & + INTEGER(KIND=i_def) ndf_adspc1_geopot_w3, undf_adspc1_geopot_w3, ndf_wtheta, undf_wtheta, ndf_adspc2_plev_geopot, & &undf_adspc2_plev_geopot INTEGER(KIND=i_def) max_halo_depth_mesh TYPE(mesh_type), pointer :: mesh => null() ! ! Initialise field and/or operator proxies ! - height_w3_proxy = height_w3%get_proxy() - height_w3_data => height_w3_proxy%data + geopot_w3_proxy = geopot_w3%get_proxy() + geopot_w3_data => geopot_w3_proxy%data exner_w3_proxy = exner_w3%get_proxy() exner_w3_data => exner_w3_proxy%data theta_wth_proxy = theta_wth%get_proxy() @@ -1034,23 +1037,23 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height ! ! Initialise number of layers ! - nlayers = height_w3_proxy%vspace%get_nlayers() + nlayers = geopot_w3_proxy%vspace%get_nlayers() ! ! Create a mesh object ! - mesh => height_w3_proxy%vspace%get_mesh() + mesh => geopot_w3_proxy%vspace%get_mesh() max_halo_depth_mesh = mesh%get_halo_depth() ! ! Look-up dofmaps for each function space ! - map_adspc1_height_w3 => height_w3_proxy%vspace%get_whole_dofmap() + map_adspc1_geopot_w3 => geopot_w3_proxy%vspace%get_whole_dofmap() map_wtheta => theta_wth_proxy%vspace%get_whole_dofmap() map_adspc2_plev_geopot => plev_geopot_proxy%vspace%get_whole_dofmap() ! - ! Initialise number of DoFs for adspc1_height_w3 + ! Initialise number of DoFs for adspc1_geopot_w3 ! - ndf_adspc1_height_w3 = height_w3_proxy%vspace%get_ndf() - undf_adspc1_height_w3 = height_w3_proxy%vspace%get_undf() + ndf_adspc1_geopot_w3 = geopot_w3_proxy%vspace%get_ndf() + undf_adspc1_geopot_w3 = geopot_w3_proxy%vspace%get_undf() ! ! Initialise number of DoFs for wtheta ! @@ -1071,10 +1074,14 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height ! DO cell=loop0_start,loop0_stop ! - CALL geo_on_pres_code(nlayers, height_w3_data, exner_w3_data, theta_wth_data, height_wth_data, exner_wth_data, nplev, & -&plevs, plev_geopot_data, p_zero, kappa, cp, gravity, planet_radius, ex_power, ndf_adspc1_height_w3, undf_adspc1_height_w3, & -&map_adspc1_height_w3(:,cell), ndf_wtheta, undf_wtheta, map_wtheta(:,cell), ndf_adspc2_plev_geopot, undf_adspc2_plev_geopot, & -&map_adspc2_plev_geopot(:,cell)) + CALL geo_on_pres_code(nlayers, geopot_w3_data, exner_w3_data, & + theta_wth_data, height_wth_data, exner_wth_data, & + nplev, plevs, plev_geopot_data, & + p_zero, kappa, cp, gravity, planet_radius, ex_power, shallow, & + ndf_adspc1_geopot_w3, undf_adspc1_geopot_w3, map_adspc1_geopot_w3(:,cell), & + ndf_wtheta, undf_wtheta, map_wtheta(:,cell), & + ndf_adspc2_plev_geopot, undf_adspc2_plev_geopot, & + map_adspc2_plev_geopot(:,cell)) END DO ! ! Set halos dirty/clean for fields modified in the above loop diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf index 52f594101..18a184a10 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf @@ -8,7 +8,7 @@ module geo_on_pres_kernel_mod_test use funit - use constants_mod, only : r_def, i_def + use constants_mod, only : r_def, i_def, l_def implicit none @@ -83,13 +83,14 @@ contains class(geo_on_pres_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-14_r_def - real(r_def), parameter :: p_zero = 100000.0_r_def - real(r_def), parameter :: kappa = 1.0_r_def - real(r_def), parameter :: ex_power = 1.0_r_def - real(r_def), parameter :: cp = 1000.0_r_def - real(r_def), parameter :: gravity = 10.0_r_def - real(r_def), parameter :: planet_radius = 6.4e6_r_def + real(r_def), parameter :: tol = 1.0e-14_r_def + real(r_def), parameter :: p_zero = 100000.0_r_def + real(r_def), parameter :: kappa = 1.0_r_def + real(r_def), parameter :: ex_power = 1.0_r_def + real(r_def), parameter :: cp = 1000.0_r_def + real(r_def), parameter :: gravity = 10.0_r_def + real(r_def), parameter :: planet_radius = 6.4e6_r_def + logical(l_def), parameter :: shallow = .true. ! Input arrays this%map_in = 1_i_def @@ -106,6 +107,8 @@ contains this%data_in(1) = 0.0_r_def this%data_in(2) = 100.0_r_def this%data_in(3) = 1000.0_r_def + ! Conver heights to geopotential + this%data_in(:) = gravity * this%data_in(:) this%height(1) = 0.0_r_def this%height(2) = 100.0_r_def this%height(3) = 3000.0_r_def @@ -132,6 +135,7 @@ contains gravity, & planet_radius, & ex_power, & + shallow, & ndf_in, & undf_in, & this%map_in, & From ebbf4f3fac004eab1654c4303e09361dc5837d86 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 21 May 2026 16:26:47 +0100 Subject: [PATCH 18/27] Improve comments. Remove redundant code. --- .../physics/map_fd_to_prognostics_alg_mod.x90 | 39 ++++--- .../regrav_geopot_kernel_mod.F90 | 101 ++++++------------ 2 files changed, 56 insertions(+), 84 deletions(-) diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 8049e99e9..bf7ecf359 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -468,9 +468,10 @@ contains ! Compute absolute temperature T = theta * exner X_times_Y( temperature, exner_in_wth, theta), & - ! Rebalance exner in vertical, using hydrostatic balance. + ! Rebalance exner in vertical, using hydrostatic balance + ! preserving absolute temperature. ! This integrates down using g (shallow) from eos_level and - ! then integrates up using phi (deep). + ! then integrates up using phi (deep). hydro_shallow_to_deep_kernel_type( exner, temperature, & c_vert_term, & moist_dyn, geopotential, & @@ -489,15 +490,18 @@ contains ! Compute absolute temperature T = theta * exner X_times_Y( temperature, exner_in_wth, theta), & - ! Rebalance exner in vertical, using hydrostatic balance. - ! This integrates down using g (shallow) from eos_level and - ! then integrates up using phi (deep). + ! Find exner at level 1 by integrating hydrostatic equation + ! downwards from eos_level, preserving theta. Gravity + ! assumed constant. hstat_cori_g_const_kernel_type( exner, rho, theta, & c_vert_term, & moist_dyn(gas_law), & moist_dyn(total_mass), & height_w3, mask, & eos_level ), & + ! Find exner at all levels by integrating hydrostatic + ! equation upwards from level 1, preserving absolute + ! temperature. Gravity is height-dependent. regrav_isotherm_kernel_type( exner, temperature, & c_vert_term, & moist_dyn(gas_law), & @@ -507,17 +511,23 @@ contains mask ) ) case( regravitate_isotherm3 ) call invoke( & + ! Find exner below eos_level by integrating hydrostatic + ! equation downwards, assuming constant gravity. hstat_cori_g_const_kernel_type( exner, rho, theta, & c_vert_term, & moist_dyn(gas_law), & moist_dyn(total_mass), & height_w3, mask, & eos_level ), & + ! Find exner in Wtheta space sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & ! Compute absolute temperature T = theta * exner X_times_Y( temperature, exner_in_wth, theta), & + ! Find exner at all levels by integrating hydrostatic + ! equation upwards from level 1, preserving absolute + ! temperature. Gravity is height-dependent. regrav_isotherm_kernel_type( exner, temperature, & c_vert_term, & moist_dyn(gas_law), & @@ -527,23 +537,26 @@ contains mask ) ) case( regravitate_geopot ) call invoke( & + ! Find exner below eos_level by integrating hydrostatic + ! equation downwards, assuming constant gravity. hstat_cori_g_const_kernel_type( exner, rho, theta, & c_vert_term, & moist_dyn(gas_law), & moist_dyn(total_mass), & height_w3, mask, & eos_level ), & + ! Find exner in Wtheta space sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - regrav_geopot_kernel_type( temperature, theta, & - exner, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - height_w3, height_wth, & - mask ), & + X_times_Y( temperature, exner_in_wth, theta), & + ! Interpolate temperature from heights corresponding to + ! geopotential height onto heights of the model grid. + regrav_geopot_kernel_type( temperature, exner, & + height_wth, mask ), & + ! Find exner at all levels by integrating hydrostatic + ! equation upwards from level 1, preserving absolute + ! temperature. Gravity is height-dependent. regrav_isotherm_kernel_type( exner, temperature, & c_vert_term, & moist_dyn(gas_law), & diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 index 00ed128fe..bb0504f04 100644 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -4,8 +4,15 @@ ! under which the code may be used. !----------------------------------------------------------------------------- -!> @brief Need to say something here -!> @details More to be said here +!> @brief Convert temperature on geopotential height to temperature on height +!> @details The input temperature is assumed to have been generated by a model +!> with constant gravity. For such a model physical height and +!> geopotential height are identical. In this kernel gravity is assumed +!> to vary correctly with height, so that physical and geopotential +!> heights are no longer identical. The physical height, height_phys, +!> corresponding to height_wth, interpreted as geopotential height, is +!> calculated. Temperature is then interpolated from height_phys onto +!> model levels height_wth. module regrav_geopot_kernel_mod use argument_mod, only: arg_type, & @@ -16,7 +23,6 @@ module regrav_geopot_kernel_mod use extrusion_config_mod, only: planet_radius use fs_continuity_mod, only: Wtheta, W3 use kernel_mod, only: kernel_type -use planet_config_mod, only: gravity, cp implicit none @@ -27,13 +33,8 @@ module regrav_geopot_kernel_mod !------------------------------------------------------------------------------- type, public, extends(kernel_type) :: regrav_geopot_kernel_type private - type(arg_type) :: meta_args(9) = (/ & + type(arg_type) :: meta_args(4) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & @@ -52,12 +53,7 @@ module regrav_geopot_kernel_mod !> @param[in] nlayers Number of layers !> @param[in,out] temperature Absolute temperature field -!> @param[in,out] theta Potential temperature field -!> @param[in,out] exner Exner pressure field -!> @param[in] coriolis_term Vertical component of the coriolis term -!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] exner Exner pressure field !> @param[in] height_wth Height coordinate in wth !> @param[in] w3_mask LBC mask or Dummy mask for w3 space !> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta @@ -68,12 +64,7 @@ module regrav_geopot_kernel_mod !> @param[in] map_w3 Dofmap for the cell at column base for w3 subroutine regrav_geopot_code( nlayers, & temperature, & - theta, & exner, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - height_w3, & height_wth, & w3_mask, & ndf_wt, & @@ -86,32 +77,23 @@ subroutine regrav_geopot_code( nlayers, & implicit none ! Arguments - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - real(kind=r_def), dimension(undf_wt), intent(inout) :: theta, temperature - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - - - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot - real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), dimension(undf_wt), intent(inout) :: temperature + real(kind=r_def), dimension(undf_w3), intent(in) :: exner + real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth + real(kind=r_def), dimension(undf_w3), intent(in) :: w3_mask ! Internal variables integer(kind=i_def) :: k - real(kind=r_def) :: temp_virt - real(kind=r_def) :: ht_wt(0:nlayers) - real(kind=r_def) :: th(0:nlayers) - real(kind=r_def) :: exner_surf - real(kind=r_def) :: weight1 + real(kind=r_def) :: height_phys(0:nlayers) + real(kind=r_def) :: temp(0:nlayers) ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) ! setting exner to 1 @@ -122,27 +104,19 @@ subroutine regrav_geopot_code( nlayers, & return end if - ! Geopotential height + ! Height corresponding to height_wth when latter represents geopotential height do k = 0, nlayers - ht_wt(k) = planet_radius * height_wth(map_wt(1)+k) / & - ( planet_radius - height_wth(map_wt(1)+k) ) + height_phys(k) = planet_radius * height_wth(map_wt(1)+k) / & + ( planet_radius - height_wth(map_wt(1)+k) ) end do - temp_virt = moist_dyn_gas( map_wt(1) ) * temperature( map_wt(1) ) / & - moist_dyn_tot( map_wt(1) ) - weight1 = height_w3( map_w3(1) ) - height_wth( map_wt(1) ) - exner_surf = exner( map_w3(1) ) * ( cp * temp_virt ) / & - ( cp * temp_virt - ( gravity - coriolis_term( map_wt(1) ) ) * weight1 ) - theta( map_wt(1) ) = temperature( map_wt(1) ) / exner_surf - - -! Reuse th as temporary storage for T +! Temporary storage for T do k = 0, nlayers - th(k) = temperature( map_wt(1) + k ) + temp(k) = temperature( map_wt(1) + k ) end do - ! Map temperature from newly computed grid back onto current model grid - call interp( nlayers+1, ht_wt, th, nlayers+1, height_wth(map_wt(1)), & + ! Interpolate temperature from newly computed grid back onto current model grid + call interp( nlayers+1, height_phys, temp, nlayers+1, height_wth(map_wt(1)), & temperature(map_wt(1)) ) end subroutine regrav_geopot_code @@ -157,7 +131,6 @@ subroutine interp( nin, zin, fin, nout, zout, fout ) integer(kind=i_def) :: kout, kin real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max -real(kind=r_def) :: eta_out, eta_1, eta_2 if ( zin(1) > zin(nin) ) then zmax = zin(1) @@ -181,13 +154,7 @@ subroutine interp( nin, zin, fin, nout, zout, fout ) do kin = 1, nin - 1 if ( ( zout(kout) - zin(kin) ) * & ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then - eta_out = sqrt(max(zout(kout), 0.0_r_def)) - eta_1 = sqrt(max(zin(kin), 0.0_r_def)) - eta_2 = sqrt(max(zin(kin+1), 0.0_r_def)) - eta_out = zout(kout) - eta_1 = zin(kin) - eta_2 = zin(kin+1) - xi = ( eta_out - eta_1 ) / ( eta_2 - eta_1 ) + xi = ( zout(kout) - zin(kin) ) / ( zin(kin+1) - zin(kin) ) fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) fout(kout) = exp(fout(kout)) end if @@ -197,12 +164,4 @@ subroutine interp( nin, zin, fin, nout, zout, fout ) end subroutine interp -! Does the input UM data, T(z), pi(z), have the correct z information? -! If T and pi have been interpolated onto the LFRic grid by UM2LFRic, why -! is there a problem with the lower boundary? I'm assuming this problem is -! that the orographic height of UM and LFRic are different. but if they are, -! then the UM and LFRic grids are different and the T(z), pi(z) are not -! valid at LFRic grid heights. Need to look at precisely what data is -! coming from um2lfric. - end module regrav_geopot_kernel_mod From 0e3512d0b16f3a84ea44a41d4edd4f0b32a735b1 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 21 May 2026 16:59:39 +0100 Subject: [PATCH 19/27] Exner is inout --- .../kernel/initialisation/regrav_geopot_kernel_mod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 index bb0504f04..ea7d1c8f4 100644 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -35,7 +35,7 @@ module regrav_geopot_kernel_mod private type(arg_type) :: meta_args(4) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & /) @@ -53,7 +53,7 @@ module regrav_geopot_kernel_mod !> @param[in] nlayers Number of layers !> @param[in,out] temperature Absolute temperature field -!> @param[in] exner Exner pressure field +!> @param[in,out] exner Exner pressure field !> @param[in] height_wth Height coordinate in wth !> @param[in] w3_mask LBC mask or Dummy mask for w3 space !> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta @@ -86,7 +86,7 @@ subroutine regrav_geopot_code( nlayers, & integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt real(kind=r_def), dimension(undf_wt), intent(inout) :: temperature - real(kind=r_def), dimension(undf_w3), intent(in) :: exner + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth real(kind=r_def), dimension(undf_w3), intent(in) :: w3_mask From 462e1bb02c14f0c444a046f580222bd537c1a2e0 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Fri, 22 May 2026 15:39:36 +0100 Subject: [PATCH 20/27] Integration must be from top --- .../algorithm/physics/map_fd_to_prognostics_alg_mod.x90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index bf7ecf359..a11674f03 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -387,7 +387,7 @@ contains type(mesh_type), pointer :: mesh type(field_type), pointer :: mask type(field_type), target :: dummy_mask - integer(i_def) :: eos_level + integer(i_def) :: eos_level, nlayers nullify(geopotential, chi, panel_id, height_w3, & height_wth, vert_coriolis, mesh, mask ) @@ -431,7 +431,8 @@ contains ! Get the index of the level closest to the height, ignoring the top eta ! level as eos_level is to be used for a W3 space rather than Wtheta eos_level = get_nearest_level( mesh, eos_height ) - if ( eos_level == (mesh%get_nlayers() + 1) ) eos_level = eos_level - 1 + nlayers = mesh%get_nlayers() + if ( eos_level == nlayers + 1 ) eos_level = nlayers if (shallow) then ! Initialize the vertical profile by starting with the equation of state at @@ -518,7 +519,7 @@ contains moist_dyn(gas_law), & moist_dyn(total_mass), & height_w3, mask, & - eos_level ), & + nlayers ), & ! Find exner in Wtheta space sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & From c4e3ec1f3abf25b54a2be3a8972b750a43365506 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Fri, 22 May 2026 16:08:37 +0100 Subject: [PATCH 21/27] Remove obsolete code --- .../physics/map_fd_to_prognostics_alg_mod.x90 | 25 --- .../hydro_shallow_to_deep_kernel_mod.F90 | 211 ------------------ 2 files changed, 236 deletions(-) delete mode 100644 science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index a11674f03..7d94d809e 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -349,8 +349,6 @@ contains only : hydrostatic_coriolis_kernel_type use hstat_cori_g_const_kernel_mod, & only : hstat_cori_g_const_kernel_type - use hydro_shallow_to_deep_kernel_mod, & - only : hydro_shallow_to_deep_kernel_type use regrav_isotherm_kernel_mod, only : regrav_isotherm_kernel_type use regrav_geopot_kernel_mod, only : regrav_geopot_kernel_type use moist_dyn_mod, only : num_moist_factors @@ -456,29 +454,6 @@ contains call theta%copy_field_properties(temperature) select case( regravitate ) - case( regravitate_isotherm1 ) - call invoke( & - ! Compute exner from EoS - sample_eos_pressure_kernel_type( exner, rho, & - theta, moist_dyn(gas_law) ), & - - ! Compute exner in Wtheta (required to calculate T) - sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & - height_wth, height_w3 ), & - - ! Compute absolute temperature T = theta * exner - X_times_Y( temperature, exner_in_wth, theta), & - - ! Rebalance exner in vertical, using hydrostatic balance - ! preserving absolute temperature. - ! This integrates down using g (shallow) from eos_level and - ! then integrates up using phi (deep). - hydro_shallow_to_deep_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn, geopotential, & - height_w3, height_wth, & - mask, gravity, & - cp, eos_level ) ) case( regravitate_isotherm2 ) call invoke( & ! Compute exner from EoS diff --git a/science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 deleted file mode 100644 index a8b733462..000000000 --- a/science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 +++ /dev/null @@ -1,211 +0,0 @@ -!----------------------------------------------------------------------------- -! (C) Crown copyright 2025 Met Office. All rights reserved. -! The file LICENCE, distributed with this code, contains details of the terms -! under which the code may be used. -!----------------------------------------------------------------------------- - -!> @brief Defines exner pressure from vertical balance with a deep gravity, -!! but using input data based on a shallow gravity. -!> @details Using the input exner pressure at level eos_index, integrate down -!! to the surface using hydrostatic balance, plus Coriolis terms, -!! using a discretisation using absolute temperature and constant -!! gravity g (shallow definition). Then integrate from the surface to -!! the top but using the model geopotential (may be defined as shallow -!! or deep depending on configuration). -module hydro_shallow_to_deep_kernel_mod - -use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & - GH_SCALAR, GH_INTEGER, & - GH_READ, GH_READWRITE, & - ANY_SPACE_9, ANY_SPACE_1, & - GH_BASIS, CELL_COLUMN, GH_EVALUATOR -use constants_mod, only : r_def, i_def -use fs_continuity_mod, only : Wtheta, W3 -use kernel_mod, only : kernel_type -use reference_element_mod, only : B - -implicit none - -private - -!------------------------------------------------------------------------------- -! Public types -!------------------------------------------------------------------------------- -type, public, extends(kernel_type) :: hydro_shallow_to_deep_kernel_type - private - type(arg_type) :: meta_args(11) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & - /) - type(func_type) :: meta_funcs(2) = (/ & - func_type(W3, GH_BASIS), & - func_type(Wtheta, GH_BASIS) & - /) - integer :: operates_on = CELL_COLUMN - integer :: gh_shape = GH_EVALUATOR -contains - procedure, nopass :: hydro_shallow_to_deep_code -end type - -!------------------------------------------------------------------------------- -! Contained functions/subroutines -!------------------------------------------------------------------------------- -public :: hydro_shallow_to_deep_code - -contains - -!! @param[in] nlayers Number of layers -!! @param[in,out] exner Exner pressure field -!! @param[in] temperature Absolute temperature field -!! @param[in] coriolis_term Vertical component of the coriolis term -!! @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!! @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!! @param[in] moist_dyn_fac Water factor -!! @param[in] phi Geopotential field -!! @param[in] height_w3 Height coordinate in w3 -!! @param[in] height_wth Height coordinate in wth -!! @param[in] w3_mask LBC mask or Dummy mask for w3 space -!! @param[in] gravity The planet gravity -!! @param[in] cp Specific heat of dry air at constant pressure -!! @param[in] eos_index Vertical level at which the equation of state -!! is satisfied -!! @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!! @param[in] undf_w3 Total number of degrees of freedom for w3 -!! @param[in] map_w3 Dofmap for the cell at column base for w3 -!! @param[in] basis_w3 Basis functions evaluated at w3 nodes -!! @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!! @param[in] undf_wt Total number of degrees of freedom for wtheta -!! @param[in] map_wt Dofmap for the cell at column base for wt -!! @param[in] basis_wt Basis functions evaluated at wt nodes -subroutine hydro_shallow_to_deep_code( & - nlayers, & - exner, & - temperature, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac, & - phi, & - height_w3, & - height_wth, & - w3_mask, & - gravity, & - cp, & - eos_index, & - ndf_w3, & - undf_w3, & - map_w3, & - basis_w3, & - ndf_wt, & - undf_wt, & - map_wt, & - basis_wt ) - - implicit none - - ! Arguments - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - integer(kind=i_def), intent(in) :: eos_index - - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask, & - phi - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac - real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & - height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term - real(kind=r_def), dimension(1,ndf_w3,ndf_w3), intent(in) :: basis_w3 - real(kind=r_def), dimension(1,ndf_wt,ndf_w3), intent(in) :: basis_wt - real(kind=r_def), intent(in) :: gravity - real(kind=r_def), intent(in) :: cp - - ! Internal variables - integer(kind=i_def) :: k, eos_index_m1 - real(kind=r_def) :: temp_moist, dz, & - weight1, weight2 - - ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) - ! setting exner to 1 - if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then - do k = 0, nlayers-1 - exner(map_w3(1) + k ) = 1.0_r_def - enddo - return - end if - - ! Levels in the kernel are numbered 0 to nlayers-1, rather than 1 to nlayers - eos_index_m1 = eos_index - 1 - - ! Integrate from eos_level to surface using shallow gravity - do k = eos_index_m1, 1, -1 - - ! Absolute temperature at which dry air would have same pressure and density - ! T_moist = T (1 + mr/eps)/(1 + sum mr) - temp_moist = moist_dyn_gas( map_wt(1) + k ) * & - temperature( map_wt(1) + k ) / & - moist_dyn_tot( map_wt(1) + k ) - - ! Discretisation weights e.g. weight 1 = w * dz where w = (z_k - z_k-1/2 ) / dz - ! and weight 2 = (1 - w) * dz, dz = z_k - z_k-1 - ! - ! z_k W3 k - ! z_k-1/2 Wtheta k - ! z_k-1 W3 k-1 - weight1 = height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) - weight2 = height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) - - ! Vertical momentum equation in hydrostatic balance (plus coriolis) with T = theta Pi - ! cp T d Pi /dz = Pi ( F - d Phi/ dz ) - ! Discretised, using d Phi /dz = g as - ! cp T_{k-1/2} ( Pi_k - Pi_k-1 ) / dz = ( w Pi_k + (1-w) Pi_k-1 ) ( F_k-1/2 - g) - ! Rearranged as - ! Pi_k-1 = Pi_k * ( cp T_{k-1/2} + (g - F_k-1/2) * (1-w) * dz ) / - ! ( cp T_{k-1/2} - (g - F_k-1/2) * w * dz ) - exner( map_w3 (1) + k-1 ) = exner( map_w3(1) + k ) * & - ( cp * temp_moist + ( gravity - coriolis_term( map_wt(1) + k ) ) * weight2 ) / & - ( cp * temp_moist - ( gravity - coriolis_term( map_wt(1) + k ) ) * weight1 ) - - end do - - ! Integrate from surface to top using model geopotential phi - do k = 1, nlayers-1 - - temp_moist = moist_dyn_gas( map_wt(1) + k ) * temperature( map_wt(1) + k ) / & - moist_dyn_tot( map_wt(1) + k ) - - dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) - weight1 = ( height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) ) / dz - weight2 = ( height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) ) /dz - - ! Pi_k = Pi_k-1 * ( cp T_k-1/2 - ( (phi_k - phi_k-1) - F_k-1/2 * dz) * (1-w) ) / - ! ( cp T_k-1/2 + ( (phi_k - phi_k-1) - F_k-1/2 * dz) * w ) - exner( map_w3(1) + k ) = exner( map_w3 (1) + k-1 ) * & - ( cp * temp_moist - ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & - - coriolis_term( map_wt(1) + k ) * dz ) * weight1 ) / & - ( cp * temp_moist + ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & - - coriolis_term( map_wt(1) + k ) * dz ) * weight2 ) - - end do - -end subroutine hydro_shallow_to_deep_code - -end module hydro_shallow_to_deep_kernel_mod From a263bec8fd3117a1048568ca5d0fccf0be33b3dd Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Fri, 22 May 2026 16:12:36 +0100 Subject: [PATCH 22/27] Contribute --- CONTRIBUTORS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 00c6e6edb..383f9863f 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -34,4 +34,5 @@ | thomasmelvin | Thomas Melvin | Met Office | 2026-01-15 | | tinyendian | Wolfgang Hayek | Earth Sciences New Zealand | 2026-02-02 | | DanStoneMO | Daniel Stone | Met Office | 2026-02-26 | -| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | \ No newline at end of file +| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | +| mo-cjsmith | Chris Smith | Met Office | 2026-05-22 | From 0a44f6c913507b4c90cbf96dfa49a9d690a768d7 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Fri, 22 May 2026 16:44:14 +0100 Subject: [PATCH 23/27] Improve metadata --- science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf index 066513127..208b7fdfa 100644 --- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf +++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf @@ -3229,7 +3229,7 @@ help=Method for switching gravity from constant to height-varying. =Used when the shallow switch is false and reading a start dump =that has been created by a constant gravity model, such as the UM. =All methods calculate an Exner field in hydrostatic balance with - =height-varying gravity. This is done by upward integr + =height-varying gravity. =Isothermodynamic: Preserve absolute temperature and pressure on = model levels by allowing the heights of levels = to change. Absolute temperature is then interpolated From 2cbcc58ef2a7a1fa91b2f4cd50a87af7ae8bd2df Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Thu, 28 May 2026 11:00:07 +0100 Subject: [PATCH 24/27] Use pre-existing kernel --- .../physics/map_fd_to_prognostics_alg_mod.x90 | 45 ++-- .../hstat_cori_g_const_kernel_mod.F90 | 194 ------------------ 2 files changed, 25 insertions(+), 214 deletions(-) delete mode 100644 science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 7d94d809e..0dbf9af6a 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -347,8 +347,6 @@ contains use sci_enforce_bc_kernel_mod, only : enforce_bc_kernel_type use hydrostatic_coriolis_kernel_mod, & only : hydrostatic_coriolis_kernel_type - use hstat_cori_g_const_kernel_mod, & - only : hstat_cori_g_const_kernel_type use regrav_isotherm_kernel_mod, only : regrav_isotherm_kernel_type use regrav_geopot_kernel_mod, only : regrav_geopot_kernel_type use moist_dyn_mod, only : num_moist_factors @@ -376,6 +374,7 @@ contains type(field_type ) :: c_vert_term, r_c type(field_type ) :: exner_in_wth type(field_type ) :: temperature + type(field_type ) :: geopotential_g_const type(field_type), pointer :: geopotential type(field_type), pointer :: chi(:) type(field_type), pointer :: panel_id @@ -452,6 +451,8 @@ contains ! by rederiving potential temperature (theta) and density (rho). call theta%copy_field_properties(exner_in_wth) call theta%copy_field_properties(temperature) + call geopotential%copy_field_properties(geopotential_g_const) + call invoke( a_times_X( geopotential_g_const, gravity, height_w3 ) ) select case( regravitate ) case( regravitate_isotherm2 ) @@ -469,12 +470,13 @@ contains ! Find exner at level 1 by integrating hydrostatic equation ! downwards from eos_level, preserving theta. Gravity ! assumed constant. - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - height_w3, mask, & - eos_level ), & + hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential_g_const, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & ! Find exner at all levels by integrating hydrostatic ! equation upwards from level 1, preserving absolute ! temperature. Gravity is height-dependent. @@ -489,12 +491,14 @@ contains call invoke( & ! Find exner below eos_level by integrating hydrostatic ! equation downwards, assuming constant gravity. - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - height_w3, mask, & - nlayers ), & + hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential_g_const, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + ! Find exner in Wtheta space sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & @@ -515,12 +519,13 @@ contains call invoke( & ! Find exner below eos_level by integrating hydrostatic ! equation downwards, assuming constant gravity. - hstat_cori_g_const_kernel_type( exner, rho, theta, & - c_vert_term, & - moist_dyn(gas_law), & - moist_dyn(total_mass), & - height_w3, mask, & - eos_level ), & + hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential_g_const, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & ! Find exner in Wtheta space sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & diff --git a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 deleted file mode 100644 index da7d29f19..000000000 --- a/science/gungho/source/kernel/initialisation/hstat_cori_g_const_kernel_mod.F90 +++ /dev/null @@ -1,194 +0,0 @@ -!----------------------------------------------------------------------------- -! (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. -!----------------------------------------------------------------------------- - -module hstat_cori_g_const_kernel_mod - -!> @brief Computes exner from equation of state and vertical balance including -!! Coriolis term. -!> @details Calculate exner on the top level, using equation of state. Then -!! use finite differencing to calculate exner on successive levels -!! below. Unlike hydrostatic_eos_kernel which balances upwards, the -!! Coriolis term is also added to the vertical balance equation. - -use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & - GH_SCALAR, GH_INTEGER, & - GH_READ, GH_WRITE, & - ANY_SPACE_9, ANY_SPACE_1, & - GH_BASIS, CELL_COLUMN, GH_EVALUATOR -use constants_mod, only : r_def, i_def -use fs_continuity_mod, only : Wtheta, W3 -use kernel_mod, only : kernel_type -use planet_config_mod, only : gravity, p_zero, kappa, rd, cp - -implicit none - -private - -!------------------------------------------------------------------------------- -! Public types -!------------------------------------------------------------------------------- -type, public, extends(kernel_type) :: hstat_cori_g_const_kernel_type - private - type(arg_type) :: meta_args(9) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & - /) - type(func_type) :: meta_funcs(2) = (/ & - func_type(W3, GH_BASIS), & - func_type(Wtheta, GH_BASIS) & - /) - integer :: operates_on = CELL_COLUMN - integer :: gh_shape = GH_EVALUATOR -contains - procedure, nopass :: hstat_cori_g_const_code -end type - -!------------------------------------------------------------------------------- -! Contained functions/subroutines -!------------------------------------------------------------------------------- -public :: hstat_cori_g_const_code - -contains - -!! @param[in] nlayers Number of layers -!! @param[in,out] exner Exner pressure field -!! @param[in] rho Density field -!! @param[in] theta Potential temperature field -!! @param[in] coriolis_term Vertical component of the coriolis term -!! @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!! @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!! @param[in] height_w3 Height coordinate in w3 -!! @param[in] w3_mask LBC mask or Dummy mask for w3 space -!! @param[in] eos_index Vertical level at which the equation of state -!! is satisfied -!! @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!! @param[in] undf_w3 Total number of degrees of freedom for w3 -!! @param[in] map_w3 Dofmap for the cell at column base for w3 -!! @param[in] basis_w3 Basis functions evaluated at w3 nodes -!! @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!! @param[in] undf_wt Total number of degrees of freedom for wtheta -!! @param[in] map_wt Dofmap for the cell at column base for wt -!! @param[in] basis_wt Basis functions evaluated at wt nodes -subroutine hstat_cori_g_const_code( nlayers, & - exner, & - rho, & - theta, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - height_w3, & - w3_mask, & - eos_index, & - ndf_w3, & - undf_w3, & - map_w3, & - basis_w3, & - ndf_wt, & - undf_wt, & - map_wt, & - basis_wt ) - - implicit none - - ! Arguments - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - integer(kind=i_def), intent(in) :: eos_index - - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - real(kind=r_def), dimension(undf_w3), intent(in) :: rho, & - height_w3, & - w3_mask - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot - real(kind=r_def), dimension(undf_wt), intent(in) :: theta - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term - real(kind=r_def), dimension(1,ndf_w3,ndf_w3), intent(in) :: basis_w3 - real(kind=r_def), dimension(1,ndf_wt,ndf_w3), intent(in) :: basis_wt - - ! Internal variables - integer(kind=i_def) :: k, df, dft, df3, eos_index_m1 - real(kind=r_def), dimension(ndf_w3) :: rho_e - real(kind=r_def), dimension(ndf_wt) :: theta_moist_e - real(kind=r_def) :: rho_cell, theta_moist, dz - - ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) - if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then - return - end if - - ! Compute exner from eqn of state at vertical-level eos_index - ! Levels in the kernel are numbered 0 to nlayers-1, rather than 1 to nlayers - - eos_index_m1 = eos_index - 1 - - k = eos_index_m1 - - do df3 = 1, ndf_w3 - rho_e( df3 ) = rho( map_w3(df3) + k) - end do - - do dft = 1, ndf_wt - theta_moist_e( dft ) = theta( map_wt(dft) + k) * moist_dyn_gas( map_wt(dft) + k ) - end do - - do df = 1, ndf_w3 - - rho_cell = 0.0_r_def - do df3 = 1, ndf_w3 - rho_cell = rho_cell + rho_e( df3 )*basis_w3( 1,df3,df ) - end do - - theta_moist = 0.0_r_def - do dft = 1, ndf_wt - theta_moist = theta_moist + theta_moist_e( dft )*basis_wt( 1,dft,df ) - end do - - exner( map_w3(df) + k ) = ( rd*rho_cell*theta_moist/p_zero ) & - **( kappa/( 1.0_r_def-kappa ) ) - end do - - ! Exner on other levels from hydrostatic balance - ! Compute exner at level k-1 from exner at level k plus other terms. - ! Add the coriolis term using the bottom face (B) from the cell at level k-1 - do k = eos_index_m1, 1, -1 - - dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) - theta_moist = moist_dyn_gas( map_wt(1) + k ) * theta( map_wt(1) + k ) / & - moist_dyn_tot( map_wt(1) + k ) - exner( map_w3(1) + k - 1 ) = exner( map_w3 (1) + k ) & - + ( gravity - coriolis_term( map_wt(1) + k ) ) * dz & - / ( cp * theta_moist ) - - end do - - do k = eos_index_m1, nlayers-2 - - dz = height_w3( map_w3(1) + k + 1 ) - height_w3( map_w3(1) + k ) - theta_moist = moist_dyn_gas( map_wt(1) + k + 1 ) * theta( map_wt(1) + k + 1) / & - moist_dyn_tot( map_wt(1) + k + 1) - exner( map_w3(1) + k + 1 ) = exner( map_w3 (1) + k ) & - - ( gravity - coriolis_term( map_wt(1) + k + 1 ) ) * dz & - / ( cp * theta_moist ) - - end do - -end subroutine hstat_cori_g_const_code - -end module hstat_cori_g_const_kernel_mod From a39118b904a6036e12749c80189d7a8e472e5c5e Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Wed, 3 Jun 2026 10:38:00 +0100 Subject: [PATCH 25/27] Remove experimental code --- .../physics/map_fd_to_prognostics_alg_mod.x90 | 5 +- .../regrav_geopot_kernel_mod.F90 | 75 ++---- .../hydro_shallow_to_deep_kernel_mod_test.pf | 238 ------------------ .../regrav_geopot_kernel_mod_test.pf | 102 ++++++++ 4 files changed, 124 insertions(+), 296 deletions(-) delete mode 100644 science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf create mode 100644 science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 0dbf9af6a..f5258d133 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -61,6 +61,7 @@ module map_fd_to_prognostics_alg_mod geometry_spherical, & topology, & topology_fully_periodic + use extrusion_config_mod, only: planet_radius use formulation_config_mod, only: rotating, shallow use finite_element_config_mod, only: coord_system, coord_system_native use initialization_config_mod, only: model_eos_height, & @@ -533,8 +534,8 @@ contains X_times_Y( temperature, exner_in_wth, theta), & ! Interpolate temperature from heights corresponding to ! geopotential height onto heights of the model grid. - regrav_geopot_kernel_type( temperature, exner, & - height_wth, mask ), & + regrav_geopot_kernel_type( temperature, height_wth, & + planet_radius ), & ! Find exner at all levels by integrating hydrostatic ! equation upwards from level 1, preserving absolute ! temperature. Gravity is height-dependent. diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 index ea7d1c8f4..504e15d8c 100644 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -18,10 +18,10 @@ module regrav_geopot_kernel_mod use argument_mod, only: arg_type, & GH_FIELD, GH_REAL, & GH_READ, GH_READWRITE, & + GH_SCALAR, & CELL_COLUMN use constants_mod, only: r_def, i_def -use extrusion_config_mod, only: planet_radius -use fs_continuity_mod, only: Wtheta, W3 +use fs_continuity_mod, only: Wtheta use kernel_mod, only: kernel_type implicit none @@ -33,11 +33,10 @@ module regrav_geopot_kernel_mod !------------------------------------------------------------------------------- type, public, extends(kernel_type) :: regrav_geopot_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) & /) integer :: operates_on = CELL_COLUMN contains @@ -53,61 +52,40 @@ module regrav_geopot_kernel_mod !> @param[in] nlayers Number of layers !> @param[in,out] temperature Absolute temperature field -!> @param[in,out] exner Exner pressure field !> @param[in] height_wth Height coordinate in wth -!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] planet_radius Radius of the planet !> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta !> @param[in] undf_wt Total number of degrees of freedom for wtheta !> @param[in] map_wt Dofmap for the cell at column base for wt -!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!> @param[in] undf_w3 Total number of degrees of freedom for w3 -!> @param[in] map_w3 Dofmap for the cell at column base for w3 subroutine regrav_geopot_code( nlayers, & temperature, & - exner, & height_wth, & - w3_mask, & + planet_radius, & ndf_wt, & undf_wt, & - map_wt, & - ndf_w3, & - undf_w3, & - map_w3 ) + map_wt ) implicit none ! Arguments integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & ndf_wt, & undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt real(kind=r_def), dimension(undf_wt), intent(inout) :: temperature - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth - real(kind=r_def), dimension(undf_w3), intent(in) :: w3_mask + real(kind=r_def), intent(in) :: planet_radius ! Internal variables integer(kind=i_def) :: k real(kind=r_def) :: height_phys(0:nlayers) real(kind=r_def) :: temp(0:nlayers) - ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) - ! setting exner to 1 - if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then - do k = 0, nlayers-1 - exner(map_w3(1) + k ) = 1.0_r_def - enddo - return - end if - ! Height corresponding to height_wth when latter represents geopotential height do k = 0, nlayers - height_phys(k) = planet_radius * height_wth(map_wt(1)+k) / & - ( planet_radius - height_wth(map_wt(1)+k) ) + height_phys(k) = planet_radius * height_wth(map_wt(1)+k) / & + ( planet_radius - height_wth(map_wt(1)+k) ) end do ! Temporary storage for T @@ -130,33 +108,18 @@ subroutine interp( nin, zin, fin, nout, zout, fout ) real(kind=r_def), intent(inout) :: fout(nout) integer(kind=i_def) :: kout, kin -real(kind=r_def) :: xi, zmin, zmax, f_at_min, f_at_max - -if ( zin(1) > zin(nin) ) then - zmax = zin(1) - zmin = zin(nin) - f_at_max = fin(1) - f_at_min = fin(nin) -else - zmax = zin(nin) - zmin = zin(1) - f_at_max = fin(nin) - f_at_min = fin(1) -end if - +real(kind=r_def) :: xi do kout = 1, nout - if ( zout(kout) <= zmin ) then - fout(kout) = f_at_min - else if ( zout(kout) >= zmax ) then - fout(kout) = f_at_max + if ( zout(kout) <= zin(1) ) then + fout(kout) = fin(1) + else if ( zout(kout) >= zin(nin) ) then + fout(kout) = fin(nin) else do kin = 1, nin - 1 - if ( ( zout(kout) - zin(kin) ) * & - ( zin(kin+1) - zout(kout) ) >= 0.0_r_def ) then + if ( zin(kin+1) > zout(kout) ) then xi = ( zout(kout) - zin(kin) ) / ( zin(kin+1) - zin(kin) ) - fout(kout) = ( 1.0_r_def - xi ) * log(fin(kin)) + xi * log(fin(kin+1)) - fout(kout) = exp(fout(kout)) + fout(kout) = ( 1.0_r_def - xi ) * fin(kin) + xi * fin(kin+1) end if end do end if diff --git a/science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf b/science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf deleted file mode 100644 index 4610f6287..000000000 --- a/science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf +++ /dev/null @@ -1,238 +0,0 @@ -!----------------------------------------------------------------------------- -! (C) Crown copyright 2025 Met Office. All rights reserved. -! The file LICENCE, distributed with this code, contains details of the terms -! under which the code may be used. -!----------------------------------------------------------------------------- - -!> Test the pressure initialization computation. -module hydro_shallow_to_deep_kernel_mod_test - - use constants_mod, only : i_def, r_def - use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & - get_w3_m3x3_q3x3x3_size, & - get_wtheta_m3x3_q3x3x3_size - use get_unit_test_w3nodal_basis_mod, only : get_wtheta_w3nodal_basis - use get_unit_test_m3x3_dofmap_mod, only : get_w0_m3x3_dofmap, & - get_w3_m3x3_dofmap, & - get_wtheta_m3x3_dofmap - use funit - - implicit none - - private - public :: hydro_shallow_to_deep_test_type, test_all - - @TestCase - type, extends(TestCase) :: hydro_shallow_to_deep_test_type - private - integer(i_def) :: ndf_w3, undf_w3, ndf_wt, undf_wt - integer(i_def), allocatable :: map_w3(:,:), map_wt(:,:) - real(r_def), allocatable :: basis_w3(:,:,:), basis_wt(:,:,:) - real(r_def), allocatable :: exner_data(:) - real(r_def), allocatable :: temperature_data(:) - real(r_def), allocatable :: coriolis_data(:) - real(r_def), allocatable :: phi_data(:) - real(r_def), allocatable :: moist_dyn_data(:,:) - real(r_def), allocatable :: height_w3_data(:) - real(r_def), allocatable :: height_wth_data(:) - real(r_def), allocatable :: mask_w3_data(:) - contains - procedure setUp - procedure tearDown - procedure test_all - end type hydro_shallow_to_deep_test_type - - real(r_def), parameter :: gravity = 10.0_r_def - real(r_def), parameter :: radius = 6000000.0_r_def - real(r_def), parameter :: omega = 8.0E-5_r_def - real(r_def), parameter :: cp = 1000.0_r_def - real(r_def), parameter :: scaling = 1.0_r_def - integer(i_def), parameter :: eos_level = 2 - -contains - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic - use finite_element_config_mod, only : cellshape_quadrilateral - use feign_config_mod, only : feign_base_mesh_config, & - feign_planet_config - - implicit none - - class(hydro_shallow_to_deep_test_type), intent(inout) :: this - - real(r_def), parameter :: dx = 6000.0_r_def, & - dy = 1000.0_r_def, & - dz = 2000.0_r_def - - integer(i_def) :: ncells, icell - integer(i_def) :: dim_space, dim_space_diff - integer(i_def) :: nqp_h, nqp_v - integer(i_def) :: nlayers, i, j, k - - call feign_base_mesh_config & - ( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_planar, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., & - f_lat_deg=0.0_r_def ) - - nlayers=3 - - call get_wtheta_w3nodal_basis(this%basis_wt) - allocate(this%basis_w3(1,1,1)) - this%basis_w3(:,:,:) = 1.0_r_def - - call get_w3_m3x3_dofmap( this%map_w3, nlayers ) - call get_w3_m3x3_q3x3x3_size( this%ndf_w3, this%undf_w3, ncells, & - dim_space, dim_space_diff, & - nqp_h, nqp_v, & - nlayers ) - - call get_wtheta_m3x3_dofmap( this%map_wt, nlayers ) - call get_wtheta_m3x3_q3x3x3_size( this%ndf_wt, this%undf_wt, ncells, & - dim_space, dim_space_diff, & - nqp_h, nqp_v, & - nlayers ) - - allocate(this%height_w3_data(this%undf_w3)) - icell = 1 - do j = 1,3 - do i = 1,3 - do k = 0,nlayers-1 - this%height_w3_data(this%map_w3(1,icell)+k) = real(k+0.5_r_def)*dz - end do - icell = icell + 1 - end do - end do - - allocate(this%height_wth_data(this%undf_wt)) - icell = 1 - do j = 1,3 - do i = 1,3 - do k = 0,nlayers - this%height_wth_data(this%map_wt(1,icell)+k) = real(k)*dz - end do - icell = icell + 1 - end do - end do - - ! Create the data - allocate(this%exner_data(this%undf_w3)) - allocate(this%phi_data(this%undf_w3)) - allocate(this%temperature_data(this%undf_wt)) - allocate(this%coriolis_data(this%undf_wt)) - allocate(this%moist_dyn_data(this%undf_wt, 3)) - allocate(this%mask_w3_data(this%undf_w3)) - - end subroutine setUp - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) - - use configuration_mod, only : final_configuration - - implicit none - - class(hydro_shallow_to_deep_test_type), intent(inout) :: this - - deallocate(this%basis_w3) - deallocate(this%basis_wt) - deallocate(this%map_w3) - deallocate(this%map_wt) - deallocate(this%exner_data) - deallocate(this%phi_data) - deallocate(this%temperature_data) - deallocate(this%coriolis_data) - deallocate(this%moist_dyn_data) - deallocate(this%height_w3_data) - deallocate(this%mask_w3_data) - - call final_configuration() - - end subroutine tearDown - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @test - subroutine test_all( this ) - - use hydro_shallow_to_deep_kernel_mod, only : hydro_shallow_to_deep_code - - implicit none - - class(hydro_shallow_to_deep_test_type), intent(inout) :: this - - real(r_def), parameter :: tol = 1.0e-6_r_def - real(r_def), parameter :: mv = 0.01_r_def ! set 10 g/kg of vapour mixing ratio - real(r_def) :: recip_epsilon - real(r_def) :: answer - integer(i_def) :: nlayers, cell, k - - nlayers=3 - cell=1 - - recip_epsilon = 461.51_r_def/300.0_r_def - ! With u=50m/s, theta=800K, rho=10^-5 - ! coriolis term = 2 \Omega u cos(lat) = 2 x 8e-5 x 50 x cos 45 = 5e-3 - this%exner_data(:) = 1.0_r_def - this%temperature_data(:) = 300.0_r_def - this%coriolis_data(:) = 5.0e-3_r_def - this%mask_w3_data(:) = 1.0_r_def - this%moist_dyn_data(:,1) = 1.0_r_def + recip_epsilon*mv - this%moist_dyn_data(:,2) = 1.0_r_def + mv - this%moist_dyn_data(:,3) = 1.0_r_def - this%phi_data(:) = 1.5_r_def * gravity * this%height_w3_data(:) - call hydro_shallow_to_deep_code( nlayers, & - this%exner_data, & - this%temperature_data, & - this%coriolis_data, & - this%moist_dyn_data(:,1), & - this%moist_dyn_data(:,2), & - this%moist_dyn_data(:,3), & - this%phi_data, & - this%height_w3_data, & - this%height_wth_data, & - this%mask_w3_data, & - gravity, & - cp, & - eos_level, & - this%ndf_w3, this%undf_w3, & - this%map_w3(:,cell), this%basis_w3, & - this%ndf_wt, this%undf_wt, & - this%map_wt(:,cell), this%basis_wt & - ) - - ! Checks hydrostatic plus coriolis balance using g (bottom level) - ! T_moist = T ( 1+ mv * recip_epsilon )/(1 + mv) - ! = 300 * (1 + (461.51/300) *0.01 )/ (1 + 0.01) - ! = 301.59910891089106 - ! exner_{k} = exner_{k+1} * ( cp T - (g - coriolis) * (1-w) ) / - ! ( cp T + (g - coriolis) * w ) - ! = 1 * ( (1000 * 301.59910891089106 ) + (10 - 0.005) * 1000 )/ - ! ( (1000 * 301.59910891089106 ) - (10 - 0.005) * 1000 ) - ! = 1.0685518461130759 - - k=0 - answer = 1.068551846113075_r_def - @assertEqual(answer, this%exner_data(this%map_w3(1, cell)+k), tol) - - ! Checks hydrostatic plus coriolis balance using phi (next level up) - ! exner_{k+1} = exner_{k} * ( cp T + (1.5g - coriolis) * (w) ) / - ! ( cp T - (1.5g - coriolis) * (1-w) ) - ! = 1.0685518461130759 * - ! ( (1000 * 301.59910891089106 ) - (15 - 0.005) * 1000 )/ - ! ( (1000 * 301.59910891089106 ) + (15 - 0.005) * 1000 ) - ! = 0.9673311696602781 - - k=1 - answer = 0.9673311696602781_r_def - @assertEqual(answer, this%exner_data(this%map_w3(1, cell)+k), tol) - - end subroutine test_all - -end module hydro_shallow_to_deep_kernel_mod_test diff --git a/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf b/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf new file mode 100644 index 000000000..5a5920f22 --- /dev/null +++ b/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf @@ -0,0 +1,102 @@ +!----------------------------------------------------------------------------- +! (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. +!----------------------------------------------------------------------------- +!> Test of interpolation from constant gravity geopotential height to true +!> geopotentail height. +!> +module regrav_geopot_kernel_mod_test + +use constants_mod, only: i_def, r_def +use, intrinsic :: iso_fortran_env, only: real64 +use funit + +implicit none + +private +public :: regrav_geopot_test_type, test_all + +@TestCase +type, extends(TestCase) :: regrav_geopot_test_type + private +contains + procedure setUp + procedure tearDown + procedure test_all +end type regrav_geopot_test_type + +integer(i_def), parameter :: nlayers = 3_i_def +integer(i_def), parameter :: ndf_wt = 2_i_def +integer(i_def), parameter :: undf_wt = nlayers + 1_i_def +real(r_def), parameter :: height_domain = 40.0e3_r_def +real(r_def), parameter :: T_surf = 300.0_r_def +real(r_def), parameter :: dTdz = -2.0e-3_r_def +real(r_def), parameter :: planet_radius = 6.4e6 + +integer(i_def), allocatable :: map_wt(:) +real(r_def), allocatable :: temperature(:) +real(r_def), allocatable :: height_wt(:) +real(r_def), allocatable :: height_phys(:) + +contains + + subroutine setUp( this ) + + implicit none + + class(regrav_geopot_test_type), intent(inout) :: this + + integer(i_def) :: k + + allocate( map_wt(ndf_wt) ) + allocate( temperature(undf_wt) ) + allocate( height_wt(undf_wt), height_phys(undf_wt) ) + + ! Create the model grid + map_wt(:) = (/ 1_i_def, 2_i_def /) + + do k = 0, nlayers + height_wt(map_wt(1) + k) = k * height_domain / nlayers + height_phys(map_wt(1) + k) = height_wt(map_wt(1) + k) / & + ( 1.0_r_def - height_wt(map_wt(1) + k) / & + planet_radius) + temperature(map_wt(1) + k) = T_surf + dTdz * height_phys(map_wt(1) + k) + end do + + end subroutine setUp + + subroutine tearDown( this ) + + implicit none + + class(regrav_geopot_test_type), intent(inout) :: this + + deallocate( map_wt ) + deallocate( temperature ) + deallocate( height_wt, height_phys ) + + end subroutine tearDown + + @Test + subroutine test_all( this ) + + use regrav_geopot_kernel_mod, only: regrav_geopot_code + + implicit none + + integer(i_def) :: k + real(r_def) :: target_value, tolerance + + call regrav_geopot_code( nlayers, temperature, height_wt, planet_radius, & + ndf_wt, undf_wt, map_wt ) + + do k = 0, nlayers + target_value = T_surf + height_wt(map_wt(1)+k) * dTdz + tolerance = 10.0_r_def * spacing(target_value) + @assertEqual(target_value, temperature(map_wt(1)+k), tolerance) + end do + + end subroutine test_all + +end module regrav_geopot_kernel_mod_test From 369a08e130e13611b03d5ddbaa8a4314fc143684 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Wed, 3 Jun 2026 13:33:06 +0100 Subject: [PATCH 26/27] Changes for developer tests --- .../algorithm/physics/map_fd_to_prognostics_alg_mod.x90 | 2 +- .../source/kernel/initialisation/regrav_geopot_kernel_mod.F90 | 2 +- .../kernel/initialisation/regrav_geopot_kernel_mod_test.pf | 4 +++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index f5258d133..e945c09af 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -468,7 +468,7 @@ contains ! Compute absolute temperature T = theta * exner X_times_Y( temperature, exner_in_wth, theta), & - ! Find exner at level 1 by integrating hydrostatic equation + ! Find exner at level 1 by integrating hydrostatic equation ! downwards from eos_level, preserving theta. Gravity ! assumed constant. hydrostatic_coriolis_kernel_type( exner, rho, theta, & diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 index 504e15d8c..498b323ab 100644 --- a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -6,7 +6,7 @@ !> @brief Convert temperature on geopotential height to temperature on height !> @details The input temperature is assumed to have been generated by a model -!> with constant gravity. For such a model physical height and +!> with constant gravity. For such a model physical height and !> geopotential height are identical. In this kernel gravity is assumed !> to vary correctly with height, so that physical and geopotential !> heights are no longer identical. The physical height, height_phys, diff --git a/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf b/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf index 5a5920f22..e950751f4 100644 --- a/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf @@ -32,7 +32,7 @@ integer(i_def), parameter :: undf_wt = nlayers + 1_i_def real(r_def), parameter :: height_domain = 40.0e3_r_def real(r_def), parameter :: T_surf = 300.0_r_def real(r_def), parameter :: dTdz = -2.0e-3_r_def -real(r_def), parameter :: planet_radius = 6.4e6 +real(r_def), parameter :: planet_radius = 6.4e6_r_def integer(i_def), allocatable :: map_wt(:) real(r_def), allocatable :: temperature(:) @@ -85,6 +85,8 @@ contains implicit none + class(regrav_geopot_test_type), intent(inout) :: this + integer(i_def) :: k real(r_def) :: target_value, tolerance From e548b9bb87ef409956ece337946d89572c924fd0 Mon Sep 17 00:00:00 2001 From: Chris Smith <256821015+mo-cjsmith@users.noreply.github.com> Date: Wed, 3 Jun 2026 16:02:45 +0100 Subject: [PATCH 27/27] Put scalar in kernel arg list --- .../physics/map_fd_to_prognostics_alg_mod.x90 | 6 ++--- .../regrav_isotherm_kernel_mod.F90 | 25 ++++++++++--------- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index e945c09af..af6cdd9c0 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -487,7 +487,7 @@ contains moist_dyn(total_mass), & geopotential, & height_w3, height_wth, & - mask ) ) + mask, cp ) ) case( regravitate_isotherm3 ) call invoke( & ! Find exner below eos_level by integrating hydrostatic @@ -515,7 +515,7 @@ contains moist_dyn(total_mass), & geopotential, & height_w3, height_wth, & - mask ) ) + mask, cp ) ) case( regravitate_geopot ) call invoke( & ! Find exner below eos_level by integrating hydrostatic @@ -545,7 +545,7 @@ contains moist_dyn(total_mass), & geopotential, & height_w3, height_wth, & - mask ) ) + mask, cp ) ) end select call invoke( & ! Recompute exner in Wtheta (required to calculate theta) diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 index 7ba70c7aa..556750dc3 100644 --- a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -32,16 +32,17 @@ module regrav_isotherm_kernel_mod !------------------------------------------------------------------------------- type, public, extends(kernel_type) :: regrav_isotherm_kernel_type private - type(arg_type) :: meta_args(9) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(10) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) /) integer :: operates_on = CELL_COLUMN contains @@ -81,6 +82,7 @@ subroutine regrav_isotherm_code( nlayers, & height_w3, & height_wth, & w3_mask, & + cp, & ndf_w3, & undf_w3, & map_w3, & @@ -88,8 +90,6 @@ subroutine regrav_isotherm_code( nlayers, & undf_wt, & map_wt ) - use planet_config_mod, only: gravity, cp - implicit none ! Arguments @@ -111,6 +111,7 @@ subroutine regrav_isotherm_code( nlayers, & real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & height_wth real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), intent(in) :: cp ! Internal variables integer(kind=i_def) :: k