Fortran bmi grid

[hellkite500]

Add gridded input and output variables to the fortran BMI test model. Unit tests covering at least two different 3D indexing orderings are included, and the following helps illustrate how gridded data is flattened and reconstructed across the Fortran/C BMI interface.

Multi Dimension Data between Fortran and C/C++ interfaces

In Fortran, the reshape intrinsic function can be used to "flatten" the array using the native indexing of the fortran grid.
For example, the following Fortran snippet will create a column major flattened array that can move across the BMI interface and be reconstructed by the framework. This creates a 2-Z dimension, 2-Y dimension, and 3-X dimension array, and then
reshapes it to a flattened array using the default (natural) order of the original array (z,y,x).

    double precision, dimension(2,2,3) :: data ! z, y, x
    double precision, dimension(12) :: flattened_data
    flattened_data = reshape(data, [size(data)])

In the example below, this flattened_data would look like [ 11, 101, 14, 104, 12, 102, 15, 105, 13, 103, 16, 106 ].

Pitfalls to be aware of in Fortran BMI

Beware that using the size intrinsic in certain BMI interface functions is problematic, since the some functions use size as a dummy variable, the dummy variable shadows the intrinsic name and size will refer to that variable in the scope of the function.

This can present as a puzzling bug when first encounted. Consider the following implementation of a simple loop for grid size

function test_grid_size(this, grid, size) result (bmi_status)
    class (bmi_test_bmi), intent(in) :: this
    integer, intent(in) :: grid
    integer, intent(out) :: size
    integer :: bmi_status
    
    integer :: i
    do i = 1, SIZE(grids)
      if ( grids(i)%id .eq. grid ) then
        size = product( grids(i)%shape )
        bmi_status = BMI_SUCCESS
        return
      end if
    end do

  end function test_grid_size

This produces compiler errors that confusing without this intrinsic aliasing context

  500 |         size = product( grids(i)%shape )
      |             1
Error: 'size' at (1) is not a variable
bmi_test_bmi_fortran.f90:480:42:

  480 |   function test_grid_size(this, grid, size) result (bmi_status)
      |                                          1
Error: PROCEDURE attribute conflicts with INTENT attribute in 'size' at (1)
bmi_test_bmi_fortran.f90:1008:41:

 1008 |       bmi_status = this%get_grid_size(1, size)
      |                                         1

There are some workarounds using block ... end block that might get around this in certain scenarios. For the grid size function above, it can be written as

function test_grid_size(this, grid, size) result (bmi_status)
    class (bmi_test_bmi), intent(in) :: this
    integer, intent(in) :: grid
    integer, intent(out) :: size
    integer :: bmi_status
    
    integer :: i, upper
    block
      intrinsic :: size
      upper = size(grids)
    end block
    do i = 1, upper
      if ( grids(i)%id .eq. grid ) then
        size = product( grids(i)%shape )
        bmi_status = BMI_SUCCESS
        return
      end if
    end do

But there are other cases where it isn't as straight forward. For example, if one needs to get the shape of a dummy variable named shape, then one must copy the dummy variable first, then apply the block trick, e.g.

integer, dimension(:) :: shape
integer, dimension(:) :: shape_
shape_ = shape
block
    intrinsic shape
    write(0,*) "Shape of shape is: ", shape(shape_)
end block

As of this writing, size, shape, and rank are all identified as victims of this shadowing in BMI 2.0.

Data Description

The test model applies this 2D grid to the 3D grid, with +10 in the Z0 dimension
and + 100 in the Z1 dimension, so the resulting grid should look like

Z0 =      x0   x1  x2
     y0: [11   12  13]
     y1: [14   15  16]

Z1 =      x0   x1  x2
     y0: [101   102  103]
     y1: [104   105  106]

Row Major Flattening given different indexing strategies

A row-major flattening of this looks like (last index changes fastest):

• assuming the index is (z,y,x):

  [ 11, 12, 13, 14, 15, 16, 101, 102, 103, 104, 105, 106 ]
  

• assuming the index is (z,x,y):

  [ 11, 14, 12, 15, 13, 16, 101, 104, 102, 105, 103, 106 ]
  

• assuming the index is (y,x,z):

  [ 11, 101, 12, 102, 13, 103, 14, 104, 15, 105, 16, 106 ]
  

• assuming the index is (y,z,x):

  [ 11, 12, 13, 101, 102, 103, 14, 15, 16, 104, 105, 106 ]
  

• assuming the index is (x,y,z):

  [ 11, 101, 14, 104, 12, 102, 15, 105, 13, 103, 16, 106 ]
  

• assuming the index is (x,z,y):

  [ 11, 14, 101, 104, 12, 15, 102, 105, 13, 16, 103, 106 ]
  

Column Major Flattening given different indexing strategies

A column-major flattening of this looks like (first index changes fastest):

• assuming the index is (z,y,x):

  [ 11, 101, 14, 104, 12, 102, 15, 105, 13, 103, 16, 106 ]
  

• assuming the index is (z,x,y):

  [ 11, 101, 12, 102, 13, 103, 14, 104, 15, 105, 16, 106 ]
  

• assuming the index is (y,x,z):

  [ 11, 14, 12, 15, 13, 16, 101, 104, 102, 105, 103, 106 ]
  

• assuming the index is (y,z,x):

  [ 11, 14, 101, 104, 12, 15, 102, 105, 13, 16, 103, 106 ]
  

• assuming the index is (x,y,z):

  [ 11, 12, 13, 14, 15, 16, 101, 102, 103, 104, 105, 106 ]
  

• assuming the index is (x,z,y):

  [ 11, 12, 13, 101, 102, 103, 14, 15, 16, 104, 105, 106 ]
  

Applying to flattened data

Row Major 3D indexing

Indexing the returned array in 3 dimensions using row-major conventions appropriately would be (assuming y,x,z indexing)

data[i][j][k] = flattened_data[ depth*(j+cols*i) + k ]

where depth is the size of the 3rd dimension (z) and cols is the size of the second second dimension (x) and

i = row iterator (first dimension)
j = column iterator (second dimension)
k = depth iterator (thrid dimension)

For example, using the grid described above, extracting the linear data into a y,x,z 3D array would look like this

data[i][j][k] = flattened_data[ 2*(j+3*i) + k ]

which implemented in loops looks like

for( i = 0; i < 3; i++){
  for( j = 0; j < 2; j++){
    for( k = 0; k < 2; k++){
         data[i][j][k] = flattened_data[ depth*(j+cols*i) + k ]
         data[i][j][k] = flattened_data[  2*3*i + 2*j + k ]
     }   
  }
}

Column Major 3D indexing

Indexing with column major conventions would look like (assuming y,x,z indexing)

data[i][j][k] = flattened_data[ rows*cols*k + rows*j + i ]

where rows is the size of the first dimension, and the rest of the of the variables are as before.

Using our example grid, this looks the following

data[i][j][k] = flattened_data[ 2*3*k + 2*j + i ]

which implemented in loops looks like

for( i = 0; i < 3; i++){
  for( j = 0; j < 2; j++){
    for( k = 0; k < 2; k++){
         data[i][j][k] = flattened_data[ rows*cols*k + rows*j + i  ]
         data[i][j][k] = arr[  2*3*k + 2*j + i  ]
     }   
  }
}  

Additions

• bmi_grid.f90 module containing a grid data structure similar to the python structured developed in Python bmi grid #493
• Gridded input and output variables in the fortran bmi test model
• suite of fortran adapter unit tests interacting with the gridded variables

Removals

Changes

• minor refactoring and upgrading of iso_c_binding module
• minor refactoring of existing unit tests

Testing

1. All existing, refactored, and new unit tests for the Fortran Adapter tests are passing locally.

Todos

• Need to follow up with an issue on the CSDMS fortran BMI repository concerning the shadowing of intrinsic functions in BMI dummy variables.
• Include any relevant documentation in the BMI conventions doc in Bmi conventions doc #507

Checklist

• PR has an informative and human-readable title
• Changes are limited to a single goal (no scope creep)
• Code can be automatically merged (no conflicts)
• Code follows project standards (link if applicable)
• Passes all existing automated tests
• Any change in functionality is tested
• New functions are documented (with a description, list of inputs, and expected output)
• Placeholder code is flagged / future todos are captured in comments
• Project documentation has been updated (including the "Unreleased" section of the CHANGELOG)
• Reviewers requested with the Reviewers tool ➡️

Target Environment support

• Linux
• MacOS

[mattw-nws]

Per discussion, I think we should actually not allow this (and related updates in this commit), and rather fail if get_var_itemsize is reporting zero... but yes this does seem to be the source of the FPE.

[hellkite500]

This needs fixed!!! sizeof will give you the number of bytes of a scalar, but when called on an array, it gives the number of bytes of the the entire array...

[mattw-nws]

Pretty sure this will pass the remaining two tests after merge... will rollback if somehow not.

[mattw-nws]

FTR: yup. https://github.com/NOAA-OWP/ngen/actions/runs/5093358554