Pass read-only arrays of scalars to kernels

[TeranIvy]

We have had repeated requests for passing arrays of scalars to kernels. The design is outlined below.

Note: The code in the examples was successfully tested with Intel 19.0.1 and Gnu 11.2.0 compilers.

Kernel layer

Full example: scalar_arrays_kernel_mod.f90.txt

Metadata and argument access

• As for scalars, only read-only user-defined scalar arrays are allowed.
• Keyword to denote the scalar arrays is ~~GH_ARRAY~~ GH_SCALAR_ARRAY. We also need to specify the number of ranks as an integer. ~~We broadly agreed on NRANKS as a keyword for ranks (rank is a Fortran intrinsic so cannot be used here).~~ The new keyword for scalar array will need to be added to argument_mod.F90 where LFRic stores metadata IDs.
• An example of ~~suggested~~ metadata design updated after the metadata review that compiles is as follows:

     type(arg_type), dimension(3) :: meta_args = (/          &
          arg_type(GH_SCALAR_ARRAY, GH_REAL,    GH_READ, 2), &
          arg_type(GH_SCALAR_ARRAY, GH_LOGICAL, GH_READ, 1), &
          arg_type(GH_SCALAR_ARRAY, GH_INTEGER, GH_READ, 4)  &
          /)

Note: We are not allowing kernels with no fields so this example is just for illustration.

Arguments

• For every scalar array, pass:
  • ~~An integer scalar of kind i_def for the number of ranks, nranks_<array_name>;~~ This is not required as it is read from the metadata.
  • A rank-1 integer array of kind i_def and size nranks_<array_name> containing the upper bounds for each rank, dims_<array_name> (lower bound is assumed to be 1 as this is how Fortran passes array slices to subroutines by default).
• Pass array of data type and kind as inferred from the metadata.

The argument list for the metadata above would look something like:

  subroutine scalar_arrays_code(dims_rarray, real_array,    &
                                dims_larray, logical_array, &
                                dims_iarray, integer_array)

and the declarations

    integer(i_def), intent(in), dimension(2) :: dims_rarray
    integer(i_def), intent(in), dimension(1) :: dims_larray
    integer(i_def), intent(in), dimension(4) :: dims_iarray
    real(r_def), intent(in), dimension(dims_rarray(1),dims_rarray(2)) :: real_array
    logical(l_def), intent(in), dimension(dims_larray(1)) :: logical_array
    integer(i_def), intent(in), dimension(dims_iarray(1),dims_iarray(2),dims_iarray(3),dims_iarray(4)) :: integer_array

PSy layer

Full example: pass_scalar_array_mod_psy.f90.txt

Note: As the array rank is known from the metadata, we can simply do

    integer(i_def), parameter :: nranks_rarray1 = 2
    integer(i_def), parameter :: nranks_larray2 = 1
    integer(i_def), parameter :: nranks_iarray3 = 4
    integer(i_def), dimension(nranks_rarray1) :: dims_rarray1
    integer(i_def), dimension(nranks_larray2) :: dims_larray2
    integer(i_def), dimension(nranks_iarray3) :: dims_iarray3
    ...
    ! Get the upper bound for each rank of each scalar array
    dims_rarray1 = shape(rarray1)
    dims_larray2 = shape(larray2)
    dims_iarray3 = shape(iarray3)

Deprecated ~~We can declare scalar arrays as assumed shape (the rank is known from the metadata), e.g.~~

    real(r_def), dimension(:,:) :: rarray
    logical(l_def), dimension(:) :: larray2
    integer(i_def), dimension(:,:,:,:) :: iarray3

~~The number of ranks and dimensions can then be retrieved using rank and size intrinsics, e.g.~~

    integer(i_def), parameter :: nranks_rarray1 = rank(rarray1)
    integer(i_def), parameter :: nranks_larray2 = rank(larray2)
    integer(i_def), parameter :: nranks_iarray3 = rank(iarray3)
    integer(i_def), dimension(nranks_rarray1) :: dims_rarray1
    integer(i_def), dimension(nranks_larray2) :: dims_larray2
    integer(i_def), dimension(nranks_iarray3) :: dims_iarray3
    ...
    ! Get the upper bound for each rank of each scalar array
    dims_rarray1 = shape(rarray1)
    dims_larray2 = shape(larray2)
    dims_iarray3 = shape(iarray3)

Algorithm layer

Full example: pass_scalar_array_alg.f90.txt

An invoke call would simply pass the name of the scalar arrays, as for other arguments, e.g.

  call invoke_0(real_array, logical_array, integer_array)

[andrewcoughtrie]

Example PSyKAl-lite code using the above implementation can be found in um_physics/source/psy/psykal_lite_phys_mod.F90 on the LFRic trunk, this is used as part of the Stochastic physics forcing pattern.

[arporter]

Priority wise, how does this compare with the 'stencils for domain (i-first) kernels' issue @andrewcoughtrie and @TeranIvy ?

[TeranIvy]

Link to PSyKAl-lite code..

@arporter, stencils are currently a priority because of optimisations that the relevant code introduces.

[andrewcoughtrie]

As far as i know we don't have an urgent need for this to be implemented, the current PSyKAl-lite code is only in one place that I know of and is unlikely to expand.

[arporter]

Thanks both :-)

[MatthewHambley]

I have discovered a case of an alternative workaround whereby global scope module variables are being used to pass arrays "around the side" of the argument list. Clearly this is bad. For those interested check out Gung Ho's init_thermo_profile_alg_mod.

[andrewcoughtrie]

Can you put a comment next to this usage with the PR reference number please @MatthewHambley, just to make it easy to search for places it needs changing when we have the functionality.

[MatthewHambley]

Can you put a comment next to this usage with the PR reference number please @MatthewHambley, just to make it easy to search for places it needs changing when we have the functionality.

It is already marked with a Doxygen "todo" comment. That's how I knew this PR existed.

[rupertford]

Thanks @MatthewHambley. It looks like the case you are talking about is in the algorithm layer. Is that right?

[mo-lottieturner]

Changing naming to distinguish NRANKS*n (previously called array_nranks, now array_size) from n extracted from said string (previously also called array_nranks, sometimes array_ranks, now called array_ndims)