New function space options for coordinate fields

CONTRIBUTORS.md

@@ -7,3 +7,4 @@
 | jennyhickson    | Jenny Hickson    | Met Office  | 2025-12-10 |
 | mo-marqh        | Mark Hedley      | Met Office  | 2025-12-11 |
 | MatthewHambley  | Matthew Hambley  | Met Office  | 2025-12-15 |
+| tommbendall     | Thomas Bendall   | Met Office  | 2026-01-23 |

applications/coupled/example/configuration_glo.nml

@@ -21,6 +21,8 @@ element_order_h=0,
 element_order_v=0,
 rehabilitate=.true.,
 coord_order=1
+coord_order_multigrid=1
+coord_space='Wchi'
 coord_system='native'
 /
 

applications/coupled/example/configuration_lam.nml

@@ -20,6 +20,8 @@ element_order_h=0,
 element_order_v=0,
 rehabilitate=.true.,
 coord_order=1
+coord_order_multigrid=1
+coord_space='Wchi'
 coord_system='native'
 /
 

applications/io_demo/example/configuration.nml

@@ -21,6 +21,8 @@ element_order_h=0,
 element_order_v=0,
 rehabilitate=.true.,
 coord_order=1
+coord_order_multigrid=1
+coord_space='Wchi'
 coord_system='native'
 /
 

applications/lbc_demo/example/configuration.nml

@@ -21,6 +21,8 @@ coord_order=1,
 coord_system='native',
 element_order_h=0,
 element_order_v=0,
+coord_order_multigrid=1
+coord_space='Wchi'
 rehabilitate=.true.,
 /
 

applications/simple_diffusion/example/configuration.nml

@@ -21,6 +21,8 @@ element_order_h=0,
 element_order_v=0,
 rehabilitate=.true.,
 coord_order=1
+coord_order_multigrid=1
+coord_space='Wtheta'
 coord_system='native'
 /
 

applications/skeleton/example/configuration.nml

@@ -21,6 +21,8 @@ element_order_h=0,
 element_order_v=0,
 rehabilitate=.true.,
 coord_order=1
+coord_order_multigrid=1
+coord_space='Wchi'
 coord_system='native'
 /
 

components/driver/rose-meta/lfric-driver/HEAD/rose-meta.conf

@@ -326,17 +326,48 @@ values='triangle', 'quadrilateral'
 
 [namelist:finite_element=coord_order]
 compulsory=true
-description=Order of the coordinate space.
-fail-if=this < 0 ;
-help=Order of the coordinate space. If 0 is chosen, this will be the continuous
-    =W0 space, with the order set by the element order. This option is only
-    =possible for certain geometries and topologies. Otherwise, this will be
-    =a discontinuous space.
+description=Horizontal polynomial order of the coordinate space.
+fail-if=this < 1 ;
+help=Horizontal order of the polynomials used to describe the coordinate space.
+    =The vertical order is determined by the coordinate system. If a native
+    =coordinates system is chosen, the vertical coordinates are linear.
+    =If Cartesian coordinates are chosen, the vertical coordinates have the same
+    =order as the horizontal.
 !kind=default
-range=0:
+range=1:
+sort-key=Panel-A03
+type=integer
+
+[namelist:finite_element=coord_order_multigrid]
+compulsory=true
+description=Horizontal polynomial order of coordinates for multigrid meshes.
+fail-if=this < 1 ;
+help=The order of the coordinate spaces can be set separately for the prime
+    =mesh and the others in a multigrid chain. For instance it may not be
+    =necessary to use higher order coordinates on coarser meshes.
+    =The vertical order is set as for the coord_order option.
+!kind=default
+range=1:
 sort-key=Panel-A03
 type=integer
 
+[namelist:finite_element=coord_space]
+compulsory=true
+description=The finite element space to use to hold the coordinate fields.
+!enumeration=true
+fail-if=
+help=Different spaces have different continuity properties:
+    =Wchi:   fully discontinuous.
+    =Wtheta: discontinuous in the horizontal, but continuous in the vertical.
+    =        This is only appropriate for native coordinate systems, which
+    =        separate the horizontal and vertical coordinates.
+    =W0:     fully continuous. This is only appropriate for non-periodic planar
+    =        domains, or when using the 'xyz' coordinate system on spherical
+    =        domains.
+sort-key=Panel-A04
+value-titles='Wchi', 'Wtheta', 'W0'
+values='Wchi', 'Wtheta', 'W0'
+
 [namelist:finite_element=coord_system]
 compulsory=true
 description=The coordinate system that will be stored in the chi coordinate

components/driver/rose-meta/lfric-driver/versions.py

@@ -31,3 +31,23 @@ def upgrade(self, config, meta_config=None):
         # Add settings
         return config, self.reports
 """
+
+class vn30_t238(MacroUpgrade):
+    """Upgrade macro for PR #238 by Thomas Bendall."""
+
+    BEFORE_TAG = "vn3.0"
+    AFTER_TAG = "vn3.0_t238"
+
+    def upgrade(self, config, meta_config=None):
+        # Commands From: rose-meta/lfric-driver
+        self.add_setting(
+            config, ["namelist:finite_element", "coord_space"], "'Wchi'"
+        )
+        coord_order = self.get_setting_value(
+            config, ["namelist:finite_element", "coord_order"]
+        )
+        self.add_setting(
+            config, ["namelist:finite_element", "coord_order_multigrid"],
+            coord_order
+        )
+        return config, self.reports

components/driver/source/driver_fem_mod.f90

@@ -12,11 +12,24 @@
 !>           * Initialises function space chains for use by the model.
 module driver_fem_mod
 
+  use base_mesh_config_mod,           only: prime_mesh_name,                   &
+                                            geometry,                          &
+                                            geometry_spherical,                &
+                                            geometry_planar,                   &
+                                            topology,                          &
+                                            topology_non_periodic
   use sci_chi_transform_mod,          only: init_chi_transforms, &
                                             final_chi_transforms
   use constants_mod,                  only: i_def, l_def, str_def
   use extrusion_mod,                  only: TWOD, PRIME_EXTRUSION
-  use finite_element_config_mod,      only: coord_order
+  use finite_element_config_mod,      only: coord_order,                       &
+                                            coord_order_multigrid,             &
+                                            coord_system,                      &
+                                            coord_system_xyz,                  &
+                                            coord_space,                       &
+                                            coord_space_W0,                    &
+                                            coord_space_Wchi,                  &
+                                            coord_space_Wtheta
   use field_mod,                      only: field_type
   use fs_continuity_mod,              only: W0, W2, W3, Wtheta, Wchi, W2v, W2h
   use function_space_mod,             only: function_space_type
@@ -68,6 +81,10 @@ subroutine init_fem( mesh_collection, chi_inventory, panel_id_inventory )
     type(field_type)                   :: panel_id
     type(function_space_type), pointer :: fs => null()
     integer(kind=i_def)                :: chi_space, coord, i
+    integer(kind=i_def)                :: coord_order_h, coord_order_v
+    integer(kind=i_def)                :: this_coord_order
+    integer(kind=i_def)                :: halo_depth
+    logical(kind=l_def)                :: is_valid
 
     character(str_def) :: mesh_name
 
@@ -98,27 +115,78 @@ subroutine init_fem( mesh_collection, chi_inventory, panel_id_inventory )
       ! Only create coordinates for 3D meshes
       if (mesh%get_extrusion_id() /= TWOD) then
 
-        ! Initialise panel ID field object ---------------------------------------
+        ! Initialise panel ID field object -------------------------------------
         twod_mesh => mesh_collection%get_mesh(mesh, TWOD)
         fs => function_space_collection%get_fs(twod_mesh, 0, 0, W3)
-        call panel_id%initialise( vector_space = fs, halo_depth = twod_mesh%get_halo_depth() )
+        halo_depth = twod_mesh%get_halo_depth()
+        call panel_id%initialise(fs, halo_depth=halo_depth)
 
-        ! Initialise chi field object --------------------------------------------
-        if ( coord_order == 0 ) then
-          chi_space = W0
-          write(log_scratch_space,'(A)') &
-              'Computing W0 coordinate fields for ' // trim(mesh_name) // 'mesh'
-          call log_event( log_scratch_space, log_level_info )
+        ! Initialise chi field object ------------------------------------------
+        ! Set coordinate order for this mesh
+        if (all_mesh_names(i) == prime_mesh_name) then
+          this_coord_order = coord_order
         else
+          this_coord_order = coord_order_multigrid
+        end if
+
+        ! Determine coordinate space
+        select case (coord_space)
+        case (coord_space_W0)
+          ! Check domain/topology is valid
+          is_valid = (                                                         &
+              geometry == geometry_spherical                                   &
+              .and. coord_system == coord_system_xyz                           &
+            ) .or. (                                                           &
+              geometry == geometry_planar                                      &
+              .and. topology == topology_non_periodic                          &
+          )
+          if (.not. is_valid) then
+            call log_event(                                                    &
+                    'Coordinate space W0 is only valid for non-periodic ' //   &
+                    'planar domains or when using the xyz coordinate system.', &
+                    LOG_LEVEL_ERROR                                            &
+            )
+          end if
+
+          ! Correct the coord_order for W0 polynomials being 1 order above W3
+          this_coord_order = this_coord_order - 1
+          chi_space = W0
+
+        case (coord_space_Wchi)
           chi_space = Wchi
-          write(log_scratch_space,'(A)') &
-              'Computing Wchi coordinate fields for ' // trim(mesh_name) // 'mesh'
-          call log_event( log_scratch_space, log_level_info )
+
+        case (coord_space_Wtheta)
+          chi_space = Wtheta
+
+        case default
+          call log_event('Invalid value for coord_space', LOG_LEVEL_ERROR)
+        end select
+
+        ! Set horizontal and vertical coordinate orders separately
+        if (coord_system == coord_system_xyz) then
+          ! Geocentric Cartesian coordinates - same order in all directions
+          coord_order_h = this_coord_order
+          coord_order_v = this_coord_order
+
+        else
+          ! For native coordinates, we separate horizontal and vertical coords
+          ! and can still accurately represent space with linear vertical coords
+          coord_order_h = this_coord_order
+
+          if (coord_space == coord_space_Wchi) then
+            coord_order_v = 1  ! Linear vertical coords for Wchi
+          else
+            coord_order_v = 0  ! Linear vertical coords for Wtheta
+          end if
         end if
-        fs => function_space_collection%get_fs(mesh, coord_order, coord_order, chi_space)
+
+        ! Create coordinate space
+        fs => function_space_collection%get_fs(                                &
+            mesh, coord_order_h, coord_order_v, chi_space                      &
+        )
 
         do coord = 1, size(chi)
-          call chi(coord)%initialise(vector_space = fs, halo_depth = twod_mesh%get_halo_depth() )
+          call chi(coord)%initialise(fs, halo_depth=halo_depth)
         end do
 
         ! Set coordinate fields --------------------------------------------------

components/driver/unit-test/assign_coordinate_xyz_mod_test.pf

@@ -34,7 +34,7 @@ contains
     use feign_config_mod,          only : feign_base_mesh_config, &
                                           feign_finite_element_config
     use finite_element_config_mod, only : cellshape_quadrilateral, &
-                                          coord_system_xyz
+                                          coord_system_xyz, coord_space_wchi
 
     implicit none
 
@@ -49,6 +49,8 @@ contains
 
     call feign_finite_element_config( cellshape=cellshape_quadrilateral, &
                                       coord_order=0_i_def,               &
+                                      coord_order_multigrid=1_i_def,     &
+                                      coord_space=coord_space_wchi,      &
                                       coord_system=coord_system_xyz,     &
                                       element_order_h=0_i_def,           &
                                       element_order_v=0_i_def,           &

components/lfric-xios/integration-test/configs/context.nml

@@ -21,6 +21,8 @@
   element_order_v = 0
   rehabilitate    =.true.
   coord_order     = 1
+  coord_order_multigrid = 1
+  coord_space     = 'Wchi'
   coord_system    ='native'
 /
 

components/lfric-xios/integration-test/configs/non_cyclic_base.nml

@@ -21,6 +21,8 @@
   element_order_v = 0
   rehabilitate    =.true.
   coord_order     = 1
+  coord_order_multigrid = 1
+  coord_space     = 'Wchi'
   coord_system    ='native'
 /
 

components/lfric-xios/integration-test/resources/configs/context.nml

@@ -21,6 +21,8 @@
   element_order_v = 0
   rehabilitate    =.true.
   coord_order     = 1
+  coord_order_multigrid = 1
+  coord_space     = 'Wchi'
   coord_system    ='native'
 /
 

components/lfric-xios/integration-test/resources/configs/non_cyclic_base.nml

@@ -21,6 +21,8 @@
   element_order_v = 0
   rehabilitate    =.true.
   coord_order     = 1
+  coord_order_multigrid = 1
+  coord_space     = 'Wchi'
   coord_system    ='native'
 /
 

components/lfric-xios/source/lfric_xios_setup_mod.x90

@@ -310,8 +310,9 @@ contains
     coord_dim_owned = proxy_coord_output(1)%vspace%get_last_dof_owned() / nfull_levels
 
     ! Obtain sample_chi, which will be used for setting up XIOS coordinates
-    if ( mesh%is_geometry_spherical() ) then
+    if ( mesh%is_geometry_spherical() .or. chi(1)%which_function_space() == Wtheta ) then
       ! Sample chi on W0 function space to prevent "unzipping" of cubed-sphere mesh
+      ! This also avoids issues when Wtheta is used as the coordinate space
       do i = 1,3
         call sample_chi(i)%initialise( vector_space = output_field_fs )
       end do
@@ -322,7 +323,10 @@ contains
       do i = 1,3
         call chi(i)%copy_field_serial(sample_chi(i))
       end do
-      ! Tell downstream output processing that the mesh is planar
+    end if
+
+    ! Tell downstream output processing that the mesh is planar
+    if ( .not. mesh%is_geometry_spherical() ) then
       call set_xios_geometry_planar()
     end if
 

components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w2_mod.F90

@@ -15,11 +15,11 @@ module sci_compute_mass_matrix_kernel_w2_mod
                                      GH_READ, GH_WRITE,         &
                                      GH_REAL, ANY_W2,           &
                                      ANY_DISCONTINUOUS_SPACE_3, &
+                                     ANY_SPACE_9,               &
                                      GH_BASIS, GH_DIFF_BASIS,   &
                                      CELL_COLUMN, GH_QUADRATURE_XYoZ
   use constants_mod,           only: i_def, r_def
   use sci_coordinate_jacobian_mod, only: coordinate_jacobian
-  use fs_continuity_mod,       only: Wchi
   use kernel_mod,              only: kernel_type
 
   implicit none
@@ -34,12 +34,12 @@ module sci_compute_mass_matrix_kernel_w2_mod
     private
     type(arg_type) :: meta_args(3) = (/                                      &
          arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, ANY_W2, ANY_W2),           &
-         arg_type(GH_FIELD*3,  GH_REAL, GH_READ,  Wchi),                     &
+         arg_type(GH_FIELD*3,  GH_REAL, GH_READ,  ANY_SPACE_9),              &
          arg_type(GH_FIELD,    GH_REAL, GH_READ,  ANY_DISCONTINUOUS_SPACE_3) &
          /)
     type(func_type) :: meta_funcs(2) = (/                                    &
-         func_type(ANY_W2, GH_BASIS),                                        &
-         func_type(Wchi,   GH_BASIS, GH_DIFF_BASIS)                          &
+         func_type(ANY_W2,      GH_BASIS),                                   &
+         func_type(ANY_SPACE_9, GH_BASIS, GH_DIFF_BASIS)                     &
          /)
     integer :: operates_on = CELL_COLUMN
     integer :: gh_shape = GH_QUADRATURE_XYoZ