Initialise with true gravity

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 |
+| ericaneininger      | Erica Neininger    | Met Office                       | 2026-03-02 |
+| mo-cjsmith          | Chris Smith        | Met Office                       | 2026-05-22 |

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
@@ -21,6 +22,8 @@ module pres_lev_diags_alg_mod
   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 extrusion_config_mod,    only: planet_radius
+  use formulation_config_mod,  only: shallow
   use planet_constants_mod,    only: ex_power
 
   implicit none
@@ -65,8 +68,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,12 +264,13 @@ 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,     &
-                                          ex_power)
+      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,        &
+                                          planet_radius, ex_power,           &
+                                          shallow)
       if (plev_geopot_flag) call plev_geopot%write_field()
 
       if (plev_geopot_clim_flag) then

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
@@ -65,7 +66,9 @@ 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]     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
@@ -85,7 +88,9 @@ subroutine geo_on_pres_code(nlayers,                         &
                               kappa,                           &
                               cp,                              &
                               gravity,                         &
+                              planet_radius,                   &
                               ex_power,                        &
+                              shallow,                         &
                               ndf_in, undf_in, map_in,         &
                               ndf_wth, undf_wth, map_wth,      &
                               ndf_out, undf_out, map_out)
@@ -110,12 +115,15 @@ 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
+    logical(kind=l_def), intent(in) :: shallow
 
     ! 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 +146,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 +158,18 @@ 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))) &
+        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) / &
                                    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 +181,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

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,14 +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, 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
+      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
@@ -1007,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()
@@ -1033,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
       !
@@ -1070,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, 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

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,12 +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 :: 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
@@ -105,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
@@ -129,7 +133,9 @@ contains
                          kappa,          &
                          cp,             &
                          gravity,        &
+                         planet_radius,  &
                          ex_power,       &
+                         shallow,        &
                          ndf_in,         &
                          undf_in,        &
                          this%map_in,    &

rose-stem/app/lfric_atm/opt/rose-app-regrav.conf

@@ -0,0 +1,5 @@
+[namelist:formulation]
+shallow=.false.
+
+[namelist:initialization]
+regravitate='geopot'

rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc

@@ -50,16 +50,16 @@
         "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,
         "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" %}