format netCDF output accodrding to UGRID conventions

[Chilipp]

Hey @koldunovn!

as discussed in #270, here is my input for the necessary changes to format the FESOM output according to the UGRID files. My suggestions are based on the file mentioned in https://github.com/FESOM/fesom2/issues/270#issue-1104723797.

---

Mesh variable

It's all about adding a mesh variable to fesom.mesh.diag.nc (I call that variable fesom_mesh here)

Here is the

original netCDF header of fesom.mesh.diag.nc

``` netcdf fesom.mesh.diag { dimensions: nod2 = 126858 ; edg_n = 371644 ; elem = 244659 ; nz = 48 ; nz1 = 47 ; n2 = 2 ; n3 = 3 ; n4 = 4 ; N = 9 ; variables: double nz(nz) ; nz:long_name = "depth of levels" ; double nz1(nz1) ; nz1:long_name = "depth of layers" ; double elem_area(elem) ; elem_area:long_name = "element areas" ; int nlevels_nod2D(nod2) ; nlevels_nod2D:long_name = "number of levels below nodes" ; int nlevels(elem) ; nlevels:long_name = "number of levels below elements" ; int nod_in_elem2D_num(nod2) ; nod_in_elem2D_num:long_name = "number of elements containing the node" ; int nod_part(nod2) ; nod_part:long_name = "nodal partitioning at the cold start" ; int elem_part(elem) ; elem_part:long_name = "element partitioning at the cold start" ; double zbar_e_bottom(elem) ; zbar_e_bottom:long_name = "element bottom depth" ; double zbar_n_bottom(nod2) ; zbar_n_bottom:long_name = "nodal bottom depth" ; double lon(nod2) ; lon:long_name = "longitude" ; double lat(nod2) ; lat:long_name = "latitude" ; double nod_area(nz, nod2) ; nod_area:long_name = "nodal areas" ; int elements(n3, elem) ; elements:long_name = "elements" ; double nodes(n2, nod2) ; nodes:long_name = "nodal geo. coordinates" ; int nod_in_elem2D(N, nod2) ; nod_in_elem2D:long_name = "elements containing the node" ; int edges(n2, edg_n) ; edges:long_name = "edges" ; int edge_tri(n2, edg_n) ; edge_tri:long_name = "edge triangles" ; double edge_cross_dxdy(n4, edg_n) ; edge_cross_dxdy:long_name = "edge cross distancess" ; double gradient_sca_x(elem, n3) ; gradient_sca_x:long_name = "x component of a gradient at nodes of an element" ; double gradient_sca_y(elem, n3) ; gradient_sca_y:long_name = "y component of a gradient at nodes of an element" ; } ```

and based on https://ugrid-conventions.github.io/ugrid-conventions/#2d-flexible-mesh-mixed-triangles-quadrilaterals-etc-topology, this is the mesh variable that I would suggest

netcdf fesom.mesh.diag {
dimensions:
    <same-as-above>
variables:
        int fesom_mesh:
        fesom_mesh:cf_role = "mesh_topology" ;
        fesom_mesh:long_name = "Topology data of 2D unstructured mesh" ;
        fesom_mesh:topology_dimension = 2 ;
        fesom_mesh:node_coordinates = "lon lat" ;
        fesom_mesh:face_node_connectivity = "elements" ;
        fesom_mesh:edge_node_connectivity = "edges" ;
        fesom_mesh:edge_face_connectivity = "edge_tri" ;
        fesom_mesh:face_dimension = "elem" ;
        fesom_mesh:edge_dimension = "edg_n" ;
        <same-as-above>

Note that I specified the face_dimension and edge_dimension here. This would not be required if the elements variable (or edge respectively) would have shape (elem, n3) (or (edg_n, n2) respectively), as the latter is the standard ordering.

---

Variables

Then, let's have a look into the variables, for instance

the header of temp.fesom.2000.nc

``` netcdf temp.fesom.2000 { dimensions: nz1 = 47 ; nod2 = 126858 ; time = UNLIMITED ; // (1 currently) variables: double nz1(nz1) ; nz1:long_name = "depth at layer midpoint" ; nz1:units = "m" ; nz1:positive = "down" ; nz1:axis = "Z" ; double time(time) ; time:long_name = "time" ; time:standard_name = "time" ; time:units = "seconds since 2000-01-01 0:0:0" ; time:axis = "T" ; time:stored_direction = "increasing" ; float temp(time, nz1, nod2) ; temp:description = "temperature" ; temp:long_name = "temperature" ; temp:units = "C" ; } ```

Here you need to add the location and the mesh attributes (see https://ugrid-conventions.github.io/ugrid-conventions/#data-variables).

You're temp variable then becomes

netcdf temp.fesom.2000 {
dimensions:
        <same-as-above>
variables:
        <same-as-above>
    float temp(time, nz1, nod2) ;
        temp:description = "temperature" ;
        temp:long_name = "temperature" ;
        temp:units = "C" ;
        temp:mesh = "fesom_mesh" ;
        temp:location = "node" ;
}

[koldunovn]

@Chilipp Thanks a lot! This seems at the first glance as something that is not hard to implement. Even change of the ordering in elem and edg_n can be done, the file is written only once and it's not that important how quickly we save it.

Really appreciate your help, as original UGRID format description is pretty exhaustive and it's not clear what would be the required minimum.

@pgierz do you want to try to implement it, or you want me to do the first attempt? :)

[pgierz]

Thanks for the brainstorming and ideas @Chilipp :-) @koldunovn, I would not mind getting my hands dirty in FESOM, but I will probably need a guide. I actually only work on the pyfesom side, mostly. The internals in FESOM of course doesn't do any dark magic (or, well, I hope not!), but just to have someone so I can ask stupid questions without getting (too many) dirty looks.

Probably not until end of next week though, I'm building out auto-scaling for CI on AWI's HPC today and probably monday/tuesday until it actually works.

[koldunovn]

@pgierz I can give you a short intro, will not take more than 15 minutes, during next week (or even today if you have time), and then in case of difficult questions we can both ask @dsidoren. To my experience the code that needed to be changed is pretty straightforward, so adding metadata and extra variable, should not be a problem :)

[koldunovn]

Iris now has experimental support for unstructured grids, based on UGRID. Should be one of our test cases. https://scitools-iris.readthedocs.io/en/stable/whatsnew/3.2.html

[Chilipp]

Very promising @koldunovn ! They also refer to Geovista that is using pyvista as a backend (what I thought about for psyplot, too). There is also an example for visualizing FESOM data (but I did not test it)

https://github.com/bjlittle/geovista/blob/main/examples/example_from_unstructured__fesom.py