Curvilinear sigma coordinate grid in Opendrift

[PabloMLo]

Dear all. We are trying to use the 3D output data from the Delft3D modelling software for particle advection in Opendrift. Delft3D uses a curvilinear sigma coordinate grid (Arakawa C-grid) and we are not sure how can we give the data so Opendrift can read it correctly. The variables we have are fixed latitude and longitude coordinates, time-varying depth coordinate of the whole water column, percentage (0-1) of ten sigma layers, and time-varying water level, salinity and eastward, northward and upward currents. First, we tried the reader_netCDF_CF_generic but it seems to be intended for Z-layers as it does not recognise the sigma levels. We have seen in previous issues of the Forum that the implementation of this correction is planned. Is there any progress on this? Then, we decided to apply the reader reader_netCDF_CF_unstructured reader that seemed to work with sigma-layers, as expressed by the following formula in the internal code: sigma * (depth + elevation). To do this, we applied a pre-processing approach to our data, to convert the regular gridded data into unstructured data. We added the variables as lists of elements, taking as reference input data examples as Nordic-4km_SLEVELS_avg_00_subset2Feb2016.ncand https://thredds.met.no/thredds/dodsC/metusers/knutfd/thredds/netcdf_unstructured_samples/AkvaplanNiva_sample_lonlat_fixed.nc. The variables we specified were geographical longitude and latitude, x and y UTM coordinates, seafloor depth above geoid, sigma layer and level, as well as sea surface height above geoid and current components.

We have now received the following errors and are concerned that OpenDrift may not be able to handle curvilinear sigma layer grids.

Traceback (most recent call last):
  .............. line 39, in <module>
    reader_current = reader_netCDF_CF_unstructured.Reader('unstructureddata.nc',proj4 = proj)                    

  File "reader_netCDF_CF_unstructured.py", line 118, in __init__
    assert self.dataset.CoordinateSystem == "Cartesian", "Only cartesian coordinate systems supported"

  File "src\netCDF4\_netCDF4.pyx", line 3141, in netCDF4._netCDF4.Dataset.__getattr__
  File "src\netCDF4\_netCDF4.pyx", line 3082, in netCDF4._netCDF4.Dataset.getncattr
  File "src\netCDF4\_netCDF4.pyx", line 1532, in netCDF4._netCDF4._get_att
  File "src\netCDF4\_netCDF4.pyx", line 2028, in netCDF4._netCDF4._ensure_nc_success
AttributeError: NetCDF: Attribute not found

Can anyone give any advice or point out what the next step might be to resolve this issue? Thank you in advance.

[knutfrode]

The proper solution would be to extend reader_netCDF_CF_generic to support also sigma-coordinates. This has been desired for some time, but has unfortunately not been prioritized yet. Thus any pull request with a solution/extension would be welcome here.

The ROMS native reader does support sigma-coordinates (and converts the returned data to z-layers, as used within OceanDrift). Thus some of that code should be possible to reuse.

Another workaround in the meantime would be to post-process the ocean data to z-depths, but ideally one would like to avoid such intermediate steps.

Alternatively, I believe OceanParcels (https://oceanparcels.org) has some better support for Delft3D native output.

[PabloMLo]

Dear knutfrode Thank you very much for your quick response and suggestions. We are very motivated to use Opendrift in our research. We have managed to convert our mesh to an unstructured one and successfully run a drift simulation! When we start simulations with a single particle and zero radius, e.g. o.seed_elements(lon=lon_value,lat=lat_value,z=-5,number=1,radius=0,time=initime), the particle moves in 3D and everything is coherent. However, if we start simulations with more than 1 particle and with a non-zero radius, errors occur. In the case of only including more than 1 particle, the model shows the following error:

ERROR   opendrift.models.basemodel:1377: Mismatch of masks
Traceback (most recent call last):
  File "C:\opendrift\opendrift\models\basemodel.py", line 1372, in get_environment
    raise ValueError('Mismatch of masks')
ValueError: Mismatch of masks

Despite the error, the model iterates until it breaks with following error:

WARNING opendrift.models.basemodel:2955: The simulation stopped before requested end time was reached.
INFO    opendrift.models.basemodel:2957: ========================
INFO    opendrift.models.basemodel:2958: End of simulation:
INFO    opendrift.models.basemodel:2959: operands could not be broadcast together with shapes (148,) (10,)
INFO    opendrift.models.basemodel:2960: Traceback (most recent call last):
  File "C:\opendrift\opendrift\models\basemodel.py", line 2893, in run
    self.get_environment(list(self.required_variables),
  File "C:\opendrift\opendrift\models\basemodel.py", line 1366, in get_environment
    np.ma.mask_or(combined_mask,
  File "C:\Users\test\.conda\envs\opendrift\Lib\site-packages\numpy\ma\core.py", line 1757, in mask_or
    return make_mask(umath.logical_or(m1, m2), copy=copy, shrink=shrink)
                     ^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: operands could not be broadcast together with shapes (148,) (10,)
INFO    opendrift.models.basemodel:2961: 'The simulation stopped before requested end time was reached.'

If, in addition to using more than one particle, the radius is greater than zero, the model does not even start to run and shows the two previous errors.

We have not found any other way to feed the particle input and it is an error we cannot understand. Why does it work with only one particle and not with more than one? Could you give us some suggestions on how to proceed?

Thank you very much in advance.

[cgrimaldi]

Hi PabloMLo I am also trying to use the reader_netCDF_CF_unstructured reader to read in outputs from Delft3D Flexible Mesh output (2D output), any insight on how you adapted the code to read the FlowFM_map.nc files in? Thank you in advance Camille