Efficient Methods for Converting ROMS Sigma Coordinates to Z Coordinates

[yj7054]

Dear all,

I am reaching out to inquire about efficient tools or methods for converting sigma-coordinate output to z-coordinates (also from z to sigma). Currently, I use MATLAB with spline interpolation at individual horizontal ocean masks; however, this approach becomes quite slow when working with high-resolution grids.

Could Fortran, Python, Ferret, or another method be more suitable for this transformation? I would greatly appreciate any suggestions or advice you might have.

Best regards,
Yi

[wilkin]

There is xroms (Python) which will do this for you https://github.com/xoceanmodel/xroms

Also, the functions stretching.m and set_depth.m in the myroms matlab repo
https://github.com/myroms/roms_matlab

And ROMS itself will write z-coordinates to the output netcdf files if you activate the Hout logical switches in roms.in ...

Code: Select all

Hout(idpthR) == F       ! z_rho              time-varying depths of RHO-points
Hout(idpthU) == F       ! z_u                time-varying depths of U-points  
Hout(idpthV) == F       ! z_v                time-varying depths of V-points  
Hout(idpthW) == F       ! z_w                time-varying depths of W-points  

[yj7054]

Thanks for your information, John.

[erikvansebille]

Hi John and others, do you also know if there is a method to calculate sigma from z (so the inverse operation of set_depth.m in roms_matlab)?

We need that for an under-the-hood implementation of ROMS (and CROCO) in our http://OceanParcels.org code. See https://github.com/OceanParcels/Parcels ... sions/1752 for the discussion of this issue.

The challenge we face is to write the vertical stretching function Cs as a function of z (instead of sigma). Do you know if that exists?

Thanks in advance for any help!

Erik

Prof dr Erik van Sebille | Professor in Oceanography & Public Engagement | Utrecht University | Department of Physics & Freudenthal Institute | www.uu.nl/staff/EvanSebille

[arango]

The stretching vertical terrain-following coordinate functions C(σ) are non-dimensional, ranging from 0 to -1. Their values are a function of the bathymetry h(lon, lat) at every grid point and the number of application vertical levels. Your question doesn't make sense to me. Please check the following information in WikiROMS: https://www.myroms.org/wiki/Vertical_S-coordinate

[wilkin]

Hello Erik,

I don't quite follow what you want the inverse relation for. Maybe you are rewriting the vertical advection for parcels in the native ROMS s-coordinate for better accuracy, as we did with the ROMSpath* offline particle tracking. But if that's the case you need ROMS Hz (layer thicknesses) which is the forward transform, not the reverse.

Computing the fractional s-coordinate from z is something we do for data assimilation, to map observation depth (meters) to the corresponding vertical k-index cell. We've only ever done that by interpolation, by knowing the discrete z(k) and finding the points that bracket the observation depth.

*Hunter, E.J., Fuchs, H.L., Wilkin, J.L., Gerbi, G.P., Chant, R.J. and Garwood, J.C., 2022. ROMSPath v1. 0: offline particle tracking for the Regional Ocean Modeling System (ROMS). Geoscientific Model Development, 15(11), pp.4297-4311, https://gmd.copernicus.org/articles/15/4297/2022/

John.

[erikvansebille]

Thanks @wilkin and @arango, for the quick responses. I had of course seen the wikiroms page and think I understand well how to transform from sigma to z. However, for our Parcels code, I need the reverse transform...

The reason for that is that virtual particles in Parcels 'live' in _physical space_ (lon, lat, depth, time), and that their movement is also controlled in this physical space. With this choice, users can combine flow fields from different models on different grids, and also add 'behaviour' on particles such as sinking, swimming, fragmentation, etc, through what we call kernels.

When the particle trajectories are computed in sigma-space, then none of the kernels built by others will work anymore. So for compatibility, it would be imperative that the particle advection is calculated in physical space. We can already do that (see https://docs.oceanparcels.org/en/latest ... co_3D.html), but only under the (erroneous) assumption that the transformation from z to sigma is linear.

Now it will be relatively straightforward to also incorporate a nonlinear mapping; but I'd need the equation for that...

Note that Parcels has been used in more than 200 articles already (https://oceanparcels.org/articles.html), on topics from water mass transport to fisheries and oil spill modelling to marine plastic litter transport; the community is growing fast. Being able to run Parcels directly on ROMS output (without interpolating to z-grid, which will lead to errors) would probably benefit everyone!

Any help is hugely appreciated!

[arango]

Have you examined the Rutgers ROMS Nonlinear/step_floats.F module? It has its 4th-order predictor/corrector time-stepping kernel for every Lagrangian particle (in the order of millions) to determine their trajectory, which is independent of the vertical coordinate system and ocean model. The background dynamics is a mathematical interpolator of ROMS state vector (ubar, vbar, u, v, W, t, s) that you can get directly from ROMS time-stepping or offline from a NetCDF file. However, you can use any other model, such as CROCO or NEMO. The interpolation is done in (lon, lat, z) space. It has random walk terms and biological behavior. We parallelize over the number of floats.

Andy Moore wrote this module's tangent linear and adjoint kernels to assimilate Lagrangian drifters into ROMS.

I don't see why you need to explore the inverse relationship for Vtransform=2 and Vstretching=4, which I believe is what CROCO uses. However, I haven't seen that code in years. The algebra is messy: S = (z - ζ) / (ζ + h) since all your source model data will need to have the same free surface (ζ) and bathymetry (h), which is very unlikely. We solved that problem for you more than 20 years ago.

I think the kernels that we have step_floats.F, tl_step_floats.F, and ad_step_floats.F is a wiser and generic approach to simulate Lagrangian particles.

[erikvansebille]

Thanks @arango; this is a useful response!

I assume you mean the code at https://github.com/myroms/roms/blob/5a1 ... #L274-L343? I'll have a closer look at this in the coming days and hopefully can repeat some of the methods you've used for the step_floats.F module

[arango]

Yes, that routine is the full-time-stepping kernel. You can also check interp_floats.F, vwalk_floats.F, and Biology/biology_floats.F. There are other packages to compute offline tracking lagrangian particles that you may examine. I am not aware of tracking done with Python. We experience that Phyton is much slower than Fortran by several orders of magnitude when dealing with I/O and large matrices.

Sasha Shcheptkin recently illustrated how slow Python is with the simple multiplication of large matrices compared with Fortran for a class with students. Of course, you need good programming skills to write efficient code that reduces the number of instructions in the assembler code. I think Phyton is inefficient for this type of problem when simulating millions of Lagrangian trajectories. Using C++ or Fortran will be more attractive because of the high level of interpolation!

[wilkin]

I don't see why your inverse transform can't be an interpolation by table lookup. You have z(k) and k, so just do a linear interpolation on k and z (both are guaranteed monotonic). This will be robust for all ROMS output files because the ROMS ocean_s_coordinate_g1 etc. have CF-Convention descriptions that I believe xarray knows about.

In the code snipet you quoted, this unnerves me ...

Note that in the code below we use the w velocity field for vertical velocity. However, it is unclear whether this is always the right choice. CROCO (and ROMS) also output an omega field, which may be more appropriate to use.

It's not a case of more appropriate, it's wrong to conflate the two. You should also plan for users having saved to output only one vertical velocity (either), but not both.