(Closes #868) implement support for kernels with field arguments that specify nlayers/ndata

doc/tutorials_and_examples/examples_intro.rst

@@ -105,14 +105,14 @@ so it can be recorded in this table.
 
 .. tabularcolumns:: |l|L|
 
-=======================   =======================================================
+========================  =======================================================
 Compiler                  Version
-=======================   =======================================================
+========================  =======================================================
 Gnu Fortran               13.3, 14.3, 15.2, 16.1
 Intel Fortran (ifort)     21
 Intel LLVM Fortran (ifx)  2025.3.1, 2026.1.0
 NVIDIA Fortran            23.5
-=======================   =======================================================
+========================  =======================================================
 
 .. _examples_dependencies:
 

doc/user_guide/lfric.rst

@@ -1744,6 +1744,8 @@ meshes cannot be on the same function space while those on the same
 mesh must also be on the same function space.
 
 
+.. _number-of-layers-metadata:
+
 Number of Layers Metadata
 """""""""""""""""""""""""
 
@@ -1774,6 +1776,7 @@ If two or more field/operator arguments are on the same function space
 and have the same number of layers (whether a literal or a name) then
 only one dofmap is passed to the kernel for those arguments.
 
+.. _ndata-metadata:
 
 Multi-Data Metadata
 """""""""""""""""""
@@ -2090,11 +2093,23 @@ conventions, are:
 
 1) If an LMA operator is passed then include the ``cells`` argument.
    ``cells`` is an ``integer`` of kind ``i_def`` and has intent ``in``.
-2) Include ``nlayers``, the number of layers in a column. ``nlayers``
-   is an ``integer`` of kind ``i_def`` and has intent ``in``. PSyclone
-   will obtain the value of ``nlayers`` to use for a particular kernel
-   from the first field or operator in the argument list.
-3) For each scalar/field/vector_field/operator/ScalarArray in the order specified by
+
+2) Include each unique value of ``nlayers``, the number of layers in a
+   column, for all of the field/operator arguments. ``nlayers``
+   is an ``integer`` of kind ``i_def`` and has intent ``in``. A typical
+   kernel will have just one value of ``nlayers`` which PSyclone will
+   obtain from the first field/operator argument in the list. For those
+   kernels which have field arguments with ``nlayers`` specified by a
+   label in the metadata (see :ref:`number-of-layers-metadata`), include
+   a value of ``nlayers`` for every unique label, in the same
+   order as they occur in the ``meta_args`` list.
+
+3) If any field/operator arguments have custom values of ``ndata`` (see
+   :ref:`ndata-metadata`) then, for each unique label, include an ``ndata``
+   argument in the same order as they appear in the ``meta_args`` list. Each
+   ``ndata`` is an ``integer`` of kind ``i_def`` and has intent ``in``.
+
+4) For each scalar/field/vector_field/operator/ScalarArray in the order specified by
    the meta_args metadata:
 
    1) If the current entry is a scalar quantity then include the Fortran
@@ -2131,16 +2146,7 @@ conventions, are:
       4) If the field entry stencil access is of type ``XORY1D`` then
          add an additional ``integer`` direction argument of kind
          ``i_def`` and with intent ``in``.
-      5) If the field is multi-data and the value of ``NDATA`` is unknown
-         (i.e. it is specified in the metadata using a label)
-         then add an additional ``integer``, scalar argument of kind ``i_def``
-         and intent ``in``. If the ``NDATA`` value is common to more than one
-         kernel argument, it is only added for the first such argument.
-      6) If the field has an unknown (i.e. specified in the metadata using a
-         label), custom number of vertical levels then pass this as an
-         additional ``integer``, scalar argument of kind ``i_def`` and
-         intent ``in``. If the number of vertical levels is common to more than
-         one kernel argument, it is only added for the first such argument.
+
    3) If the current entry is a field vector then for each dimension
       of the vector, include a field array. The field array name is
       specified as
@@ -2166,7 +2172,7 @@ conventions, are:
       the data type and kind specified in the metadata. The ScalarArray
       must be denoted with intent ``in`` to match its read-only nature.
 
-4) DoF maps for function spaces are handled in the order they appear in the
+5) DoF maps for function spaces are handled in the order they appear in the
    metadata arguments (the ``to`` function space of an operator is considered
    to be before the ``from`` function space of the same operator as it appears
    first in lexicographic order). Note that if two fields on a given function
@@ -2247,7 +2253,7 @@ conventions, are:
       |               |           | W2Vtrace, W3, Wtheta, Wchi         |
       +---------------+-----------+------------------------------------+
 
-5) If either the ``normals_to_horizontal_faces`` or
+6) If either the ``normals_to_horizontal_faces`` or
    ``outward_normals_to_horizontal_faces`` properties of the reference
    element are required then pass the number of horizontal faces of the
    reference element (``nfaces_re_h``). Similarly, if either the
@@ -2269,7 +2275,7 @@ conventions, are:
       a rank-2 ``integer`` array of kind ``i_def`` with dimensions
       ``(3, nfaces_re)``.
 
-6) If the ``adjacent_face`` mesh property is required then:
+7) If the ``adjacent_face`` mesh property is required then:
 
    1) If the number of horizontal cell faces obtained from the reference
       element (``nfaces_re_h``) is not already being passed to the kernel (due
@@ -2278,7 +2284,7 @@ conventions, are:
    2) Pass a rank-1, ``integer`` array of kind ``i_def`` and extent
       ``nfaces_re_h``.
 
-7) If Quadrature is required (``gh_shape = gh_quadrature_*``) then, for
+8) If Quadrature is required (``gh_shape = gh_quadrature_*``) then, for
    each shape in the order specified in the ``gh_shape`` metadata:
 
    1) Include ``integer``, scalar arguments of kind ``i_def`` with intent

examples/common.mk

@@ -94,6 +94,12 @@ ${NOTEBOOK_FILES}:
 
 notebook: ${NOTEBOOK_FILES}
 
+%.o: %.F90
+	$(F90) $(F90FLAGS) -c $<
+
+%.o: %.f90
+	$(F90) $(F90FLAGS) -c $<
+
 # By default we clean-up emacs backup files and generated Jupyter notebooks
 clean:
 	${RM} ./*~ ./*.nbconvert.ipynb ${GENERATED_FILES}

examples/lfric/README.md

@@ -390,6 +390,13 @@ the algorithm subroutine generated by this example:
  5. Run ``../bin/skeleton ./configuration.nml`` in the ``miniapps/skeleton/example`` directory.
 
 
+## Example 21: Kernels with arguments that specify ndata and/or nlayers
+
+Shows the use of kernels that have field arguments that specify non-default
+values of ndata and/or nlayers. Although the code compiles, it cannot be
+executed because, for simplicity, it lacks the necessary set-up of the LFRic
+field objects.
+
 ## Code
 
 Location of LFRic algorithm and kernel code that is used by two or

examples/lfric/code/testkern_nlayers_ndata2_mod.F90

@@ -0,0 +1,69 @@
+! -----------------------------------------------------------------------------
+! SPDX-FileCopyrightText: Copyright (c) 2026 Science and Technology
+!                         Facilities Council
+! SPDX-License-Identifier: BSD-3-Clause
+! See the full LICENSE file in the project root for details.
+! -----------------------------------------------------------------------------
+
+!> An example LFRic kernel which has arguments with non-default values
+!! of NLAYERS and NDATA.
+module testkern_nlayers_ndata2_mod
+
+  use argument_mod
+  use fs_continuity_mod
+  use kernel_mod
+  use constants_mod
+
+  implicit none
+
+  type, extends(kernel_type) :: testkern_nlayers_ndata2_type
+     type(arg_type), dimension(5) :: meta_args = (/                           &
+             arg_type(gh_field,  gh_real, gh_inc,  w1),                       &
+             ! Non-default number of layers.
+             arg_type(gh_field,  gh_real, gh_read, w2, nlayers="shallow"),    &
+             ! Non-default number of data values per dof.
+             arg_type(gh_field,  gh_real, gh_read, w2, ndata="precip"),       &
+             ! Although we have seen the values of ndata and nlayers before
+             ! (individually), the dof map will be different.
+             arg_type(gh_field,  gh_real, gh_read, w2, ndata="precip",        &
+                      nlayers="shallow"),                                     &
+             arg_type(gh_field, gh_real, gh_read, w2, ndata="3", nlayers="1") &
+           /)
+     integer :: operates_on = cell_column
+   contains
+     procedure, nopass :: code => testkern_nlayers_ndata2_code
+  end type testkern_nlayers_ndata2_type
+
+contains
+
+  subroutine testkern_nlayers_ndata2_code(          &
+       nlayers, nlayers_shallow, ndata_precip,      &
+       fld1, fld2, fld3, fld4, fld5,                &
+       ndf_w1, undf_w1, map_w1,                     &
+       ndf_fld2, undf_fld2, map_w2_fld2,            &
+       ndf_fld3, undf_fld3, map_w2_fld3,            &
+       ndf_fld4, undf_fld4, map_w2_fld4,            &
+       ndf_fld5, undf_fld5, map_w2_fld5)
+    implicit none
+
+    integer(kind=i_def), intent(in) :: nlayers
+    integer(kind=i_def), intent(in) :: nlayers_shallow
+    integer(kind=i_def), intent(in) :: ndata_precip
+    integer(kind=i_def), intent(in) :: ndf_w1, ndf_fld2, ndf_fld3, ndf_fld4, &
+                                       ndf_fld5
+    integer(kind=i_def), intent(in) :: undf_w1, undf_fld2, undf_fld3, &
+                                       undf_fld4, undf_fld5
+    integer(kind=i_def), intent(in), dimension(ndf_w1)   :: map_w1
+    integer(kind=i_def), intent(in), dimension(ndf_fld2) :: map_w2_fld2
+    integer(kind=i_def), intent(in), dimension(ndf_fld3) :: map_w2_fld3
+    integer(kind=i_def), intent(in), dimension(ndf_fld4) :: map_w2_fld4
+    integer(kind=i_def), intent(in), dimension(ndf_fld5) :: map_w2_fld5
+    real(kind=r_def), intent(inout), dimension(undf_w1) :: fld1
+    real(kind=r_def), intent(in), dimension(undf_fld2)  :: fld2
+    real(kind=r_def), intent(in), dimension(undf_fld3)  :: fld3
+    real(kind=r_def), intent(in), dimension(undf_fld4)  :: fld4
+    real(kind=r_def), intent(in), dimension(undf_fld5)  :: fld5
+
+  end subroutine testkern_nlayers_ndata2_code
+
+end module testkern_nlayers_ndata2_mod

examples/lfric/code/testkern_nlayers_ndata_mod.F90

@@ -0,0 +1,75 @@
+! -----------------------------------------------------------------------------
+! SPDX-FileCopyrightText: Copyright (c) 2026 Science and Technology
+!                         Facilities Council
+! SPDX-License-Identifier: BSD-3-Clause
+! See the full LICENSE file in the project root for details.
+! -----------------------------------------------------------------------------
+
+!> An example LFRic kernel which has arguments with non-default values
+!! of NLAYERS and NDATA.
+module testkern_nlayers_ndata_mod
+
+  use argument_mod
+  use fs_continuity_mod
+  use kernel_mod
+  use constants_mod
+
+  implicit none
+
+  type, extends(kernel_type) :: testkern_nlayers_ndata_type
+     type(arg_type), dimension(7) :: meta_args =                           &
+          (/ arg_type(gh_scalar, gh_real, gh_read),                        &
+             arg_type(gh_field,  gh_real, gh_inc,  w1),                    &
+             arg_type(gh_field,  gh_real, gh_read, w2),                    &
+             ! Non-default number of layers.
+             arg_type(gh_field,  gh_real, gh_read, w2, nlayers="shallow"), &
+             ! Non-default number of layers but same as previous arg. so
+             ! has same dof map.
+             arg_type(gh_field,  gh_real, gh_read, w2, nlayers="shallow"), &
+             ! Non-default number of data values per dof.
+             arg_type(gh_field,  gh_real, gh_read, w2, ndata="precip"),    &
+             ! Although we have seen the values of ndata and nlayers before
+             ! (individually), the dof map will be different.
+             arg_type(gh_field,  gh_real, gh_read, w2, ndata="precip",     &
+                      nlayers="shallow")                                   &
+           /)
+     integer :: operates_on = cell_column
+   contains
+     procedure, nopass :: code => testkern_nlayers_ndata_code
+  end type testkern_nlayers_ndata_type
+
+contains
+
+  subroutine testkern_nlayers_ndata_code(           &
+       nlayers, nlayers_shallow, ndata_precip,      &
+       ascalar, fld1, fld2, fld3, fld4, fld5, fld6, &
+       ndf_w1, undf_w1, map_w1,                     &
+       ndf_w2, undf_w2, map_w2,                     &
+       ndf_fld3, undf_fld3, map_w2_fld3,            &
+       ndf_fld5, undf_fld5, map_w2_fld5,            &
+       ndf_fld6, undf_fld6, map_w2_fld6)
+    implicit none
+
+    integer(kind=i_def), intent(in) :: nlayers
+    integer(kind=i_def), intent(in) :: nlayers_shallow
+    integer(kind=i_def), intent(in) :: ndata_precip
+    integer(kind=i_def), intent(in) :: ndf_w1, ndf_w2, ndf_fld3, ndf_fld5, &
+                                       ndf_fld6
+    integer(kind=i_def), intent(in) :: undf_w1, undf_w2, undf_fld3, &
+                                       undf_fld5, undf_fld6
+    integer(kind=i_def), intent(in), dimension(ndf_w1)   :: map_w1
+    integer(kind=i_def), intent(in), dimension(ndf_w2)   :: map_w2
+    integer(kind=i_def), intent(in), dimension(ndf_fld3) :: map_w2_fld3
+    integer(kind=i_def), intent(in), dimension(ndf_fld5) :: map_w2_fld5
+    integer(kind=i_def), intent(in), dimension(ndf_fld6) :: map_w2_fld6
+    real(kind=r_def), intent(in) :: ascalar
+    real(kind=r_def), intent(inout), dimension(undf_w1) :: fld1
+    real(kind=r_def), intent(in), dimension(undf_w2)    :: fld2
+    real(kind=r_def), intent(in), dimension(undf_fld3)  :: fld3
+    real(kind=r_def), intent(in), dimension(undf_fld3)  :: fld4
+    real(kind=r_def), intent(in), dimension(undf_fld5)  :: fld5
+    real(kind=r_def), intent(in), dimension(undf_fld6)  :: fld6
+
+  end subroutine testkern_nlayers_ndata_code
+
+end module testkern_nlayers_ndata_mod

examples/lfric/eg21/Makefile

@@ -0,0 +1,42 @@
+# -----------------------------------------------------------------------------
+# SPDX-FileCopyrightText: Copyright (c) 2026 Science and Technology
+#                         Facilities Council
+# SPDX-License-Identifier: BSD-3-Clause
+# See the full LICENSE file in the project root for details.
+# -----------------------------------------------------------------------------
+
+include ../lfric_common.mk
+
+ALG_NAME = nlayers_ndata_eg_alg
+
+OBJ = testkern_nlayers_ndata_mod.o testkern_nlayers_ndata2_mod.o ${ALG_NAME}_psy.o ${ALG_NAME}.o
+GENERATED_FILES = ${ALG_NAME}.f90 ${ALG_NAME}_psy.f90 \
+    testkern_nlayers_ndata_mod.F90 testkern_nlayers_ndata2_mod.F90 *.o \
+    *.mod will_not_run.exe
+
+# Although we link an executable, the example is too simplistic to
+# actually execute.
+run: compile
+	@echo "lfric/eg21 is not executable."
+
+# Dependencies
+%.F90:
+	ln -s ../code/$@ .
+
+# Keep the generated psy and alg files
+.precious: ${ALG_NAME}_psy.f90 ${ALG_NAME}.f90
+
+transform: ${ALG_NAME}_psy.f90
+
+# Process the algorithm layer with PSyclone to generate a compilable
+# algorithm file and a suitable PSy-layer subroutine.
+${ALG_NAME}_psy.f90 ${ALG_NAME}.f90: ${ALG_NAME}.x90
+	${PSYCLONE} -api lfric -d ../code -opsy ${ALG_NAME}_psy.f90 -oalg ${ALG_NAME}.f90 $<
+
+# Compile and link purely to test correctness - the resulting executable
+# won't actually work.
+compile: transform ${LFRIC_LIB} ${OBJ}
+	${F90} ${F90FLAGS} -o will_not_run.exe ${OBJ} ${LDFLAGS}
+
+allclean: clean
+	make -C $(LFRIC_PATH) allclean

examples/lfric/eg21/nlayers_ndata_eg_alg.x90

@@ -0,0 +1,25 @@
+! -----------------------------------------------------------------------------
+! SPDX-FileCopyrightText: Copyright (c) 2026 Science and Technology
+!                         Facilities Council
+! SPDX-License-Identifier: BSD-3-Clause
+! See the full LICENSE file in the project root for details.
+! -----------------------------------------------------------------------------
+
+!> A (non-executable) example LFRic algorithm layer that invokes two different
+!! kernels that both have field arguments with custom values of NLAYERS and NDATA.
+!! The example is not executable because it omits all the initialisation
+!! required to construct fully valid field objects.
+program my_alg
+  use constants_mod, only: r_def
+  use field_mod, only: field_type
+  use testkern_nlayers_ndata2_mod, only: testkern_nlayers_ndata2_type
+  use testkern_nlayers_ndata_mod, only: testkern_nlayers_ndata_type
+  implicit none
+  real(kind=r_def) :: scalar
+  type(field_type) :: fld1, fld2, fld3, fld4, fld5, fld6, fld_2d
+
+  call invoke( &
+     testkern_nlayers_ndata_type(scalar, fld1, fld2, fld3, fld4, fld5, fld6), &
+     testkern_nlayers_ndata2_type(fld1, fld2, fld3, fld4, fld_2d) )
+
+end program

external/lfric_infrastructure/src/kernel_metadata/argument_mod.f90

@@ -231,10 +231,10 @@ module argument_mod
      !> Optional metadata (fields only) specifying how many data values
      !! there are at each DoF. Defaults to 1.
      character(len=METADATA_CHAR_LENGTH) :: ndata = "1"
-     !> Optional metadata specifying how many vertical levels
+     !> Optional metadata specifying how many vertical layers
      !! the field has. Default is to use the number associated with the
      !! first field/operator kernel argument.
-     character(len=METADATA_CHAR_LENGTH) :: nlevels = ""
+     character(len=METADATA_CHAR_LENGTH) :: nlayers = ""
   end type arg_type
 
   !> Optional metadata for the basis/differential basis functions required with