Temperature/Salinity cannot be read from FVCOM4.4.6 output via reader_netCDF_CF_unstructured in oil spill modelling.

[Sailor-Yao]

What I use: FVCOM 4.4.6; IDE Pycharm

Root Causes:

• Temperature/salinity data failed to pass correctly due to incorrect nodes variable naming

• Even when passed successfully, sigma-coordinate (siglev siglay) issues prevented proper data reading

Implemented Fixes:

I modified the nodes variable and added code to read ’node_variables‘ with (time, siglay, node) dimensional structure.

1. Required Variables The oil spill model's requested_variables requires these four variables: ['x_sea_water_velocity', 'y_sea_water_velocity', 'sea_water_temperature', 'sea_water_salinity']

2.Unstructured Grid Data Issue For unstructured grids, temperature (sea_water_temperature) and salinity (sea_water_salinity) are stored as node variables. The original code failed to read these because: These variables use siglay (sigma layers) rather than siglev (sigma levels) for vertical coordinates

3.Implemented Solution Added siglay support to enable proper reading of these variables. Below is the netCDF structure from my data same as AkvaplanNiva_sample_lonlat_fixed.nc:

float32 temp(time, siglay, node)
    long_name: temperature
    standard_name: sea_water_temperature
    units: degrees_C
    grid: fvcom_grid
    coordinates: time siglay lat lon
    type: data
    mesh: fvcom_mesh
    location: node
unlimited dimensions: time
current shape = (24, 34, 152842)
 
float32 salinity(time, siglay, node)
    long_name: salinity
    standard_name: sea_water_salinity
    units: 1e-3
    grid: fvcom_grid
    coordinates: time siglay lat lon
    type: data
    mesh: fvcom_mesh
location: node

Here are my modified

a. changed the standard_name in nodes variable, we can noticed the standard_name in temp and Salinity are standard_name: sea_water_temperature, standard_name: sea_water_salinity.

Without this fix, the code will raise exceptions"missing variables requested" because the len(node_variables) is Null.

    node_variables = [
        #---changed by cyao -----
        # 'salinity',
        # 'temperature',
        'sea_water_salinity',
        'sea_water_temperature',
        #--------END-----------
        'sea_floor_depth_below_sea_level',
        'sea_floor_depth_below_geoid',
        'sea_surface_height_above_geoid',
    ]

b. added the parts to read nodes Variables with (time, siglay, node) dimension. I added the codes below in this part.

# ----- added by cyao 0418 2025-----------
                elif 'siglay' in dvar.dimensions:
                    sigma_ind = self.__nearest_face_sigma__(dvar, nodes, z)
                    assert len(sigma_ind) == len(z)

                    # Reading the smallest profile slice covering the actual data
                    block = dvar[indx_nearest,
                                 slice(sigma_ind.min(),
                                       sigma_ind.max() + 1),
                                 slice(nodes.min(),
                                       nodes.max() + 1)]
                    # Picking the nearest value
                    variables[var] = block[sigma_ind - sigma_ind.min(),
                                           nodes - nodes.min()]
 #  -----------END ------------------------------------

C. Modified the part about h_center The data I obtained is from FVCOM_4.4.6, and the output variables don't include 'h_center'. This issue has been reported by other in https://github.com/OpenDrift/opendrift/issues/1385.

I took a straightforward approach—I modified reader_netCDF_CF_unstructured by directly calculating the water depth at the triangular mesh elements as the average of the depths at the three nodes.

Regarding this issue, This isn't an issue with Opendrift itself—users should adapt it to their specific needs.

Here’s the code I added:

if self.ocean_depth_nele is None:
    logger.info('Reading ocean depth center into memory...')
    #--Changed by Cyao------
    if 'h_center' in self.dataset.variables:
        self.ocean_depth_nele = self.dataset['h_center'][:]
    else:
        # Compute h_center from nodal depths and element connectivity
        h_node = self.dataset["h"][:]
        nv = self.dataset["nv"][:]
        if nv.shape[0] == 3:  # Shape (3, nele) - standard FVCOM format
            # Convert 1-based to 0-based indexing and average
            h_center = np.mean(h_node[nv - 1], axis=0)
        self.ocean_depth_nele = h_center
    #----END-------------------```

[knutfrode]

Hi, I am not using this reader myself, but if you would like to provide a pull request, I can merge it, as long as all tests and examples are passing.