Functions for grid related meta-data declaration and initialization

Remark # 1:
We have to make the distinction here between the rank of a variable (i.e. the number of informatic dimensions) and the number of physical dimensions covered by the variable. Hereafter, I will respectively use the words "rank" and "physical dimensions". Each physical dimension will be described by one coordinate array.

For example, an unstructured grid covering a 2D physical space will be a grid of physical dimension 2 but the localisation of the grid points will be expressed with arrays of rank 1, e.g. longitude(nbr_pts), latitude(nbr_pts).

Sophie> Do you agree?
Rene> Yes and no. I think that this information is not sufficient to describe unstructured grids. To cover unstructured grids in general the MpCCI specification may serve as an example. E.g. connectivity has to be specified for unstructured grids, which is known implicitely for our blockstructured grids.

Additional information is needed for unstructured grid, like the connectivity. This could be given by an additional routine such as PRISM_def_connectivity.
 
PRISM_def_grid  ( grid_id
comp_id
grid_valid_shape
grid_nodims
grid_type
ierror )

 
grid_id [OUT]  unique identifier of the grid
comp_id [IN] component Id as given by PRISM_init_comp
grid_valid_shape [IN] ilow,ihigh,[jlow,jhigh],[klow,khigh] of valid range, without halo regions
grid_nodims  [IN] to indicate whether the grid is 1-dim, 2-dim or 3-dim
grid_type [IN} PRISM integer named parameter
default: to be decided
z-regular: 1-dim array for z-coordinate is sufficient, x(i,j), y(i,j), z(k)
y-regular:
x-regular
to be continued

>Sophie: gridGenerationType: would you have only one value per grid? Which means e.g. that there would be one gridGenerationType for a 3D grid being regular in x and y and z, and a different gridGenerationType for a 3D grid being regular in x and y but irregular in z?

>Rene: We could allow for input like xz-regular, xyz-regular etc. or fixed integer parameters PRISM_xz_regular, ...

> Sophie: The gridGenerationType will define the number of physical dimensions covered by the grid (= number of coordinate arrays) but also the rank of each coordinate array.
Some examples:

Did I get this right?

Rene> Yes.

> Sophie: If this is the case, then I would say that this could also be used for unstructured grid:

Rene> See my remark about unstructured grids above.

> Sophie: If this is the case, the indication on "whether the grid is 1-dim, 2-dim or 3-dim" is included in gridGenerationType and noDimensions is no longer needed...?

Rene> The cases that have to be covered are (ignoring the unstructured grids which I think need extra calls):

1D
2D-regular 2D-irregular
3D-regular  3D-irregular
3D-regular-1st_dim
3D-regular-2nd_dim
3D-regular-3rd_dim

I am still in favor for two parameters noDimensions (which is either 1, 2, or 3) indication the number of physical dimensions and a gridGenerationType (which will be regular, irregular, regular-1st_dim, regular-2nd_dim or regular 3rd_dim )
Are there more cases to be covered?
 
 

To Do:
Sophie> As discussed this morning and by e-mail, we keep both the grid_nodims (no_dimensions) and the grid_type as arguments. We still have to agree on whether grid_type is a grid_type or an axis_type, and we still have to agree on the
related list of integer named parameters.
 
 
 

Sophie> Obviously, more work is required to define exactly the gridGenerationType supported but this way of doing allows us to follow an incremental approach which is good! We do not need to precise all these details before submitting this to IPSL.

Rene> If we get some information about what is needed for unstructured grids we could start working on something like PRISM_defineUnstructured ().

Sophie> I suppose that validBlockShape indicate the range of informatic index valid range.

Rene> Yes.

Sophie> If this is the case,  I would say that it is still possible, at this point, that validBlockShape describes the (multi-block or multi-segment) local partition in the global grid (i.e the vector "shape()" in my mail from Mon, 17 Jun).
Last remark: as the definition of the grid point location must be done in all cases, couldn't the information given here in PRISM_DefineBlock be simply added as arguments to the PRISM_SetLocation_... calls? For me it is very unnatural to define locationShape in  PRISM_DefineBlock and validBlockShape in PRISM_SetLocation_.... I would prefer to have them in only one call (as validCornerShape and CornerShape in PRISM_SetCellCorners_...). What do you think?

Rene> In principle this should be possible. But I would feel more comfortable if there is some specific grid-related information already included in the DefineBlock call to alert the code developer that this definition is related to the later calls. Otherwise the developers will make several calls to DefineBlock to get a bunch of grid Ids and use incorrect gridIds in the later calls. Also the validBlockShape info can be used to preallocate memory.
 



 
PRISM_set_center_3D ( grid_id
center_actual_shape
center_1stcoord_array
center_2ndcoord_array
center_3rdcoord_array
ierror )

 
grid_id [IN] value received from PRISM_def_grid
center_actual_shape [IN] array containing minimum and maximum index values for each informatic dimension (1, 2 or 3). This must cover at least the grid_valid_shape region but can also include halo regions.
The shape of the coord_arrays depends on the grid_type. See example for z-regular above.

>Sophie: the locations here are the locations of the centers?
With PRISM_SetLocation_Block* the grid points are defined on which the data a located. Any other suggestions for more intuitive and self explaining names are more than welcome.

>Sophie: This routine is called to describe the location of the points of a grid of physical dimension 3. But 1st-coord-array, 2nd-coord-array and 3rd-coord-array could still be of rank 1, 2, or 3 depending on the gridGenerationType? 1st-coord-array would describe the longitude, 2nd-coord-array the latitude and 3rd-coord-array the depth or height.
Did I get this right?

>Rene: Yes

>Sophie: If this is the case, why did you write "All arrays are declared in the same (i,j,k) space." This is contradictory to the
definition of gridGenerationType in which you wrote e.g. "z-regular: 1-dim array for z-coordinate is sufficient" ?

>Rene: What I wanted to say here that the space should be the same. The shape of the data array should be included for each dimension in the sahpe of the coordinate array.
Or in other words: x,y,z (1st, 2nd, 3rd dimension) associated to i,j,k should be the same for all calls/declarations.
If the valid range is from 3 to 97 than a 1 point halo is expected to range from 2 to 98 but not from 102 to 198 nor from 1 to 97, etc.



 
PRISM_set_center_2D ( grid_id
center_actual_shape
center_1stcoord_array
center_2ndcoord_array
ierror )

 
grid_id [IN] value received from PRISM_def_grid
center_actual_shape [IN] Array containing the minimum and maximum index values for each informatic dimension (1, or 2). This must cover at least the grid_valid_shape region but can also include halo regions.
The shape of the coord_arrays depends on the grid_type.

>Sophie:
-This routine is called to describe the location of the points of a grid of physical dimension 2. But 1st-coord-array, 2nd-coord-array could still be of rank 1, 2 depending on the gridGenerationType. Let's suppose the grid is of gridGenerationType=4 as defined above (unstructured horizontal grid); then locationShape would have only two elements indicating nbr_ptslow, nbr_ptshigh ?
-Why did you write "All arrays are declared in the same (i,j,k) space."
-Note 1: For horizontal logically-rectangular (i,j) grid, it would be simpler to define them with 1st-coord-array(i,j) and 2nd-coord-array(i,j) in ALL cases (even for regular lat-lon grid that could be expressed as 1st-coord-array(i) and 2nd-coord-array(j)). The other type of horizontal grid would be unstructured grid expressed as 1st-coord-array(nbr_pts) and 2nd-coord-array(nbr_pts). See Note 1 follow on.

Rene>  I personally find it to restrictive. People with pure regular grids would be forced to construct 2-d arrays. (Something like regular-z would be the only remainder.)
Again, based on my experience calls for unstructured grids in general will look much more different/require more information.

>Sophie: We agreed that regular lat-lon grid could be expressed as 1stcoord_array(i) and 2ndcoord_array(j).



PRISM_set_center_1D (

grid_id
center_actual_shape
center_coord_array
ierror )

 
grid_id [IN] value received from PRISM_def_grid
center_actual_shape [IN] Array containing the minimum and maximum index values for the informatic dimension. This must cover at least the grid_valid_shape region but can also include halo regions.

>Sophie: Why did you write "All arrays are declared in the same (i,j,k) space."

Rene> See remark above.



 
PRISM_set_4corners_2D  ( grid_id
corner_valid_shape
corner_actual_shape
corner_1stcoord_array
corner_2ndcoord_array
ierror )

 
corner_valid_shape [IN] Array containing the minimum and maximum index values for the informatic dimension of the valid range, without halo regions. 
This is needed as it may be different from grid_valid_shape
corner_actual_shape [IN] Array containing the minimum and maximum index values for the informatic dimension  of the actual range. This must cover at least corner_valid_shape, can also include halo regions

>Sophie: I do not understand the use of defining within the same call a valid range and an actual range.
There may be the case that the coordinate array containing the corners of the cell may cover more than just the valid range.
The general idea why we made a distinction between validShape and Shape (including halos) was merely to allow for putting the arrays directly into the subroutine without previous copies to some buffer.

>Sophie: Are their cases where cornerShape would be different from locationShape ?

Rene> Yes. Assuming that e.g. validLocationShape ranges from 1 to 100 validCornerShape may either range from 0 to 100 or from 1 to 101.

>Sophie: Is it supposed here that x-coord and y-coord may be of different rank depending on the gridGenerationType of that grid? I would suppose so.

Rene> Yes

>Sophie: As for FMS, if we suppose that there are only 4 vertices per cell, then these vertices can simply be expressed by x-coord and y-coord being of same rank then 1st-coord-array and 2nd-coord-array with an extent+1 for each dimension.
We could on the contrary decide that x-coord and y-coord have an additional informatic dimension of extent equal to number of vertices. If so it implies a lot of redundant information as a vertex is necessarily common to at least two meshes. But we may make this choice as it makes the definition of the vertex location easier, and allows for an easy definition of meshes having a number of vertices different from 4.
-How to define the vertices of a regular lat-lon grid? If the grid is made of rectangular meshes then it is possible to define the vertices by two vectors x-coord(i+1), y-coord(j+1). The grid vertex locations could then be reconstructed.

Rene> The default method will make such assumptions, i.e. if no call to SetCellCorner is implemented for the specific locations.
Nevertheless there is more than none way to obtain the cell corners, depending on the difference scheme which only the calling code "knows". E.g. MOM provided to different way to compute the cell corners. In a non uniform staggered grid the pressure can be located central to the velocity , or the velocity will be located central to the pressure.

>Sophie: Note 1 follow on: If we follow Note 1 suggestion, we would have on the contrary only two possibilities for x-coord and y-coord:
=> for horizontal logically-rectangular (i,j) grid  (even for regular lat-lon grid): x-coord(i,j,noCorners) and y-coord(i,j,noCorners)
=> for unstructured horizontal grid: x-coord(nbr_pts, noCorners) and y-coord(nbr_pts,noCorners).

Rene> Agree but this would mean to duplicate information and to use extra unnecessary memory. Furthermore we introduced this ...coord-regular parameters because at an ealier stage it was said that we do not want the developers to specify 2-d coord arrays x(1:idim,1:jdim) if the code works with 1-d coord arrays x(1:idim).

>Sophie: -Obviously more work is required here to define x-coord and y-coord but we do not need to precise all these details before submitting this to IPSL.

-I suppose that we should have an equivalent PRISM_SetCellCorners_Block3D(gridId, noCorners, validCornerShape, cornerShape, x-coord, y-coord, z-coord, ierror) ?

Rene> Yes. Since this was originally mend to underline some of our comments the list of calls was not complete in that way that it could be used directly for the interface specification.

>Sophie: We agreed that we would have one call to define 4 corners for each cell. In this case, the each dimension of corner_1st-coord-array and corner_2nd-coord-array has one more element as compared to the dimensions of center_1st-coord-array and center_2nd-coord-array of the PRISM_set_location_2D.

For 1D grid, we will then need:
 
PRISM_set_2corners_1D  ( grid_id
corner_valid_shape
corner_actual_shape
corner_1stcoord_array
ierror )
For 3D grid, we will need
 
PRISM_set_8corners_3D  ( grid_id
corner_valid_shape
corner_actual_shape
corner_1stcoord_array
corner_2ndcoord_array
corner_3rdcoord_array
ierror )

We would define other routines if we ever need to specify more than 4 corners for 2D grids, or more than 8 corners for 3D grids.



 
PRISM_set_mask grid_id
mask_id
mask_actual_shape
mask_array
ierror )

 
grid_id  [IN] value received from PRISM_DefineBlock
mask_id [OUT] handle to the defined mask
mask_actual_shape [IN] actual range of coordinates, must cover at least grid_valid_shape but can also include halo regions

>Sophie: What is centerShape? Are you talking here about locationShape? yes, corrected.

>Sophie: With this call, you associate one mask to one grid. This is not flexible enough as a user may want to associate different masks to one grid depending on the variable. So I think that there should be a maskId as [OUT] and that this maskId should be indicated in the variable declaration call. Added here and as input for _DefineField.

>Sophie: OK for maskId but then gridId can be removed as argument.

Rene> For my understanding the mask belongs to a grid as well. Furthermore: see below.

>Sophie: Why don't we have also a valid range (as for PRISM_SetCellCorners_Block2D)?

Rene>  I guess we do not need it. The mask describes which points of the data are valid. Since we will operate on the valid data only (and this information is already provided it is not needed again.

>Sophie: The rank of mask can be 1, 2, or 3 (and it will be possible to know it in the PSMILe thanks to the gridId)?

Rene> Yes.



 
PRISM_set_scalefactor grid_id
scalefactor_actual_shape
scalefactor_array
ierror )

 
grid_id  [IN] value received from PRISM_def_grid
scalefactor_actual_shape [IN] Array containing the minimum and maximum index values for each informatic dimension (1, 2 or 3). This must cover at least the grid_valid_shape region but can also include halo regions.

>Sophie: Why don't we have also a valid range (as for PRISM_SetCellCorners_Block2D)?

Rene> See remark for SetMask.

>Sophie: Will we have 3 scaleFactors for 3D grids, or would we have a separate PRISM_SetCellVolume?

Rene> Good point. I think SetCellAreas can cover 3D grids as well. In this case we should go for another name.
(SetCellVolume thinking of an area as a degenerated volume? )

>Sophie: For 2D grids, we need 2 scalefactors; for 3D grids, we need 3 scalefactors. I do not remember if we said that we should have two different calls
PRISM_set_scalefactor_2D (grid_id, scalefactor_actual_shape, 1st_scalefactor_array, 2nd_scalefactor_array, ierror)
and
PRISM_set_scalefactor_3D (grid_id, scalefactor_actual_shape, 1st_scalefactor_array, 2nd_scalefactor_array, 3rd_scalefactor_array, ierror)
or only one call as above with scalefactor_array having one extra dimension (as compared to the shape of center_array), this extra dimension being of extent 2 for 2D grids or of extent 3 for 3D grids. ???

>Rene: One routine PRISM_set_scalefactor should do it. How we have to interpret the scale_factor_shape array can be derived from grid_id.



 
PRISM_set_angle grid_id
angle_actual_shape
angle_array
ierror )

 
grid_id  [IN] value received from PRISM_DefineBlock
angle_actual_shape [IN] Array containing the minimum and maximum index values for each informatic dimension (1, 2 or 3). This must cover at least the grid_valid_shape region but can also include halo regions.

>Sophie: Why don't we have also a valid range as for PRISM_SetCellCorners_Block2D?

Rene> See remark for SetMask.



 
PRISM_def_var var_id
name
grid_id
mask_id
comp_id
var_nodims
var_shape
ierror )

 
var_id [OUT] return handle later use for calls to PRISM_Put/Get
name [IN] character identifier for the variable
grid_id  [IN] value received from PRISM_def_grid
mask_id [IN] mask ID as given by PRISM_def_mask
comp_id [IN] component Id as provided by PRISM_init_comp
var_nodims [IN] var_nodims(1): the overall number of dimensions of array that is going to be transferred with var_id, i.e. same than grid_nodims except for bundle variables for which it is one more.
var_nodims(2): number of bundles
var_actual_shape [IN] as for center_actual_shape but for the data arrays that are going to be transferred with PRISM_put or PRISM_get. 

Further remark for bundles: If an array is defined as a bundle only one handleId will be returned. This implies that the complete bundle will be transferred (send) to the remote component. The remote component has to be able to treat these bundles. Both components have to agree on the the precise sequence of the physical variables contained in this bundle.

There is no need to define the number of vector component as each component will be transferred through a separate call (the number of vector component will come from the PMIOD/SMIOC).

>Sophie: How do you allow here for vectors which components are in different grids? This is not covered here as gridId is only one integer?
The same problem arises if both vector components are located on the same grid. I think that the coupler or transformer has to get "knowledge" about which fields are components of a vector from the SMIOC(?) file.

>Sophie: I suggested that the number of bundles or vector component should be indicated as noDimensions(2). For me, fieldShape should describe the spatial partition and therefore should not indicates the number of bundles or vector components. This does not really make a difference however.

>rene: noDimensions(2) is perhaps are more consistent way to define the interface, since in this case all shape arrays cover only
spatial information while the distinction between scalar or vector/bundles is done with the help of a separate parameter.

>Sophie: noDimensions: isn't it necessarily the same than the number of dimensions of array used to define the grid (except for bundle variables for which it is one more) ?

Rene> noDimensions(1) Yes.

>Sophie: Do we agree that the number for bundles should be in noDimensions and not in fieldShape?

Rene> noDimension(2)

>Sophie: Do we agree that there is no need to define the number of vector component as each component will be transferred through a separate call (the number of vector component will come from the PMIOD/SMIOC).

Rene> Yes. I do not see any other way to cover a Arakawa C grid.

>Sophie: At this point, I would say that the door is still open to the possibility that fieldShape describes the (multi-block or multi-segment) local partition in the global grid (i.e. the vector "shape()" in my mail from Mon, 17 Jun).


e.) Gather/Scatter

Specifying this info e.g. with the PRISM_Put call works for the I/O. For coupling such will not work in general. Information behind this gather/scatter array is nothing else than a mask. These information is needed for neighborhood search (in case of parallel communication) and for the interpolation of data.

For coupling this information has to be specified at the same stage were we foresee the specification of the mask. If the gather/scatter changes in time interpolations matrices have to be recalculated and neighborhood search needs to be redone. All this is not required for pure I/O.

>Sophie: I am not sure I understand your remark. Do you consider (like I do) that this should be treated like any other masking transformation defined by the user in the SCC or SMIOC?

Rene> I am not completely sure about the functionality/purpose behind a gather/scatter. Currently our thinking is that we need a precise definition of the coupling fields (including the location of the data in the physical space) to perform the neighborhood calculation and calculation of the interpolation matrices. All this is not needed for I/O for which this gather/scatter may work without any problems.
Another view of this problem would be to describe such data on an unstructured grid. For I/O we may have a special I/O routine that can handle gathered or scattered data.