From d0281ac38193fed94467668024d8a7fc76ff0e34 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Mon, 13 Jul 2026 11:40:25 +0100 Subject: [PATCH 1/9] Add consistent metrics for the wind computation --- .../opt/rose-app-consistent_metric.conf | 2 + .../gungho_model/tasks_gungho_model.cylc | 16 ++ .../site/meto/groups/groups_gungho_model.cylc | 5 + .../control/consistent_wind_alg_mod.X90 | 137 ++++++++++++++ .../gungho_transport_control_alg_mod.X90 | 13 +- .../mol/wt_advective_update_alg_mod.x90 | 51 +----- .../configuration/check_configuration_mod.F90 | 5 - .../mol/consistent_wind_kernel_mod.F90 | 63 +++---- .../horizontal_cubic_sl_metric_kernel_mod.F90 | 171 ++++++++++++++++++ 9 files changed, 382 insertions(+), 81 deletions(-) create mode 100644 rose-stem/app/gungho_model/opt/rose-app-consistent_metric.conf create mode 100644 science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 create mode 100644 science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 diff --git a/rose-stem/app/gungho_model/opt/rose-app-consistent_metric.conf b/rose-stem/app/gungho_model/opt/rose-app-consistent_metric.conf new file mode 100644 index 0000000000..c297b61e41 --- /dev/null +++ b/rose-stem/app/gungho_model/opt/rose-app-consistent_metric.conf @@ -0,0 +1,2 @@ +[namelist:transport] +consistent_metric=.true. diff --git a/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc b/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc index 727dfdfd92..9b6fc5d3bd 100644 --- a/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc +++ b/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc @@ -399,6 +399,22 @@ "plot_str": "cartesian_mountain_plot.py diagGungho $NODAL_DATA_DIR u 3 T000900 2 2 8 schar zoom_1 lines $PLOT_DIR", }) %} +{% elif task_ns.conf_name == "schar_cart-consistent-BiP200x8-500x500" %} + + {% do task_dict.update({ + "opt_confs": ["schar_cart", "consistent_metric"], + "resolution": "BiP200x8-500x500", + "DT": 20, + "mpi_parts": 5, + "tsteps": 900, + "threads": 3, + "panel_decomp": "row", + "wallclock": 30, + "use_xios": false, + "plot_str": "cartesian_mountain_plot.py diagGungho $NODAL_DATA_DIR u 3 T000900 2 2 8 schar zoom_1 lines $PLOT_DIR", + }) %} + + {% elif task_ns.conf_name == "schar_cart-alt2-BiP100x4-1000x1000" %} {% do task_dict.update({ diff --git a/rose-stem/site/meto/groups/groups_gungho_model.cylc b/rose-stem/site/meto/groups/groups_gungho_model.cylc index 0704f11adc..840a72951e 100644 --- a/rose-stem/site/meto/groups/groups_gungho_model.cylc +++ b/rose-stem/site/meto/groups/groups_gungho_model.cylc @@ -7,6 +7,11 @@ {# Azspice Groups #} {% do site_groups.update({ + "gungho_model_schar": [ + "gungho_model_schar_cart-BiP200x8-500x500_azspice_gnu_fast-debug-64bit", + "gungho_model_schar_cart-consistent-BiP200x8-500x500_azspice_gnu_fast-debug-64bit", + "gungho_model_schar_cart-consistent-BiP200x8-500x500_azspice_gnu_full-debug-64bit", + ], "gungho_model_azspice_developer": [ "gungho_model_agnesi_hyd_cart-BiP120x8-2000x2000_azspice_gnu_fast-debug-64bit", "gungho_model_baroclinic-C24_MG_azspice_gnu_fast-debug-64bit", diff --git a/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 new file mode 100644 index 0000000000..dc1fdf939a --- /dev/null +++ b/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 @@ -0,0 +1,137 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 2026 +!----------------------------------------------------------------------------- +!> @brief Compute the consistent vertical wind + +module consistent_wind_alg_mod + + use constants_mod, only: i_def, l_def, r_def + use log_mod, only: log_event, LOG_LEVEL_INFO + use transport_metadata_collection_mod, only: transport_metadata_collection + use transport_metadata_mod, only: transport_metadata_type + use transport_controller_mod, only: transport_controller_type + use sci_geometric_constants_mod, only: get_coordinates, & + get_height_fv, & + get_face_selector_ew, & + get_face_selector_ns + use consistent_wind_kernel_mod, only: consistent_wind_kernel_type + use fs_continuity_mod, only: Wtheta, W2v, W2h + use model_clock_mod, only: model_clock_type + use combine_w2_field_kernel_mod, only: combine_w2_field_kernel_type + use split_w2_field_kernel_mod, only: split_w2_field_kernel_type + use integer_field_mod, only: integer_field_type + use field_mod, only: field_type + use function_space_mod, only: function_space_type + use function_space_collection_mod, only: function_space_collection + use wt_advective_update_alg_mod, only: hori_wt_update_alg + use transport_enumerated_types_mod, only: split_method_mol + use horizontal_cubic_sl_metric_kernel_mod, & + only: horizontal_cubic_sl_metric_kernel_type + + + implicit none + + private + + ! Contained suroutines + public :: consistent_wind_alg + + +contains + + subroutine consistent_wind_alg(transport_controller, & + rho_d_latest, & + wind_n, & + wind_np1, & + model_clock, & + outer, & + cheap_update_step) + + implicit none + + integer(kind=i_def), intent(in) :: outer + logical(kind=l_def), intent(in) :: cheap_update_step + class(model_clock_type), intent(in) :: model_clock + type(field_type), intent(in) :: rho_d_latest, wind_n, wind_np1 + type(transport_controller_type), intent(inout) :: transport_controller + + + type(field_type) :: chi3_inc, chi3 + type(field_type), pointer :: chi3_default + type(transport_metadata_type), pointer :: transport_metadata + + + integer(kind=i_def) :: mesh_id + + type(field_type), pointer :: chi(:) + type(field_type), pointer :: wind_adv_ptr + type(field_type), target :: wind_adv + type(field_type) :: wind_v, wind_h + + type(integer_field_type), pointer :: face_selector_ew + type(integer_field_type), pointer :: face_selector_ns + + type(function_space_type), pointer :: w2h_fs, w2v_fs + + integer(kind=i_def), parameter :: stencil_depth = 2 + + ! Get the height field with a big halo + mesh_id = rho_d_latest%get_mesh_id() + chi3_default => get_height_fv(Wtheta, mesh_id) + call chi3%initialise( chi3_default%get_function_space(), halo_depth = stencil_depth ) + call invoke( setval_X( chi3, chi3_default ) ) + + + ! Transport the height field + call log_event( "Computing consistent vertical wind...", LOG_LEVEL_INFO) + transport_metadata => & + transport_metadata_collection%get_transport_metadata('potential_temperature') + call chi3_inc%initialise( chi3%get_function_space() ) + + ! Split out the horzontal and vertical components of the wind + w2h_fs => function_space_collection%get_fs( wind_n%get_mesh(), & + wind_n%get_element_order_h(), & + wind_n%get_element_order_v(), W2h ) + w2v_fs => function_space_collection%get_fs( wind_n%get_mesh(), & + wind_n%get_element_order_h(), & + wind_n%get_element_order_v(), W2v ) + call wind_h%initialise( w2h_fs ) + call wind_v%initialise( w2v_fs ) + call wind_adv%initialise( wind_n%get_function_space() ) + face_selector_ew => get_face_selector_ew( mesh_id ) + face_selector_ns => get_face_selector_ns( mesh_id ) + call invoke( aX_plus_bY( wind_adv, 0.5_r_def, wind_n, 0.5_r_def, wind_np1 ), & + split_w2_field_kernel_type( wind_h, wind_v, wind_adv, & + face_selector_ew, & + face_selector_ns ) ) + + ! Compute u.grad(chi3) + if ( transport_metadata%get_horizontal_method() == split_method_mol ) then + ! With MoL + call hori_wt_update_alg(chi3_inc, wind_h, wind_h, chi3, transport_metadata) + else + ! With SL + call invoke(horizontal_cubic_sl_metric_kernel_type(chi3_inc, chi3, stencil_depth, wind_h) ) + end if + + ! Modify the vertical componenet of the transporting wind + chi => get_coordinates(mesh_id) + call invoke( consistent_wind_kernel_type(wind_v, wind_adv, & + chi3_inc, chi(3)) ) + + ! Combine the modified vertical wind with the horizontal wind to form a new consistent wind + call invoke( combine_w2_field_kernel_type( wind_adv, wind_h, wind_v, & + face_selector_ew, & + face_selector_ns ) ) + + ! Recompute controller object (departure points etc) using new consistent wind + call log_event( "Recomputing transport controller with consistent wind...", LOG_LEVEL_INFO) + wind_adv_ptr => wind_adv + call transport_controller%finalise() + call transport_controller%initialise( & + model_clock, rho_d_latest, wind_adv, wind_adv_ptr, outer, & + cheap_update_step & + ) + end subroutine consistent_wind_alg + +end module consistent_wind_alg_mod diff --git a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 index c9b1f6f50a..75f20d9ce5 100644 --- a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 +++ b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 @@ -26,7 +26,8 @@ module gungho_transport_control_alg_mod use multires_coupling_config_mod, only: coarse_aerosol_transport, & aerosol_mesh_name use r_tran_field_mod, only: r_tran_field_type - use transport_config_mod, only: transport_ageofair + use transport_config_mod, only: transport_ageofair, & + consistent_metric use transport_controller_mod, only: transport_controller_type use transport_metadata_mod, only: transport_metadata_type use transport_metadata_collection_mod, only: transport_metadata_collection_type, & @@ -211,6 +212,7 @@ contains coarse_collection_transport_alg use transport_field_mod, only: transport_field use wind_transport_alg_mod, only: wind_transport_alg + use consistent_wind_alg_mod, only: consistent_wind_alg implicit none @@ -319,6 +321,15 @@ contains ! TRANSPORT OF PROGNOSTICS: EVERY OUTER LOOP ! ======================================================================== ! + ! Modify winds to use consistent metrics if required + if (consistent_metric) call consistent_wind_alg(transport_controller, & + rho_d_latest_ptr, & + wind_for_adv, & + wind_np1, & + model_clock, & + outer, & + cheap_update_step) + ! ------------------------------------------------------------------------ ! ! Transport dry density call log_event( "Transporting density...", LOG_LEVEL_DEBUG) diff --git a/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 b/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 index 46b6b5fbce..88b2bf709b 100644 --- a/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 +++ b/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 @@ -18,39 +18,30 @@ !! size as an integer argument. module wt_advective_update_alg_mod - use constants_mod, only: r_def, i_def, l_def, str_def, & - r_tran - use sci_geometric_constants_mod, only: get_coordinates, & - get_extended_coordinates, & - get_height_fe, get_height_fv + use constants_mod, only: r_def, i_def, l_def + use sci_geometric_constants_mod, only: get_height_fv use sci_fem_constants_mod, only: get_inverse_lumped_mass_matrix_fe, & get_inverse_lumped_mass_matrix_fv use field_mod, only: field_type use r_tran_field_mod, only: r_tran_field_type - use fs_continuity_mod, only: W3, Wtheta, W2H, W2V + use fs_continuity_mod, only: Wtheta, W2H, W2V use function_space_mod, only: function_space_type use function_space_collection_mod, only: function_space_collection - use inventory_by_mesh_mod, only: inventory_by_mesh_type - use quadrature_rule_gauss_lobatto_mod, only: quadrature_rule_gauss_lobatto_type use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type use quadrature_xyoz_mod, only: quadrature_xyoz_type use quadrature_face_mod, only: quadrature_face_type - use quadrature_edge_mod, only: quadrature_edge_type use transport_config_mod, only: operators, & fv_horizontal_order, & fv_vertical_order, & operators_fv, & operators_fem, & oned_reconstruction, & - consistent_metric, & panel_edge_treatment, & panel_edge_treatment_none use transport_constants_mod, only: get_element_order_h_transport, & get_element_order_v_transport use log_mod, only: log_event, & - LOG_LEVEL_ERROR, & - LOG_LEVEL_INFO - use mesh_collection_mod, only: mesh_collection + LOG_LEVEL_ERROR use mesh_mod, only: mesh_type use poly1d_reconstruction_kernel_mod, only: poly1d_reconstruction_kernel_type use poly2d_reconstruction_kernel_mod, only: poly2d_reconstruction_kernel_type @@ -286,12 +277,9 @@ contains subroutine wt_advective_update_alg(rhs, adv_field, wind, direction, & transport_metadata, final_rk_stage) - use consistent_wind_kernel_mod, only: consistent_wind_kernel_type - use poly1d_vert_adv_kernel_mod, only: poly1d_vert_adv_kernel_type use rtheta_kernel_mod, only: rtheta_kernel_type use rtheta_bd_kernel_mod, only: rtheta_bd_kernel_type use reference_element_mod, only: reference_element_type - use mesh_mod, only: mesh_type use split_w2_field_kernel_mod, only: split_w2_field_kernel_type use integer_field_mod, only: integer_field_type use sci_geometric_constants_mod, only: get_face_selector_ew, & @@ -312,13 +300,8 @@ contains type(function_space_type), pointer :: w2h_fs, w2v_fs type(integer_field_type), pointer :: face_selector_ew type(integer_field_type), pointer :: face_selector_ns - type(field_type), pointer :: chi(:) type(field_type), pointer :: mt_lumped_inv_r_def type(r_tran_field_type) :: mt_lumped_inv_r_tran - type(r_tran_field_type) :: consistent_wind - type(r_tran_field_type) :: u_grad_z - type(field_type), pointer :: height_rdef - type(r_tran_field_type) :: height class(reference_element_type), pointer :: reference_element type(mesh_type), pointer :: mesh type(quadrature_xyoz_type) :: qr @@ -379,29 +362,9 @@ contains face_selector_ns) ) ! Horizontal update call hori_wt_update_alg(u_grad_f, wind_h, wind_h, adv_field, transport_metadata) - - ! TODO #3325: Make consistent wind be pre-computed - if ( consistent_metric ) then - ! Consistent wind will be vertical only -- initialise as wind_v - height_rdef => get_height_fv(Wtheta, mesh%get_id()) - call height%initialise( vector_space = height_rdef%get_function_space() ) - call copy_field( height_rdef, height ) - chi => get_coordinates(mesh%get_id()) - call u_grad_z%initialise( u_grad_f%get_function_space() ) - call consistent_wind%initialise( wind_v%get_function_space() ) - call invoke( setval_c(u_grad_z, 0.0_r_def), & - setval_X(consistent_wind, wind_v) ) - call hori_wt_update_alg(u_grad_z, wind_h, wind_h, height, transport_metadata) - call invoke( consistent_wind_kernel_type(consistent_wind, wind, & - u_grad_z, chi(3)) ) - ! Then vertical update uses consistent wind - call vert_wt_update_alg(u_grad_f, consistent_wind, adv_field, & - transport_metadata, final_rk_stage) - else - ! Not using consistent wind, so just do normal vertical update - call vert_wt_update_alg(u_grad_f, wind_v, adv_field, & - transport_metadata, final_rk_stage) - end if + ! Vertical update + call vert_wt_update_alg(u_grad_f, wind_v, adv_field, & + transport_metadata, final_rk_stage) case ( direction_h ) call hori_wt_update_alg(u_grad_f, wind, wind, adv_field, transport_metadata) diff --git a/science/gungho/source/configuration/check_configuration_mod.F90 b/science/gungho/source/configuration/check_configuration_mod.F90 index e2a1e8e6ec..6c1d62070d 100644 --- a/science/gungho/source/configuration/check_configuration_mod.F90 +++ b/science/gungho/source/configuration/check_configuration_mod.F90 @@ -12,7 +12,6 @@ module check_configuration_mod use transport_config_mod, only: operators, & operators_fv, & operators_fem, & - consistent_metric, & fv_horizontal_order, & fv_vertical_order, & cheap_update, & @@ -277,10 +276,6 @@ subroutine check_configuration(modeldb) end if ! Check the transport namelist - if ( geometry == geometry_spherical .and. consistent_metric) then - write( log_scratch_space, '(A)' ) 'Consistent metric option only valid for planar geometries' - call log_event( log_scratch_space, LOG_LEVEL_ERROR ) - end if any_scheme_mol = check_any_scheme_mol() if (any_scheme_mol) then ! Check that flux orders are even diff --git a/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 b/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 index fed6b05c83..c77d3e253e 100644 --- a/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 @@ -13,13 +13,14 @@ module consistent_wind_kernel_mod -use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & - GH_READWRITE, GH_READ, & - GH_BASIS, GH_DIFF_BASIS, & - CELL_COLUMN, GH_EVALUATOR -use constants_mod, only : r_def, i_def, r_tran -use fs_continuity_mod, only : Wtheta, W2, W2v, Wchi +use argument_mod, only : arg_type, func_type, & + GH_FIELD, GH_REAL, & + GH_READWRITE, GH_READ, & + GH_BASIS, GH_DIFF_BASIS, & + CELL_COLUMN, GH_EVALUATOR, & + ANY_SPACE_9 +use constants_mod, only : r_def, i_def +use fs_continuity_mod, only : Wtheta, W2, W2v use kernel_mod, only : kernel_type implicit none @@ -31,15 +32,15 @@ module consistent_wind_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the PSy layer type, public, extends(kernel_type) :: consistent_wind_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W2v), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wchi) & + type(arg_type) :: meta_args(4) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W2v), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_9) & /) - type(func_type) :: meta_funcs(2) = (/ & - func_type(W2, GH_BASIS), & - func_type(Wchi, GH_DIFF_BASIS) & + type(func_type) :: meta_funcs(2) = (/ & + func_type(W2, GH_BASIS), & + func_type(ANY_SPACE_9, GH_DIFF_BASIS) & /) integer :: operates_on = CELL_COLUMN integer :: gh_shape = GH_EVALUATOR @@ -112,10 +113,10 @@ subroutine consistent_wind_code(nlayers, & integer(kind=i_def), dimension(ndf_w2v), intent(in) :: map_w2v integer(kind=i_def), dimension(ndf_wx), intent(in) :: map_wx - real(kind=r_tran), dimension(undf_w2v), intent(inout) :: consistent_wind - real(kind=r_tran), dimension(undf_wt), intent(in) :: theta_metrics - real(kind=r_tran), dimension(undf_w2), intent(in) :: wind - real(kind=r_def), dimension(undf_wx), intent(in) :: height + real(kind=r_def), dimension(undf_w2v), intent(inout) :: consistent_wind + real(kind=r_def), dimension(undf_wt), intent(in) :: theta_metrics + real(kind=r_def), dimension(undf_w2), intent(in) :: wind + real(kind=r_def), dimension(undf_wx), intent(in) :: height real(kind=r_def), dimension(3,ndf_wx,ndf_w2v), intent(in) :: diff_basis_wx real(kind=r_def), dimension(3,ndf_w2,ndf_w2v), intent(in) :: basis_w2 @@ -123,30 +124,30 @@ subroutine consistent_wind_code(nlayers, & ! Local variables integer(kind=i_def) :: k, df, w2v_df - real(kind=r_tran) :: dz, dzdx, dzdy - real(kind=r_tran), dimension(3,0:nlayers) :: u_av + real(kind=r_def) :: dzdx, dzdy, dz + real(kind=r_def), dimension(2,0:nlayers) :: u_av ! Compute u and v averaged to w points - u_av = 0.0_r_tran + u_av = 0.0_r_def do k = 0, nlayers-1 do df = 1,4 ! Basis W2 has shape (dim, ndf_w2, ndf_w2v) = (3,6,2) - u_av(:,k) = u_av(:,k) + wind(map_w2(df)+k)*real(basis_w2(:,df,1), r_tran) - u_av(:,k+1) = u_av(:,k+1) + wind(map_w2(df)+k)*real(basis_w2(:,df,2), r_tran) + u_av(:,k) = u_av(:,k) + wind(map_w2(df)+k)*basis_w2(1:2,df,1) + u_av(:,k+1) = u_av(:,k+1) + wind(map_w2(df)+k)*basis_w2(1:2,df,2) end do end do - u_av(:,1:nlayers-1) = 0.5_r_tran*u_av(:,1:nlayers-1) + u_av(:,1:nlayers-1) = 0.5_r_def*u_av(:,1:nlayers-1) layer_loop: do k = 1, nlayers-1 ! Compute dz/dx & dz/dy on w/theta points w2v_df = 1 - dz = 0.0_r_tran - dzdx = 0.0_r_tran - dzdy = 0.0_r_tran + dz = 0.0_r_def + dzdx = 0.0_r_def + dzdy = 0.0_r_def do df = 1,ndf_wx - dzdx = dzdx + real(height(map_wx(df)+k)*diff_basis_wx(1,df,w2v_df), r_tran) - dzdy = dzdy + real(height(map_wx(df)+k)*diff_basis_wx(2,df,w2v_df), r_tran) - dz = dz + real(height(map_wx(df)+k)*diff_basis_wx(3,df,w2v_df), r_tran) + dzdx = dzdx + height(map_wx(df)+k)*diff_basis_wx(1,df,w2v_df) + dzdy = dzdy + height(map_wx(df)+k)*diff_basis_wx(2,df,w2v_df) + dz = dz + height(map_wx(df)+k)*diff_basis_wx(3,df,w2v_df) end do consistent_wind(map_w2v(w2v_df)+k) = consistent_wind(map_w2v(w2v_df)+k) & diff --git a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 new file mode 100644 index 0000000000..513ab940f1 --- /dev/null +++ b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 @@ -0,0 +1,171 @@ +!------------------------------------------------------------------------------- +! (c) Crown copyright 2023 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 Calculates the advective increments in x and y at time n+1 using +!! cubic semi-Lagrangian transport. +!> @details This kernel using cubic interpolation to solve the one-dimensional +!! advection equation in both x and y, giving advective increments +!! in both directions. This is the second part of the COSMIC splitting, +!! so the x increment works on the field previously advected in the +!! y-direction(and vice versa). +!! +!> @note This kernel only works when field is a W3/Wtheta field at lowest order. + +module horizontal_cubic_sl_metric_kernel_mod + + use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + CELL_COLUMN, GH_WRITE, & + GH_READ, GH_SCALAR, & + STENCIL, CROSS2D, GH_INTEGER, & + ANY_DISCONTINUOUS_SPACE_1 + use constants_mod, only: r_tran, i_def, l_def + use fs_continuity_mod, only: W2H + use kernel_mod, only: kernel_type + + implicit none + + private + + !----------------------------------------------------------------------------- + ! Public types + !----------------------------------------------------------------------------- + !> The type declaration for the kernel. Contains the metadata needed by the PSy layer + type, public, extends(kernel_type) :: horizontal_cubic_sl_metric_kernel_type + private + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! increment + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1, & + STENCIL(CROSS2D)), & ! z + arg_type(GH_FIELD, GH_REAL, GH_READ, W2H) & ! wind + /) + integer :: operates_on = CELL_COLUMN + contains + procedure, nopass :: horizontal_cubic_sl_metric_code + end type + + !----------------------------------------------------------------------------- + ! Contained functions/subroutines + !----------------------------------------------------------------------------- + public :: horizontal_cubic_sl_metric_code + +contains + + !> @brief Compute advective transport in x and y directions using 1D + !! Semi-Lagrangian schemes, with a cubic reconstruction. This is the + !! "outer" step of a COSMIC splitting scheme. + !> @param[in] nlayers Number of layers + !> @param[in,out] increment_x Advective increment in x direction + !> @param[in,out] increment_y Advective increment in y direction + !> @param[in] field_x Field from x direction + !> @param[in] stencil_sizes_x Sizes of the branches of the cross stencil + !> @param[in] stencil_max_x Maximum size of a cross stencil branch + !> @param[in] stencil_map_x Dofmap for the field_x stencil + !> @param[in] field_y Field from y direction + !> @param[in] stencil_sizes_y Sizes of the branches of the cross stencil + !> @param[in] stencil_max_y Maximum size of a cross stencil branch + !> @param[in] stencil_map_y Dofmap for the field_y stencil + !> @param[in] dep_pts Departure points + !> @param[in] monotone Horizontal monotone option for cubic SL + !> @param[in] ndf_wf Num of DoFs for field per cell + !> @param[in] undf_wf Num of DoFs for this partition for field + !> @param[in] map_wf Map for Wf + !> @param[in] ndf_w2h Num of DoFs for W2H per cell + !> @param[in] undf_w2h Num of DoFs for this partition for W2H + !> @param[in] map_w2h Map for W2H + subroutine horizontal_cubic_sl_metric_code( nlayers, & + increment, & + z, & + stencil_sizes, & + stencil_max, & + stencil_map, & + wind, & + ndf_wf, & + undf_wf, & + map_wf, & + ndf_w2h, & + undf_w2h, & + map_w2h ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: undf_wf + integer(kind=i_def), intent(in) :: ndf_wf + integer(kind=i_def), intent(in) :: undf_w2h + integer(kind=i_def), intent(in) :: ndf_w2h + integer(kind=i_def), intent(in) :: stencil_max + integer(kind=i_def), intent(in) :: stencil_sizes(4) + + ! Arguments: Maps + integer(kind=i_def), intent(in) :: map_wf(ndf_wf) + integer(kind=i_def), intent(in) :: map_w2h(ndf_w2h) + integer(kind=i_def), intent(in) :: stencil_map(ndf_wf,stencil_max,4) + + ! Arguments: Fields + real(kind=r_tran), intent(inout) :: increment(undf_wf) + real(kind=r_tran), intent(in) :: z(undf_wf) + real(kind=r_tran), intent(in) :: wind(undf_w2h) + + ! Local scalars + integer(kind=i_def) :: k, kp, km + + real(kind=r_tran) :: dzdx, dzdy + real(kind=r_tran) :: z_l, z_r, up, um, vp, vm + + ! Interpolation coefficients + real(kind=r_tran), parameter :: b0 = -1.0/6.0 + real(kind=r_tran), parameter :: b1 = 5.0/6.0 + real(kind=r_tran), parameter :: b2 = 2.0/6.0 + real(kind=r_tran), parameter :: c1 = 2.0/6.0 + real(kind=r_tran), parameter :: c2 = 5.0/6.0 + real(kind=r_tran), parameter :: c3 = -1.0/6.0 + + do k = 0, nlayers + km = max(0, k-1) + kp = min(nlayers-1, k) + + um = 0.5_r_tran*(wind(map_w2h(1)+km) + wind(map_w2h(1)+kp)) + up = 0.5_r_tran*(wind(map_w2h(3)+km) + wind(map_w2h(3)+kp)) + vm = -0.5_r_tran*(wind(map_w2h(2)+km) + wind(map_w2h(2)+kp)) + vp = -0.5_r_tran*(wind(map_w2h(4)+km) + wind(map_w2h(4)+kp)) + + ! dzdx + ! Compute upwind Z on the left and right sides of the cell + if ( um > 0.0_r_tran ) then + z_l = b0*z(stencil_map(1,3,1)+k) + b1*z(stencil_map(1,2,1)+k) + b2*z(stencil_map(1,1,1)+k) + else + z_l = c1*z(stencil_map(1,2,1)+k) + c2*z(stencil_map(1,1,1)+k) + c3*z(stencil_map(1,2,3)+k) + end if + if ( up > 0.0_r_tran ) then + z_r = b0*z(stencil_map(1,2,1)+k) + b1*z(stencil_map(1,1,1)+k) + b2*z(stencil_map(1,2,3)+k) + else + z_r = c1*z(stencil_map(1,1,3)+k) + c2*z(stencil_map(1,2,3)+k) + c3*z(stencil_map(1,3,3)+k) + end if + + dzdx = (z_r- z_l) + + ! dzdy + ! Compute upwind Z on the left and right sides of the cell + if ( vm > 0.0_r_tran ) then + z_l = b0*z(stencil_map(1,3,2)+k) + b1*z(stencil_map(1,2,2)+k) + b2*z(stencil_map(1,1,2)+k) + else + z_l = c1*z(stencil_map(1,2,2)+k) + c2*z(stencil_map(1,1,2)+k) + c3*z(stencil_map(1,2,4)+k) + end if + if ( vp > 0.0_r_tran ) then + z_r = b0*z(stencil_map(1,2,2)+k) + b1*z(stencil_map(1,1,2)+k) + b2*z(stencil_map(1,2,4)+k) + else + z_r = c1*z(stencil_map(1,1,4)+k) + c2*z(stencil_map(1,2,4)+k) + c3*z(stencil_map(1,3,4)+k) + end if + + dzdy = (z_r- z_l) + + increment(map_wf(1)+k) = 0.5_r_tran*(um+up)*dzdx + 0.5_r_tran*(vm+vp)*dzdy + end do + + end subroutine horizontal_cubic_sl_metric_code + +end module horizontal_cubic_sl_metric_kernel_mod From 0879e92c724b67122193381893233c6de6fb6431 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Mon, 13 Jul 2026 15:38:01 +0100 Subject: [PATCH 2/9] Undo changes to fs argument obtained from main --- .../mol/consistent_wind_kernel_mod.F90 | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 b/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 index c77d3e253e..70af8e39c5 100644 --- a/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 @@ -17,10 +17,9 @@ module consistent_wind_kernel_mod GH_FIELD, GH_REAL, & GH_READWRITE, GH_READ, & GH_BASIS, GH_DIFF_BASIS, & - CELL_COLUMN, GH_EVALUATOR, & - ANY_SPACE_9 + CELL_COLUMN, GH_EVALUATOR use constants_mod, only : r_def, i_def -use fs_continuity_mod, only : Wtheta, W2, W2v +use fs_continuity_mod, only : Wtheta, W2, W2v, Wchi use kernel_mod, only : kernel_type implicit none @@ -32,15 +31,15 @@ module consistent_wind_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the PSy layer type, public, extends(kernel_type) :: consistent_wind_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W2v), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_9) & + type(arg_type) :: meta_args(4) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W2v), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wchi) & /) - type(func_type) :: meta_funcs(2) = (/ & - func_type(W2, GH_BASIS), & - func_type(ANY_SPACE_9, GH_DIFF_BASIS) & + type(func_type) :: meta_funcs(2) = (/ & + func_type(W2, GH_BASIS), & + func_type(Wchi, GH_DIFF_BASIS) & /) integer :: operates_on = CELL_COLUMN integer :: gh_shape = GH_EVALUATOR From ee531e8e4507b0b948353a614bc3aa880dfb4161 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Mon, 13 Jul 2026 15:39:21 +0100 Subject: [PATCH 3/9] Undo changes to fs argument obtained from main --- .../transport/mol/consistent_wind_kernel_mod.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 b/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 index 70af8e39c5..6c8e7c1a34 100644 --- a/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/mol/consistent_wind_kernel_mod.F90 @@ -13,12 +13,12 @@ module consistent_wind_kernel_mod -use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & - GH_READWRITE, GH_READ, & - GH_BASIS, GH_DIFF_BASIS, & +use argument_mod, only : arg_type, func_type, & + GH_FIELD, GH_REAL, & + GH_READWRITE, GH_READ, & + GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_EVALUATOR -use constants_mod, only : r_def, i_def +use constants_mod, only : r_def, i_def, r_tran use fs_continuity_mod, only : Wtheta, W2, W2v, Wchi use kernel_mod, only : kernel_type From e42a9055832d4433c59a54dd73293759d180465d Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Mon, 17 Aug 2026 10:57:46 +0100 Subject: [PATCH 4/9] Update to metric kernel for regions near stencil boundaries --- rose-stem/app/lfric_atm/rose-app.conf | 2 +- .../site/meto/groups/groups_gungho_model.cylc | 6 +- .../site/meto/groups/groups_lfric_atm.cylc | 2 + .../control/consistent_wind_alg_mod.X90 | 15 ++ .../mol/wt_advective_update_alg_mod.x90 | 2 +- .../sl/horizontal_cubic_sl_kernel_mod.F90 | 2 +- .../horizontal_cubic_sl_metric_kernel_mod.F90 | 68 ++++----- ...izontal_cubic_sl_metric_kernel_mod_test.pf | 134 ++++++++++++++++++ 8 files changed, 190 insertions(+), 41 deletions(-) create mode 100644 science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf diff --git a/rose-stem/app/lfric_atm/rose-app.conf b/rose-stem/app/lfric_atm/rose-app.conf index bd907a5518..34a4cf5d9b 100644 --- a/rose-stem/app/lfric_atm/rose-app.conf +++ b/rose-stem/app/lfric_atm/rose-app.conf @@ -1380,7 +1380,7 @@ cfl_mol_1d_stab=1.0 cfl_mol_2d_stab=1.0 cfl_mol_3d_stab=1.0 cheap_update=.false. -consistent_metric=.false. +consistent_metric=.true. dep_pt_stencil_extent=7 dry_field_name='density' enforce_min_value=.false.,.false.,.false.,.true.,.true.,.true. diff --git a/rose-stem/site/meto/groups/groups_gungho_model.cylc b/rose-stem/site/meto/groups/groups_gungho_model.cylc index 840a72951e..d05e190e9b 100644 --- a/rose-stem/site/meto/groups/groups_gungho_model.cylc +++ b/rose-stem/site/meto/groups/groups_gungho_model.cylc @@ -7,11 +7,6 @@ {# Azspice Groups #} {% do site_groups.update({ - "gungho_model_schar": [ - "gungho_model_schar_cart-BiP200x8-500x500_azspice_gnu_fast-debug-64bit", - "gungho_model_schar_cart-consistent-BiP200x8-500x500_azspice_gnu_fast-debug-64bit", - "gungho_model_schar_cart-consistent-BiP200x8-500x500_azspice_gnu_full-debug-64bit", - ], "gungho_model_azspice_developer": [ "gungho_model_agnesi_hyd_cart-BiP120x8-2000x2000_azspice_gnu_fast-debug-64bit", "gungho_model_baroclinic-C24_MG_azspice_gnu_fast-debug-64bit", @@ -44,6 +39,7 @@ "gungho_model_sbr_lam-n96_MG_lam_azspice_gnu_fast-debug-64bit", "gungho_model_sbr_lam-n96_MG_lam_rotate_azspice_gnu_fast-debug-64bit", "gungho_model_schar_cart-BiP200x8-500x500_azspice_gnu_fast-debug-64bit", + "gungho_model_schar_cart-consistent-BiP200x8-500x500_azspice_gnu_fast-debug-64bit", "gungho_model_schar_cart-alt2-BiP100x4-1000x1000_azspice_gnu_fast-debug-64bit", "gungho_model_semi-implicit-for-linear-C12_azspice_gnu_fast-debug-64bit", "gungho_model_shallow-hot-jupiter-C24_MG_azspice_gnu_fast-debug-64bit", diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index ea319d080d..fbaba0be9c 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -37,6 +37,8 @@ "lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit", + "lfric_atm_ral3-seuk_MG_azspice_gnu_full-debug-32bit", + "lfric_atm_ral3_ens-seuk_MG_azspice_gnu_full-debug-32bit", ], "lfric_atm_scm_azspice_developer": [ "lfric_atm_scm_gal9_bomex-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit", diff --git a/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 index dc1fdf939a..d2a2f8677d 100644 --- a/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 +++ b/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 @@ -39,6 +39,21 @@ module consistent_wind_alg_mod contains + !============================================================================= + !> @brief Compute a consistent vertical wind for the current transport step. + !> @details Modifies the vertical advecting wind so that the component that comes + !! from the horizontal wind over terrain uses the transport scheme reconstruction + !! to compute the orographic gradient instead of the finite element scheme. This + !! results in the correct cancellation with the horizontal transport of the field + !! over terrain. + !> @param[in,out] transport_controller Encapsulating object containing the + !! transport counter and precomputations + !> @param[in] rho_d_latest Dry density at the latest time level + !> @param[in] wind_n Wind at the start of the step + !> @param[in] wind_np1 Wind at the end of the step + !> @param[in] model_clock Model time information + !> @param[in] outer Outer iteration index + !> @param[in] cheap_update_step Whether a cheap controller update is used subroutine consistent_wind_alg(transport_controller, & rho_d_latest, & wind_n, & diff --git a/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 b/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 index 88b2bf709b..00e740dde9 100644 --- a/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 +++ b/science/gungho/source/algorithm/transport/mol/wt_advective_update_alg_mod.x90 @@ -271,7 +271,7 @@ contains !> @param[in,out] rhs Advective increment !> @param[in] adv_field Field to advect !> @param[in] wind Advecting wind field - !> @param[in] reversible Use the reversible reconstruction + !> @param[in] direction Direction for split transport update !> @param[in] transport_metadata Contains transport configuration options !> @param[in] final_rk_stage Whether this is the last Runge-Kutta stage subroutine wt_advective_update_alg(rhs, adv_field, wind, direction, & diff --git a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_kernel_mod.F90 b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_kernel_mod.F90 index 69a7dc26ec..1f61916b4c 100644 --- a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_kernel_mod.F90 @@ -1,5 +1,5 @@ !------------------------------------------------------------------------------- -! (c) Crown copyright 2023 Met Office. All rights reserved. +! (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. !------------------------------------------------------------------------------- diff --git a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 index 513ab940f1..f0c310541f 100644 --- a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 @@ -1,15 +1,13 @@ !------------------------------------------------------------------------------- -! (c) Crown copyright 2023 Met Office. All rights reserved. +! (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 Calculates the advective increments in x and y at time n+1 using -!! cubic semi-Lagrangian transport. -!> @details This kernel using cubic interpolation to solve the one-dimensional -!! advection equation in both x and y, giving advective increments -!! in both directions. This is the second part of the COSMIC splitting, -!! so the x increment works on the field previously advected in the -!! y-direction(and vice versa). +!> @brief Compute advective increment u.grad(z) using the Nivana scheme used in +!! the ffsl and sl transport schemes. +!> @details Compute the advective increment u.grad(z) = u*dz/dx + v*dz/dy where +!! the horizontal gradients are computed using the same reconstruction method as +!! the transport scheme (in this case Nivana). !! !> @note This kernel only works when field is a W3/Wtheta field at lowest order. @@ -53,22 +51,15 @@ module horizontal_cubic_sl_metric_kernel_mod contains - !> @brief Compute advective transport in x and y directions using 1D - !! Semi-Lagrangian schemes, with a cubic reconstruction. This is the - !! "outer" step of a COSMIC splitting scheme. + !> @brief Compute advective increment u.grad(z) using the Nivana scheme used in + !! the ffsl and sl transport schemes. !> @param[in] nlayers Number of layers - !> @param[in,out] increment_x Advective increment in x direction - !> @param[in,out] increment_y Advective increment in y direction - !> @param[in] field_x Field from x direction - !> @param[in] stencil_sizes_x Sizes of the branches of the cross stencil - !> @param[in] stencil_max_x Maximum size of a cross stencil branch - !> @param[in] stencil_map_x Dofmap for the field_x stencil - !> @param[in] field_y Field from y direction - !> @param[in] stencil_sizes_y Sizes of the branches of the cross stencil - !> @param[in] stencil_max_y Maximum size of a cross stencil branch - !> @param[in] stencil_map_y Dofmap for the field_y stencil - !> @param[in] dep_pts Departure points - !> @param[in] monotone Horizontal monotone option for cubic SL + !> @param[in,out] increment Horizontal metric increment + !> @param[in] z Height field values used in reconstruction + !> @param[in] stencil_sizes Sizes of the branches of the cross stencil + !> @param[in] stencil_max Maximum size of a cross stencil branch + !> @param[in] stencil_map Dofmap for the z stencil + !> @param[in] wind Horizontal wind/departure field on W2H !> @param[in] ndf_wf Num of DoFs for field per cell !> @param[in] undf_wf Num of DoFs for this partition for field !> @param[in] map_wf Map for Wf @@ -111,7 +102,7 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & real(kind=r_tran), intent(in) :: wind(undf_w2h) ! Local scalars - integer(kind=i_def) :: k, kp, km + integer(kind=i_def) :: k, kp, km, d, d2(4), d3(4) real(kind=r_tran) :: dzdx, dzdy real(kind=r_tran) :: z_l, z_r, up, um, vp, vm @@ -123,7 +114,18 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & real(kind=r_tran), parameter :: c1 = 2.0/6.0 real(kind=r_tran), parameter :: c2 = 5.0/6.0 real(kind=r_tran), parameter :: c3 = -1.0/6.0 - + + ! Ensure that we don't do out of the domain and if there are not enough points + ! then revert to constant reconstruction + d2(:) = 1 + d3(:) = 1 + do d = 1, 4 + if ( stencil_sizes(d) == stencil_max ) then + d2(d) = 2 + d3(d) = 3 + end if + end do + do k = 0, nlayers km = max(0, k-1) kp = min(nlayers-1, k) @@ -136,14 +138,14 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & ! dzdx ! Compute upwind Z on the left and right sides of the cell if ( um > 0.0_r_tran ) then - z_l = b0*z(stencil_map(1,3,1)+k) + b1*z(stencil_map(1,2,1)+k) + b2*z(stencil_map(1,1,1)+k) + z_l = b0*z(stencil_map(1,d3(1),1)+k) + b1*z(stencil_map(1,d2(1),1)+k) + b2*z(stencil_map(1,1,1)+k) else - z_l = c1*z(stencil_map(1,2,1)+k) + c2*z(stencil_map(1,1,1)+k) + c3*z(stencil_map(1,2,3)+k) + z_l = c1*z(stencil_map(1,d2(1),1)+k) + c2*z(stencil_map(1,1,1)+k) + c3*z(stencil_map(1,d2(3),3)+k) end if if ( up > 0.0_r_tran ) then - z_r = b0*z(stencil_map(1,2,1)+k) + b1*z(stencil_map(1,1,1)+k) + b2*z(stencil_map(1,2,3)+k) + z_r = b0*z(stencil_map(1,d2(1),1)+k) + b1*z(stencil_map(1,1,1)+k) + b2*z(stencil_map(1,d2(3),3)+k) else - z_r = c1*z(stencil_map(1,1,3)+k) + c2*z(stencil_map(1,2,3)+k) + c3*z(stencil_map(1,3,3)+k) + z_r = c1*z(stencil_map(1,1,3)+k) + c2*z(stencil_map(1,d2(3),3)+k) + c3*z(stencil_map(1,d3(3),3)+k) end if dzdx = (z_r- z_l) @@ -151,14 +153,14 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & ! dzdy ! Compute upwind Z on the left and right sides of the cell if ( vm > 0.0_r_tran ) then - z_l = b0*z(stencil_map(1,3,2)+k) + b1*z(stencil_map(1,2,2)+k) + b2*z(stencil_map(1,1,2)+k) + z_l = b0*z(stencil_map(1,d3(2),2)+k) + b1*z(stencil_map(1,d2(2),2)+k) + b2*z(stencil_map(1,1,2)+k) else - z_l = c1*z(stencil_map(1,2,2)+k) + c2*z(stencil_map(1,1,2)+k) + c3*z(stencil_map(1,2,4)+k) + z_l = c1*z(stencil_map(1,d2(2),2)+k) + c2*z(stencil_map(1,1,2)+k) + c3*z(stencil_map(1,d2(4),4)+k) end if if ( vp > 0.0_r_tran ) then - z_r = b0*z(stencil_map(1,2,2)+k) + b1*z(stencil_map(1,1,2)+k) + b2*z(stencil_map(1,2,4)+k) + z_r = b0*z(stencil_map(1,d2(2),2)+k) + b1*z(stencil_map(1,1,2)+k) + b2*z(stencil_map(1,d2(4),4)+k) else - z_r = c1*z(stencil_map(1,1,4)+k) + c2*z(stencil_map(1,2,4)+k) + c3*z(stencil_map(1,3,4)+k) + z_r = c1*z(stencil_map(1,1,4)+k) + c2*z(stencil_map(1,d2(4),4)+k) + c3*z(stencil_map(1,d3(4),4)+k) end if dzdy = (z_r- z_l) diff --git a/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf b/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf new file mode 100644 index 0000000000..ed59f8abb4 --- /dev/null +++ b/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf @@ -0,0 +1,134 @@ +!----------------------------------------------------------------------------- +! (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 Test the horizontal cubic semi-Lagrangian metric kernel + +module horizontal_cubic_sl_metric_kernel_mod_test + + use constants_mod, only : i_def, r_tran + + implicit none + +contains + + !------------------------------------------------------------------ + + @Test + subroutine horizontal_cubic_sl_metric_kernel_test( ) + + use, intrinsic :: iso_fortran_env, only: real64 + use funit + use horizontal_cubic_sl_metric_kernel_mod, only: horizontal_cubic_sl_metric_code + + implicit none + + real(kind=r_tran), parameter :: tol = 1.0e-12_r_tran + real(kind=r_tran) :: use_tol + + integer(kind=i_def), parameter :: nlayers = 1 + integer(kind=i_def), parameter :: ndf_wf = 1 + integer(kind=i_def), parameter :: ndf_w2h = 4 + integer(kind=i_def), parameter :: undf_wf = 18 + integer(kind=i_def), parameter :: undf_w2h = 4 + integer(kind=i_def), parameter :: stencil_max = 3 + + integer(kind=i_def), dimension(ndf_wf) :: map_wf + integer(kind=i_def), dimension(ndf_w2h) :: map_w2h + integer(kind=i_def), dimension(4) :: stencil_sizes + integer(kind=i_def), dimension(ndf_wf,stencil_max,4) :: stencil_map + + real(kind=r_tran), dimension(undf_wf) :: increment + real(kind=r_tran), dimension(undf_wf) :: z + real(kind=r_tran), dimension(undf_w2h) :: wind + + ! Base indices of each horizontal point for 2 levels (k=0,1): + ! west: 3->5, 2->3, center->1, east: 6->11, 7->13, + ! south: 4->7, 5->9, north: 8->15, 9->17. + integer(kind=i_def), parameter :: idx_c = 1 + integer(kind=i_def), parameter :: idx_w1 = 3 + integer(kind=i_def), parameter :: idx_w2 = 5 + integer(kind=i_def), parameter :: idx_e1 = 11 + integer(kind=i_def), parameter :: idx_e2 = 13 + integer(kind=i_def), parameter :: idx_s1 = 7 + integer(kind=i_def), parameter :: idx_s2 = 9 + integer(kind=i_def), parameter :: idx_n1 = 15 + integer(kind=i_def), parameter :: idx_n2 = 17 + + ! Get correct tolerance + if ( r_tran == real64 ) then + use_tol = tol + else + use_tol = 1.0e-6_r_tran + end if + + map_wf(1) = 1 + map_w2h(:) = (/ 1, 2, 3, 4 /) + + stencil_sizes(:) = stencil_max + stencil_map(:,:,:) = reshape( & + (/ idx_c, idx_w1, idx_w2, & ! W branch + idx_c, idx_s1, idx_s2, & ! S branch + idx_c, idx_e1, idx_e2, & ! E branch + idx_c, idx_n1, idx_n2 /), & ! N branch + [ndf_wf, stencil_max, 4] & + ) + + z(:) = 0.0_r_tran + ! x-like profile that is linear along W/E and constant along S/N. + z(idx_w2:idx_w2+1) = -2.0_r_tran + z(idx_w1:idx_w1+1) = -1.0_r_tran + z(idx_c:idx_c+1) = 0.0_r_tran + z(idx_e1:idx_e1+1) = 1.0_r_tran + z(idx_e2:idx_e2+1) = 2.0_r_tran + z(idx_s1:idx_s1+1) = 0.0_r_tran + z(idx_s2:idx_s2+1) = 0.0_r_tran + z(idx_n1:idx_n1+1) = 0.0_r_tran + z(idx_n2:idx_n2+1) = 0.0_r_tran + + ! TEST 1: Positive zonal wind => positive increment on both levels. + increment(:) = 0.0_r_tran + wind(:) = (/ 2.0_r_tran, 0.0_r_tran, 2.0_r_tran, 0.0_r_tran /) + + call horizontal_cubic_sl_metric_code( nlayers, & + increment, & + z, & + stencil_sizes, & + stencil_max, & + stencil_map, & + wind, & + ndf_wf, & + undf_wf, & + map_wf, & + ndf_w2h, & + undf_w2h, & + map_w2h ) + + @assertEqual(2.0_r_tran, increment(1), use_tol) + @assertEqual(2.0_r_tran, increment(2), use_tol) + + ! TEST 2: Negative zonal wind => negative increment on both levels. + increment(:) = 0.0_r_tran + wind(:) = (/ -2.0_r_tran, 0.0_r_tran, -2.0_r_tran, 0.0_r_tran /) + + call horizontal_cubic_sl_metric_code( nlayers, & + increment, & + z, & + stencil_sizes, & + stencil_max, & + stencil_map, & + wind, & + ndf_wf, & + undf_wf, & + map_wf, & + ndf_w2h, & + undf_w2h, & + map_w2h ) + + @assertEqual(-2.0_r_tran, increment(1), use_tol) + @assertEqual(-2.0_r_tran, increment(2), use_tol) + + end subroutine horizontal_cubic_sl_metric_kernel_test + +end module horizontal_cubic_sl_metric_kernel_mod_test From 70453095bf0bba88e6a79a95bcdbaf0f448c8907 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Wed, 26 Aug 2026 13:24:54 +0100 Subject: [PATCH 5/9] Fix code to use correct r_tran types --- .../control/consistent_wind_alg_mod.X90 | 22 +++++++--- .../horizontal_cubic_sl_metric_kernel_mod.F90 | 44 +++++++++---------- 2 files changed, 37 insertions(+), 29 deletions(-) diff --git a/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 index d2a2f8677d..3505bdaa57 100644 --- a/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 +++ b/science/gungho/source/algorithm/transport/control/consistent_wind_alg_mod.X90 @@ -21,13 +21,13 @@ module consistent_wind_alg_mod use split_w2_field_kernel_mod, only: split_w2_field_kernel_type use integer_field_mod, only: integer_field_type use field_mod, only: field_type + use r_tran_field_mod, only: r_tran_field_type use function_space_mod, only: function_space_type use function_space_collection_mod, only: function_space_collection use wt_advective_update_alg_mod, only: hori_wt_update_alg use transport_enumerated_types_mod, only: split_method_mol use horizontal_cubic_sl_metric_kernel_mod, & only: horizontal_cubic_sl_metric_kernel_type - implicit none @@ -69,13 +69,11 @@ contains class(model_clock_type), intent(in) :: model_clock type(field_type), intent(in) :: rho_d_latest, wind_n, wind_np1 type(transport_controller_type), intent(inout) :: transport_controller - type(field_type) :: chi3_inc, chi3 type(field_type), pointer :: chi3_default type(transport_metadata_type), pointer :: transport_metadata - integer(kind=i_def) :: mesh_id type(field_type), pointer :: chi(:) @@ -83,6 +81,8 @@ contains type(field_type), target :: wind_adv type(field_type) :: wind_v, wind_h + type(r_tran_field_type) :: chi3_inc_r_tran, wind_h_r_tran, chi3_r_tran + type(integer_field_type), pointer :: face_selector_ew type(integer_field_type), pointer :: face_selector_ns @@ -111,7 +111,7 @@ contains wind_n%get_element_order_h(), & wind_n%get_element_order_v(), W2v ) call wind_h%initialise( w2h_fs ) - call wind_v%initialise( w2v_fs ) + call wind_v%initialise( w2v_fs ) call wind_adv%initialise( wind_n%get_function_space() ) face_selector_ew => get_face_selector_ew( mesh_id ) face_selector_ns => get_face_selector_ns( mesh_id ) @@ -119,11 +119,19 @@ contains split_w2_field_kernel_type( wind_h, wind_v, wind_adv, & face_selector_ew, & face_selector_ns ) ) - + ! Compute u.grad(chi3) if ( transport_metadata%get_horizontal_method() == split_method_mol ) then ! With MoL - call hori_wt_update_alg(chi3_inc, wind_h, wind_h, chi3, transport_metadata) + call chi3_inc_r_tran%initialise( chi3_inc%get_function_space() ) + call wind_h_r_tran%initialise( wind_h%get_function_space() ) + call chi3_r_tran%initialise( chi3%get_function_space() ) + call invoke( real_to_real_X( chi3_inc_r_tran, chi3_inc ), & + real_to_real_X( wind_h_r_tran, wind_h ), & + real_to_real_X( chi3_r_tran, chi3 ) ) + call hori_wt_update_alg(chi3_inc_r_tran, wind_h_r_tran, wind_h_r_tran, & + chi3_r_tran, transport_metadata) + call invoke( real_to_real_X( chi3_inc, chi3_inc_r_tran ) ) else ! With SL call invoke(horizontal_cubic_sl_metric_kernel_type(chi3_inc, chi3, stencil_depth, wind_h) ) @@ -141,7 +149,7 @@ contains ! Recompute controller object (departure points etc) using new consistent wind call log_event( "Recomputing transport controller with consistent wind...", LOG_LEVEL_INFO) - wind_adv_ptr => wind_adv + wind_adv_ptr => wind_adv call transport_controller%finalise() call transport_controller%initialise( & model_clock, rho_d_latest, wind_adv, wind_adv_ptr, outer, & diff --git a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 index f0c310541f..f4b1fd3757 100644 --- a/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod.F90 @@ -19,7 +19,7 @@ module horizontal_cubic_sl_metric_kernel_mod GH_READ, GH_SCALAR, & STENCIL, CROSS2D, GH_INTEGER, & ANY_DISCONTINUOUS_SPACE_1 - use constants_mod, only: r_tran, i_def, l_def + use constants_mod, only: r_def, i_def, l_def use fs_continuity_mod, only: W2H use kernel_mod, only: kernel_type @@ -97,23 +97,23 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & integer(kind=i_def), intent(in) :: stencil_map(ndf_wf,stencil_max,4) ! Arguments: Fields - real(kind=r_tran), intent(inout) :: increment(undf_wf) - real(kind=r_tran), intent(in) :: z(undf_wf) - real(kind=r_tran), intent(in) :: wind(undf_w2h) + real(kind=r_def), intent(inout) :: increment(undf_wf) + real(kind=r_def), intent(in) :: z(undf_wf) + real(kind=r_def), intent(in) :: wind(undf_w2h) ! Local scalars integer(kind=i_def) :: k, kp, km, d, d2(4), d3(4) - real(kind=r_tran) :: dzdx, dzdy - real(kind=r_tran) :: z_l, z_r, up, um, vp, vm + real(kind=r_def) :: dzdx, dzdy + real(kind=r_def) :: z_l, z_r, up, um, vp, vm ! Interpolation coefficients - real(kind=r_tran), parameter :: b0 = -1.0/6.0 - real(kind=r_tran), parameter :: b1 = 5.0/6.0 - real(kind=r_tran), parameter :: b2 = 2.0/6.0 - real(kind=r_tran), parameter :: c1 = 2.0/6.0 - real(kind=r_tran), parameter :: c2 = 5.0/6.0 - real(kind=r_tran), parameter :: c3 = -1.0/6.0 + real(kind=r_def), parameter :: b0 = -1.0_r_def/6.0_r_def + real(kind=r_def), parameter :: b1 = 5.0_r_def/6.0_r_def + real(kind=r_def), parameter :: b2 = 2.0_r_def/6.0_r_def + real(kind=r_def), parameter :: c1 = 2.0_r_def/6.0_r_def + real(kind=r_def), parameter :: c2 = 5.0_r_def/6.0_r_def + real(kind=r_def), parameter :: c3 = -1.0_r_def/6.0_r_def ! Ensure that we don't do out of the domain and if there are not enough points ! then revert to constant reconstruction @@ -130,19 +130,19 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & km = max(0, k-1) kp = min(nlayers-1, k) - um = 0.5_r_tran*(wind(map_w2h(1)+km) + wind(map_w2h(1)+kp)) - up = 0.5_r_tran*(wind(map_w2h(3)+km) + wind(map_w2h(3)+kp)) - vm = -0.5_r_tran*(wind(map_w2h(2)+km) + wind(map_w2h(2)+kp)) - vp = -0.5_r_tran*(wind(map_w2h(4)+km) + wind(map_w2h(4)+kp)) + um = 0.5_r_def*(wind(map_w2h(1)+km) + wind(map_w2h(1)+kp)) + up = 0.5_r_def*(wind(map_w2h(3)+km) + wind(map_w2h(3)+kp)) + vm = -0.5_r_def*(wind(map_w2h(2)+km) + wind(map_w2h(2)+kp)) + vp = -0.5_r_def*(wind(map_w2h(4)+km) + wind(map_w2h(4)+kp)) ! dzdx ! Compute upwind Z on the left and right sides of the cell - if ( um > 0.0_r_tran ) then + if ( um > 0.0_r_def ) then z_l = b0*z(stencil_map(1,d3(1),1)+k) + b1*z(stencil_map(1,d2(1),1)+k) + b2*z(stencil_map(1,1,1)+k) else z_l = c1*z(stencil_map(1,d2(1),1)+k) + c2*z(stencil_map(1,1,1)+k) + c3*z(stencil_map(1,d2(3),3)+k) end if - if ( up > 0.0_r_tran ) then + if ( up > 0.0_r_def ) then z_r = b0*z(stencil_map(1,d2(1),1)+k) + b1*z(stencil_map(1,1,1)+k) + b2*z(stencil_map(1,d2(3),3)+k) else z_r = c1*z(stencil_map(1,1,3)+k) + c2*z(stencil_map(1,d2(3),3)+k) + c3*z(stencil_map(1,d3(3),3)+k) @@ -152,12 +152,12 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & ! dzdy ! Compute upwind Z on the left and right sides of the cell - if ( vm > 0.0_r_tran ) then + if ( vm > 0.0_r_def ) then z_l = b0*z(stencil_map(1,d3(2),2)+k) + b1*z(stencil_map(1,d2(2),2)+k) + b2*z(stencil_map(1,1,2)+k) else z_l = c1*z(stencil_map(1,d2(2),2)+k) + c2*z(stencil_map(1,1,2)+k) + c3*z(stencil_map(1,d2(4),4)+k) end if - if ( vp > 0.0_r_tran ) then + if ( vp > 0.0_r_def ) then z_r = b0*z(stencil_map(1,d2(2),2)+k) + b1*z(stencil_map(1,1,2)+k) + b2*z(stencil_map(1,d2(4),4)+k) else z_r = c1*z(stencil_map(1,1,4)+k) + c2*z(stencil_map(1,d2(4),4)+k) + c3*z(stencil_map(1,d3(4),4)+k) @@ -165,8 +165,8 @@ subroutine horizontal_cubic_sl_metric_code( nlayers, & dzdy = (z_r- z_l) - increment(map_wf(1)+k) = 0.5_r_tran*(um+up)*dzdx + 0.5_r_tran*(vm+vp)*dzdy - end do + increment(map_wf(1)+k) = 0.5_r_def*(um+up)*dzdx + 0.5_r_def*(vm+vp)*dzdy + end do end subroutine horizontal_cubic_sl_metric_code From 89611bb91eb341a8394d8ad05ac2223de18c7e15 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Thu, 27 Aug 2026 11:38:02 +0100 Subject: [PATCH 6/9] KGO update --- rose-stem/app/lfric_atm/rose-app.conf | 2 +- ...alt2-BiP100x4-1000x1000_azspice_gnu_fast-debug-64bit.txt | 6 +++--- ...lt2-BiP256x16-200x50_op_azspice_gnu_fast-debug-64bit.txt | 6 +++--- ...rt-alt2-BiP100x4-1000x1000_ex1a_gnu_fast-debug-64bit.txt | 6 +++--- ...m-alt2-BiP256x16-200x50_op_ex1a_gnu_fast-debug-64bit.txt | 6 +++--- 5 files changed, 13 insertions(+), 13 deletions(-) diff --git a/rose-stem/app/lfric_atm/rose-app.conf b/rose-stem/app/lfric_atm/rose-app.conf index 34a4cf5d9b..bd907a5518 100644 --- a/rose-stem/app/lfric_atm/rose-app.conf +++ b/rose-stem/app/lfric_atm/rose-app.conf @@ -1380,7 +1380,7 @@ cfl_mol_1d_stab=1.0 cfl_mol_2d_stab=1.0 cfl_mol_3d_stab=1.0 cheap_update=.false. -consistent_metric=.true. +consistent_metric=.false. dep_pt_stencil_extent=7 dry_field_name='density' enforce_min_value=.false.,.false.,.false.,.true.,.true.,.true. diff --git a/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_azspice_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_azspice_gnu_fast-debug-64bit.txt index b81805f523..5c4466503c 100644 --- a/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_azspice_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_azspice_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 40C1E430BB50B00C -Inner product checksum theta = 41F13953C045587E -Inner product checksum u = 4393B8DB189E7930 +Inner product checksum rho = 40C1E4603A077978 +Inner product checksum theta = 41F13945B83E316E +Inner product checksum u = 4393BCD0E5EC2194 diff --git a/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_azspice_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_azspice_gnu_fast-debug-64bit.txt index e2553efc31..7b2a0f8063 100644 --- a/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_azspice_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_azspice_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 40FA16E969D9ADBC -Inner product checksum theta = 42069BF1BB6B3AAE -Inner product checksum u = 42FFC30F50D36E92 +Inner product checksum rho = 40FA16E969D95AD8 +Inner product checksum theta = 42069BF1BB6C33BA +Inner product checksum u = 42FFC30F50FCC7DE diff --git a/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_ex1a_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_ex1a_gnu_fast-debug-64bit.txt index adf57f0b7a..0bf36ff95e 100644 --- a/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_ex1a_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_schar_cart-alt2-BiP100x4-1000x1000_ex1a_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 40C1E430BB50B1BB -Inner product checksum theta = 41F13953C0455941 -Inner product checksum u = 4393B8DB189EAE30 +Inner product checksum rho = 40C1E4603A075C2E +Inner product checksum theta = 41F13945B83E30E2 +Inner product checksum u = 4393BCD0E5EE0793 diff --git a/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_ex1a_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_ex1a_gnu_fast-debug-64bit.txt index 64502ae9f9..c48d141996 100644 --- a/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_ex1a_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_straka_200m-alt2-BiP256x16-200x50_op_ex1a_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 40FA16E969D78B0B -Inner product checksum theta = 42069BF1BB70B7A9 -Inner product checksum u = 42FFC30F4F0FCDA2 +Inner product checksum rho = 40FA16E969D75BB2 +Inner product checksum theta = 42069BF1BB705BB7 +Inner product checksum u = 42FFC30F4F1E8D4C From 303c68c6d37da764cbbc2e485d6b8b47f113ed54 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Fri, 28 Aug 2026 11:58:46 +0100 Subject: [PATCH 7/9] Remove added debug tests --- rose-stem/site/meto/groups/groups_lfric_atm.cylc | 2 -- 1 file changed, 2 deletions(-) diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index fbaba0be9c..ea319d080d 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -37,8 +37,6 @@ "lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit", - "lfric_atm_ral3-seuk_MG_azspice_gnu_full-debug-32bit", - "lfric_atm_ral3_ens-seuk_MG_azspice_gnu_full-debug-32bit", ], "lfric_atm_scm_azspice_developer": [ "lfric_atm_scm_gal9_bomex-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit", From 07ab27b31ad8b22be641d12ef2ca4a7e33adc183 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Fri, 18 Sep 2026 09:28:12 +0100 Subject: [PATCH 8/9] Response to coppilot reviews --- .../meto/common/suite_config_azspice.cylc | 5 +- .../gungho_transport_control_alg_mod.X90 | 27 ++++++--- ...izontal_cubic_sl_metric_kernel_mod_test.pf | 55 +++++++++++++++++++ 3 files changed, 79 insertions(+), 8 deletions(-) diff --git a/rose-stem/site/meto/common/suite_config_azspice.cylc b/rose-stem/site/meto/common/suite_config_azspice.cylc index 945d6a636b..2bc10b202c 100644 --- a/rose-stem/site/meto/common/suite_config_azspice.cylc +++ b/rose-stem/site/meto/common/suite_config_azspice.cylc @@ -7,6 +7,7 @@ {% set azspice_base = 'umask 0022 ; '~ 'module purge ; '~ + 'export PYTHONNOUSERSITE=1 ; '~ 'module use /home/users/lfricadmin/lmod' %} {% set azspice_compiler_gnu = 'module load lfric/vn3.2' %} @@ -16,7 +17,8 @@ {% set azspice_run = 'ulimit -s unlimited' %} -{% set azspice_scitools = 'module load scitools/production-os47-1' %} +{% set azspice_scitools = 'export PYTHONNOUSERSITE=1 ; '~ + 'module load scitools/production-os47-1' %} {% set azspice_tech = 'module load lfric/vn3.2' %} @@ -117,6 +119,7 @@ inherit = AZSPICE_BASE pre-script = """ module purge + export PYTHONNOUSERSITE=1 module load scitools/production-os47-1 module load um_tools/2025.10.1/openmp module load csc/2025.3.20 diff --git a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 index 4548dc3d3c..4085085715 100644 --- a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 +++ b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 @@ -455,13 +455,26 @@ contains ! ---------------------------------------------------------------------- ! if (cheap_update) then ! Transport is now with transporting wind for the whole time step, not - ! the difference in winds from the cheap update - call transport_controller%finalise() - call transport_controller%initialise( & - config, & - model_clock, rho_d_n_ptr, wind_n, wind_np1, & - outer=outer, cheap_update_step=.false. & - ) + ! the difference in winds from the cheap update so the controller needs + ! to be updated. If consistent metric is used if will be done in + ! consistent_wind_alg, otherwise it is recomputed here + if (consistent_metric) then + call consistent_wind_alg(config, & + transport_controller, & + rho_d_n_ptr, & + wind_n, & + wind_np1, & + model_clock, & + outer, & + cheap_update_step) + else + call transport_controller%finalise() + call transport_controller%initialise( & + config, & + model_clock, rho_d_n_ptr, wind_n, wind_np1, & + outer=outer, cheap_update_step=.false. & + ) + end if ! Copy dry density flux for the whole timestep into flux precomputations call transport_controller%initialise_flux_precomputations(total_dry_flux) diff --git a/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf b/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf index ed59f8abb4..8fcab89175 100644 --- a/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/transport/sl/horizontal_cubic_sl_metric_kernel_mod_test.pf @@ -129,6 +129,61 @@ contains @assertEqual(-2.0_r_tran, increment(1), use_tol) @assertEqual(-2.0_r_tran, increment(2), use_tol) + z(:) = 0.0_r_tran + ! x-like profile that is linear along S/N and constant along E/W. + z(idx_n2:idx_n2+1) = -2.0_r_tran + z(idx_n1:idx_n1+1) = -1.0_r_tran + z(idx_c:idx_c+1) = 0.0_r_tran + z(idx_s1:idx_s1+1) = 1.0_r_tran + z(idx_s2:idx_s2+1) = 2.0_r_tran + z(idx_w1:idx_w1+1) = 0.0_r_tran + z(idx_w2:idx_w2+1) = 0.0_r_tran + z(idx_e1:idx_e1+1) = 0.0_r_tran + z(idx_e2:idx_e2+1) = 0.0_r_tran + + ! TEST 3: Positive meridional wind => positive increment on both levels. + increment(:) = 0.0_r_tran + wind(:) = (/ 0.0_r_tran, 2.0_r_tran, 0.0_r_tran, 2.0_r_tran /) + + call horizontal_cubic_sl_metric_code( nlayers, & + increment, & + z, & + stencil_sizes, & + stencil_max, & + stencil_map, & + wind, & + ndf_wf, & + undf_wf, & + map_wf, & + ndf_w2h, & + undf_w2h, & + map_w2h ) + + @assertEqual(2.0_r_tran, increment(1), use_tol) + @assertEqual(2.0_r_tran, increment(2), use_tol) + + ! TEST 4: Negative meridional wind => negative increment on both levels. + increment(:) = 0.0_r_tran + wind(:) = (/ 0.0_r_tran, -2.0_r_tran, 0.0_r_tran, -2.0_r_tran /) + + call horizontal_cubic_sl_metric_code( nlayers, & + increment, & + z, & + stencil_sizes, & + stencil_max, & + stencil_map, & + wind, & + ndf_wf, & + undf_wf, & + map_wf, & + ndf_w2h, & + undf_w2h, & + map_w2h ) + + @assertEqual(-2.0_r_tran, increment(1), use_tol) + @assertEqual(-2.0_r_tran, increment(2), use_tol) + + end subroutine horizontal_cubic_sl_metric_kernel_test end module horizontal_cubic_sl_metric_kernel_mod_test From 56b3ada8791c16a16abfbfb1c9975b2df2ad4f08 Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Fri, 18 Sep 2026 09:49:49 +0100 Subject: [PATCH 9/9] Revert accidental commit of local changes to azspice config --- rose-stem/site/meto/common/suite_config_azspice.cylc | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/rose-stem/site/meto/common/suite_config_azspice.cylc b/rose-stem/site/meto/common/suite_config_azspice.cylc index 2bc10b202c..945d6a636b 100644 --- a/rose-stem/site/meto/common/suite_config_azspice.cylc +++ b/rose-stem/site/meto/common/suite_config_azspice.cylc @@ -7,7 +7,6 @@ {% set azspice_base = 'umask 0022 ; '~ 'module purge ; '~ - 'export PYTHONNOUSERSITE=1 ; '~ 'module use /home/users/lfricadmin/lmod' %} {% set azspice_compiler_gnu = 'module load lfric/vn3.2' %} @@ -17,8 +16,7 @@ {% set azspice_run = 'ulimit -s unlimited' %} -{% set azspice_scitools = 'export PYTHONNOUSERSITE=1 ; '~ - 'module load scitools/production-os47-1' %} +{% set azspice_scitools = 'module load scitools/production-os47-1' %} {% set azspice_tech = 'module load lfric/vn3.2' %} @@ -119,7 +117,6 @@ inherit = AZSPICE_BASE pre-script = """ module purge - export PYTHONNOUSERSITE=1 module load scitools/production-os47-1 module load um_tools/2025.10.1/openmp module load csc/2025.3.20