Clarify rules for determining nlayers in LFRic kernels

[TeranIvy]

LFRic API seems to return different nlayers depending on whether invoke calls to setval_c built-in and a kernel are joint or not. This was commented on in the LFRic Apps ticket 21, pres_lev_diags_alg_mod.x90. The example pruned-down code is below.

Mock algorithm code

program test_invokes

  ! Description: test single and joint invokes
  use constants_mod, only: r_def
  use field_mod,     only: field_type
  use testkern_mod,  only: testkern_type

  implicit none

  type(field_type) :: plev_heaviside, exner_wth

  ! Case 1: Separate invokes, setval_c before the kernel call
  call invoke( setval_c(plev_heaviside, 1.0_r_def) )
  call invoke( &
       testkern_type(exner_wth, plev_heaviside) )

  ! Case 2: Separate invokes, setval_c after the kernel call
  call invoke( &
       testkern_type(exner_wth, plev_heaviside) )
  call invoke( setval_c(plev_heaviside, 1.0_r_def) )

  ! Case 3: Joint invoke, setval_c before the kernel call
  call invoke( setval_c(plev_heaviside, 1.0_r_def), &
               testkern_type(exner_wth, plev_heaviside) )

  ! Case 3: Joint invoke, setval_c after the kernel call
  call invoke( testkern_type(exner_wth, plev_heaviside), &
               setval_c(plev_heaviside, 1.0_r_def) )

end program test_invokes

Mock kernel code

module testkern_mod

  use argument_mod
  use fs_continuity_mod
  use kernel_mod
  use constants_mod

  implicit none

  type, extends(kernel_type) :: testkern_type
     type(arg_type), dimension(2) :: meta_args =                                   &
          (/ arg_type(GH_FIELD, GH_REAL, GH_READ,      ANY_DISCONTINUOUS_SPACE_1), &
             arg_type(GH_FIELD, GH_REAL, GH_READWRITE, ANY_DISCONTINUOUS_SPACE_2)  &
           /)
     integer :: operates_on = cell_column
   contains
     procedure, nopass :: code => testkern_code
  end type testkern_type

contains

  subroutine testkern_code(nlayers,                &
                          exner,                   &
                          plev_heaviside,          &
                          ndf_in, undf_in, map_in, &
                          ndf_out, undf_out, map_out)

    implicit none

    integer(kind=i_def), intent(in) :: nlayers
    integer(kind=i_def), intent(in) :: ndf_in
    integer(kind=i_def), intent(in), dimension(ndf_in) :: map_in
    integer(kind=i_def), intent(in) :: ndf_out
    integer(kind=i_def), intent(in), dimension(ndf_out) :: map_out
    integer(kind=i_def), intent(in) :: undf_in, undf_out
    real(kind=r_def), intent(in), dimension(undf_in) :: exner
    real(kind=r_def), intent(inout), dimension(undf_out) :: plev_heaviside

  end subroutine testkern_code

end module testkern_mod

Generated PSy-layer code for nlayers

Note that setval_c built-in defines the field access as GH_WRITE (see https://github.com/stfc/PSyclone/blob/master/src/psyclone/parse/lfric_builtins_mod.f90).

When invokes are separate (Case 1 and Case 2), nlayers is determined from the first field in the kernel argument list, exner, which has GH_READ access: nlayers = exner_wth_proxy%vspace%get_nlayers().

When invokes are joint, nlayers is determined from the first field in the first kernel or built-in call in the joint invoke. In Case 3 this is nlayers = plev_heaviside_proxy%vspace%get_nlayers(), and in Case 4 this is nlayers = exner_wth_proxy%vspace%get_nlayers().

Questions/discussion for @arporter, @sergisiso and @hiker

• Is this intended behaviour and is it documented? I grepped for nlayers in both user and developer guide and I was not able to find a relevant record.
• Would it make sense that the nlayers is determined from fields that are updated, regardless of the order of metadata? This is already done for loop limits over cells (with continuous writers having priority regardless of the metadata order).
• There is issue #868 for custom function spaces to pass different number of layers. Since this issue is quite encompassing, would there be a way in the meantime to signal PSyclone what field it would extract the desired number of layers from?