S-coordinates - issues with top and bottom cells

[nvogtvincent]

Firstly, thank you so much for producing OceanParcels (and all of the documentation) - fabulous work!

I'm trying to use fields from Croco (a version of ROMS) which uses an s-coordinate system. Following this discussion and using the s-coordinate functions described on this page (currently using the old ROMS function for legacy reasons), I've produced a sigma-depth file using the code at the bottom of this post (adapted from the code Philippe wrote in the discussion above).

The problem is that whilst this is running, I'm getting "OutOfBoundsError" when I initialise shallow (<5m) particles. I've been sloppy and for testing purposes I'm currently using a snapshot as my zeta field rather than a time-mean or time-evolving zeta, but this doesn't explain why particles initialised at -4m are 'out of bounds', and it isn't related to the vertical velocity because I have the same issue with 2D advection.

I'm pretty certain that the reason why this isn't working is because I'm currently using the sigma/stretching function values at the rho points rather than the zeta points, so the shallowest values in my depth-sigma field is as deep as -3.5m rather than zeta which is what I assume it's supposed to be. But I'm not sure what to do about this because Croco returns u/v fields at the rho depth level, and if I change the coordinates I'm inputting into the dimensions list to use the n+1 w depths rather than the n u/v/rho depths, I get a "index X is out of bounds..." error.

I'm a bit perplexed at the moment, assistance would be greatly appreciated!

from parcels import FieldSet, ParticleSet, Field, JITParticle, AdvectionRK4_3D, AdvectionRK4, CurvilinearSGrid, RectilinearSGrid, plotTrajectoriesFile
from netCDF4 import Dataset
from datetime import timedelta
from glob import glob 
import numpy as np
from scipy.interpolate import interp2d, RectBivariateSpline
from datetime import timedelta as delta
from os import path

# Import files
data_path = 'XXX'
his = data_path + 'history_file.nc'
grd = data_path + 'grid_file.nc' # Not used

# Load in coordinates
fh = Dataset(his, mode='r+')
lon_rho = fh.variables['lon_rho'][:]
lat_rho = fh.variables['lat_rho'][:]
lon_u = fh.variables['lon_u'][:]
lat_u = fh.variables['lat_u'][:]
lon_v = fh.variables['lon_v'][:]
lat_v = fh.variables['lat_v'][:]
s_rho = fh.dimensions['s_rho'] # Sigma rho
s_rho_d = fh.variables['s_rho'] # Sigma rho (variable-ified)
Cs_r = fh.variables['Cs_r'][:] # Stretching function (of sigma) rho
t = fh.variables['time'][:]
h = fh.variables['h'][:]
zeta = fh.variables['zeta'][:]
h_c = fh.variables['hc'] # Vertical stretching parameter
xi_rho = fh.dimensions['xi_rho']
eta_rho = fh.dimensions['eta_rho']
zdim = s_rho.size

fh.close

# Construct psi grid
[lon_psi_, lat_psi_] = np.meshgrid(lon_u[0,:], lat_v[:,0])
depth_psi = np.zeros((Cs_r.size, lon_psi_.shape[0], lon_psi_.shape[1]))

for k in range(zdim):
    #depth_rho = zeta[0, :, :]*(1+s_w_d[k]) + h_c*s_w_d[k] + (h - h_c)*Cs_w[k] # Convert sigma to depth - check
    S_ = h_c*s_rho_d[k] + (h - h_c)*Cs_r[k]
    depth_rho = S_ + zeta[0,:,:]*(1+(S_/h))
    depth_interp = RectBivariateSpline(lat_rho[:,0], lon_rho[0,:], depth_rho)
    depth_psi[k,:,:] = depth_interp(lat_v[:,0], lon_u[0, :],)
    # depth_psi[k,:,:] = depth_interp(lon_psi1, lat_psi1)

# Question still remains of how to apply this to u/v
if 'sigma_psi' in fh.variables:
    pass
else:
    sigma_psi = fh.createVariable('sigma_psi', 'f', dimensions=('s_rho', 'eta_v', 'xi_u'))
    lat_psi = fh.createVariable('lat_psi', 'f', dimensions=('eta_v', 'xi_u'))
    lon_psi = fh.createVariable('lon_psi', 'f', dimensions=('eta_v', 'xi_u'))

    lat_psi[:] = lat_psi_ # Write depth to file
    lon_psi[:] = lon_psi_ # Write depth to file
    sigma_psi[:] = depth_psi # Write depth to file

filenames = {'U': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'V': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'W': {'lon': grd, 'lat': grd, 'depth': his, 'data': his}}

variables = {'U': 'u',
             'V': 'v',
             'W': 'w'}

dimensions = {'U': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'V': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'W': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'}}

fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions)

[erikvansebille]

Thanks for reporting this Issue. I'd be happy to have a closer look into what's going on, but it would be much faster if you can also give me access to the netcdf files you're using. Can you send them to me, e.g. via a wetransfer link or similar? My email address is e.vansebille@uu.nl

[nvogtvincent]

Thanks Erik, I'm out of the office at the moment so haven't got access to my files but I'll send them over tomorrow!

[nvogtvincent]

Further to the email I sent you, I should have mentioned that I had to make a slight change to parcels/plotting.py to get the script to run. I had to change line 301 to:

return field.grid.lon.shape[0], 0, field.grid.lat.shape[0], 0

Because field.grid.lon.shape and field.grid.lat.shape appear to be 1D when this line executes. Not sure if that's a bug or just something I've messed up!

[erikvansebille]

Thanks for sending the file. Having had a quick look, I think that indeed you should defined your velocities on sigma_w points, see the discussion of C-grid interpolation in our new paper at https://www.geosci-model-dev-discuss.net/gmd-2018-339/gmd-2018-339.pdf. Section 3.2 deals with sigma-grids too.

However, you'll then run into the Index out of bound error. I think we've had this before for POP B-grids, and then added some code here https://github.com/OceanParcels/parcels/blob/master/parcels/field.py#L1491-L1510 I've been trying to implemented that for your file too, but without success yet. Perhaps you want to give it a try?

Note that @delandmeterp clearly is the expert on all matters related to Grids, but he's away from email until mid-August. So if we haven't fixed it by then, he might be able to help you out

[nvogtvincent]

Sorry for the slow reply. I've also unfortunately not managed to get very far with that code; it indeed gives the data the correct dimensions so the index error disappears, but the processed fieldset is still not suitable for particle tracking and I'm really not sure why! What confuses me is that here @delandmeterp reports his code (which I've more or less copied) working smoothly with ROMS output, which uses the same grid and an almost-identical output format as Croco. So I'm not sure why it suddenly isn't working here. In any case, if I get the time then I'll keep trying but I'll probably have to wait for @delandmeterp. Thank you for your advice though!

[jaseeverett]

Thanks @Plagioclase for getting in touch. I don't know how Croco differs from the ROMS output we used, but here is the code that we used to get parcels running with ROMS output. There is a bunch of extra code that could probably be stripped down, but perhaps you see what is different to yours. Let me know how you go.

#!/srv/scratch/z9902002/anaconda3/envs/py3_parcels/bin/python

from parcels import FieldSet, Field, AdvectionRK4, ParticleSet, JITParticle, plotTrajectoriesFile, Variable, BrownianMotion2D, random
from parcels import ErrorCode
import numpy as np
from glob import glob
import time as timelib
from datetime import timedelta as delta
from datetime import datetime as datetime
import cartopy
import os
from operator import attrgetter

array_ref = int(os.environ['PBS_ARRAYID'])
#array_ref = 0

data_dir = ' /srv/scratch/z3097808/20year_run/20year_freerun_output_NEWnci/'
out_dir = (os.environ['TMPDIR'])  # number of particles to be released
print(out_dir)

repeatdt = delta(days=1)  # release from the same set of locations every X day
npart = 1000  # number of particles to be released

# Forward: 9
temp_lon_array = np.array([153.8072, 153.5873, 153.5460, 153.6929, 153.7817, 153.7955, 153.7790, 153.7062, 153.5131])
temp_lat_array = np.array([-26.0, -26.5, -27.0, -27.5, -28.0, -28.5, -29.0, -29.50, -30.00])
temp_year_array = np.arange(1994, 2016, 1)

lon_array = np.repeat(temp_lon_array, temp_year_array.size)
lat_array = np.repeat(temp_lat_array, temp_year_array.size)
year_array = np.tile(temp_year_array, temp_lat_array.size)

lon = np.repeat(lon_array[array_ref],npart)
lat = np.repeat(lat_array[array_ref],npart)

start_time = datetime(year_array[array_ref],8, 1)
end_time = datetime(year_array[array_ref]+1,5, 15)  #year, month, day,

runtime = end_time-start_time + delta(days=1)

ufiles = sorted(glob('/srv/scratch/z3097808/20year_run/20year_freerun_output_NEWnci/outer_avg_*'))
vfiles = ufiles
tfiles = ufiles
bfiles = 'EACouter_mesh_srho.nc'
mesh_mask = 'EACouter_mesh_srho.nc'

# Set diffusion constants.
Kh_zonal = 100
Kh_meridional = 100

## ######
filenames = {'U': ufiles,
             'V': vfiles,
             'temp': tfiles,
             'bathy': bfiles,
             'mesh_mask': mesh_mask}

variables = {'U': 'u',
             'V': 'v',
             'temp': 'temp',
             'bathy': 'h'}

out_file = str(out_dir)+'/'+str(year_array[array_ref])+'_Lat'+str(lat_array[array_ref])+'_For.nc'

dimensions = {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 's_rho', 'time': 'ocean_time'}
dimensions['bathy'] = {'lon': 'lon_rho', 'lat': 'lat_rho'}

indices = {'depth': [29]}

if os.path.exists(out_file):
    os.remove(out_file)

def DeleteParticle(particle, fieldset, time):
    particle.delete()

fieldset = FieldSet.from_nemo(filenames, variables, dimensions, indices, allow_time_extrapolation=True)#, transpose=True)
fieldset.add_constant('maxage', 40.*86400)
fieldset.temp.interp_method = 'nearest'

# Create field of Kh_zonal and Kh_meridional, using same grid as U
#[time, depth, particle.lon, particle.lat] # Think this order is correct for here
size4D = (30,30,fieldset.U.grid.ydim, fieldset.U.grid.xdim)
fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size4D), grid=fieldset.temp.grid))
fieldset.add_field(Field('Kh_meridional', Kh_meridional*np.ones(size4D), grid=fieldset.temp.grid))

random.seed(123456) # Set random seed

class SampleParticle(JITParticle):         # Define a new particle class
    age = Variable('age', dtype=np.float32, initial=0.) # initialise age
    temp = Variable('temp', dtype=np.float32, initial=fieldset.temp)  # initialise temperature
    bathy = Variable('bathy', dtype=np.float32, initial=fieldset.bathy)  # initialise bathy
    distance = Variable('distance', initial=0., dtype=np.float32)  # the distance travelled
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
                        initial=attrgetter('lon'))  # the previous longitude
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
                        initial=attrgetter('lat'))  # the previous latitude.

def SampleDistance(particle, fieldset, time):
    # Calculate the distance in latitudinal direction (using 1.11e2 kilometer per degree latitude)
    lat_dist = (particle.lat - particle.prev_lat) * 1.11e2
    # Calculate the distance in longitudinal direction, using cosine(latitude) - spherical earth
    lon_dist = (particle.lon - particle.prev_lon) * 1.11e2 * math.cos(particle.lat * math.pi / 180)
    # Calculate the total Euclidean distance travelled by the particle
    particle.distance += math.sqrt(math.pow(lon_dist, 2) + math.pow(lat_dist, 2))
    particle.prev_lon = particle.lon  # Set the stored values for next iteration.
    particle.prev_lat = particle.lat

def SampleAge(particle, fieldset, time):
    particle.age = particle.age + math.fabs(dt)
    if particle.age > fieldset.maxage:
        particle.delete()

def SampleTemp(particle, fieldset, time):
    particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon]

def SampleBathy(particle, fieldset, time):
    particle.bathy = fieldset.bathy[0, 0, particle.lat, particle.lon]

start_time = np.repeat(start_time,len(lon))

pset = ParticleSet.from_list(fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=start_time, repeatdt=repeatdt)

pfile = pset.ParticleFile(out_file, outputdt=delta(days=1))

kernels = pset.Kernel(AdvectionRK4) + SampleAge + SampleTemp + SampleBathy + BrownianMotion2D + SampleDistance

pset.execute(kernels, 
             dt=delta(minutes=5), 
             output_file=pfile, 
             verbose_progress=False,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle},
             endtime=end_time)

[nvogtvincent]

Thanks for your reply @jaseeverett. If I understand your code correctly, are you extracting and essentially flattening the uppermost s-layer and using that as a 2D grid?

[jaseeverett]

Do you mean here:

indices = {'depth': [29]}
....
fieldset = FieldSet.from_nemo(filenames, variables, dimensions, indices, allow_time_extrapolation=True)#, transpose=True)

We are actually only using the surface layer (29th layer) for our dispersion. So we aren't squashing it. We are only using 1 layer.

[nvogtvincent]

Perhaps I'm completely confused (which is possible) but if you're using the surface s-layer for dispersion, does the depth this corresponds to not depend on the bathymetry depth? If the u/v velocities are defined at the sigma-rho coordinates (which I thought they were), then the shallowest layer has sigma != 0 and will therefore have a dependence on h? That's what I meant by flattening. If you are focusing on surface dispersion as I believe you are then that may be a reasonable approximation but I unfortunately need the depth dimension as well!

[nvogtvincent]

Using a modified script, it appears to be working now:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Script to import a Croco history file to be used in OceanParcels

from parcels import FieldSet, Field, ParticleSet, JITParticle, AdvectionRK4_3D, plotTrajectoriesFile
from netCDF4 import Dataset
from datetime import timedelta
import numpy as np
from scipy.interpolate import RectBivariateSpline
import os

# Import files
pwd = os.path.dirname(os.path.realpath(__file__)) + '/'
his = pwd + 'croco_avg.nc.1'
grd = pwd + 'croco_grd.nc.1' 

# Load in coordinates
fh = Dataset(grd, mode='r')
lon_psi = fh.variables['lon_psi'][:]
lat_psi = fh.variables['lat_psi'][:]
fh.close

fh = Dataset(his, mode='r+')
lon_rho = fh.variables['lon_rho'][:]
lat_rho = fh.variables['lat_rho'][:]
lon_u = fh.variables['lon_u'][:]
lat_u = fh.variables['lat_u'][:]
lon_v = fh.variables['lon_v'][:]
lat_v = fh.variables['lat_v'][:]
Cs_w = fh.variables['Cs_w'][:] # Stretching function of sigma-w C(s_w)
s_w = fh.dimensions['s_w'] # Sigma coords (w)
s_w_v = fh.variables['s_w'] # Sigma coords (w) (variable-ified)
t = fh.variables['time'][:]
h = fh.variables['h'][:] # Depth
zeta = fh.variables['zeta'][:]
h_c = fh.variables['hc'] # Vertical stretching parameter
xi_rho = fh.dimensions['xi_rho'] # (x)
eta_rho = fh.dimensions['eta_rho'] # (y)
zdim = s_w.size # Number of vertical points

# Construct psi grid
depth_psi = np.zeros((zdim, lon_psi.shape[0], lon_psi.shape[1]))

for k in range(zdim):
    # This generates a field showing the depth at each point in the sigma-rho grid
    S_ = (h_c*s_w_v[k] + h*Cs_w[k])/(h_c+h)
    depth_w = zeta[0,:,:]+(zeta[0,:,:]+h)*S_

    # Interpolate the sigma-rho depth field to psi locations for OceanParcels
    depth_interp = RectBivariateSpline(lat_rho[:,0], lon_rho[0,:], depth_w)
    depth_psi[k,:,:] = depth_interp(lat_psi[:,0], lon_psi[0, :])

# Write these fields to the nc file if not done already
if 'sigma_psi' in fh.variables:
    pass
else:
    sigma_psi = fh.createVariable('sigma_psi', 'f', dimensions=('s_w', 'eta_v', 'xi_u'))
    sigma_psi[:] = depth_psi # Write sigma_psi to file

fh.close

# Here using the history file for lat-lon-depth since it contains the same
# info as the grid file

filenames = {'U': {'lon': grd, 'lat': grd, 'depth': his, 'data': his}, 
             'V': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'W': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'Zeta': {'lon': grd, 'lat': grd, 'data': his}}

variables = {'U': 'u',
             'V': 'v',
             'W': 'w',
             'Zeta': 'zeta'}

dimensions = {'U': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'V': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'W': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'Zeta': {'lon': 'lon_psi', 'lat': 'lat_psi', 'time': 'time'}}

fieldset = FieldSet.from_nemo(filenames, variables, dimensions)
lon = np.linspace(122,126,num=50)
lat = np.ones([50,])*22.5
depth= np.ones([50,])*-100

pset= ParticleSet(fieldset, JITParticle, lon=lon, lat=lat, depth=depth)
pset.show(field=fieldset.Zeta)
output_file = pset.ParticleFile(name='traj.nc', outputdt=timedelta(hours=1))
pset.execute(AdvectionRK4_3D, runtime=timedelta(days=20), 
             dt=timedelta(seconds=60), output_file=output_file)
pset.show(field=fieldset.Zeta, time = timedelta(days=20))

Modifications were required to lines 1491 and 1502 of field.py to apply these routines to C-grids. Other than that it seems to be working very well, the initial outputs at least look qualitatively reasonable.

[delandmeterp]

Well done @Plagioclase ! Could you also share the modified lines of field.py, or even opening a PR with your modifications?

Thanks!

[nvogtvincent]

1491: if self.indices['depth'][-1] == self.data_full_zdim-1 and data.shape[0] == self.data_full_zdim-1 and self.interp_method in ['bgrid_velocity', 'bgrid_w_velocity', 'bgrid_tracer', 'cgrid_velocity', 'cgrid_w_velocity', 'cgrid_tracer']:

1502:

if self.indices['depth'][-1] == self.data_full_zdim-1 and data.shape[1] == self.data_full_zdim-1 and self.interp_method in ['bgrid_velocity', 'bgrid_w_velocity', 'bgrid_tracer', 'cgrid_velocity', 'cgrid_w_velocity', 'cgrid_tracer']:

I will put in the disclaimer that whilst testing the particle tracking I have come across some "out of bounds errors" with particles emerging at positive depths. I suspect that there's probably a simple explanation for this and that the integration of the CROCO output into OceanParcels is working fine, but I haven't got the time at the moment to test this. Will update when I do!

[CKehl]

@Plagioclase Have you had any time at your dispersal to re-test this procedure in new parcels versions (>=2.1.5) ? Cheers, Christian

[nvogtvincent]

Hi @CKehl, I've not worked on this since I made these posts. I will be returning to particle tracking this summer, but that will probably not be for another few months. If you happen to find anything out in the mean time, I'd be very interested to hear. Otherwise I will report back when I give it a go!

[HugoDENISFR]

Hello, I am a bit late for this discussion but my current issues with reading ROMS files in parcels is very likely to be linked to this topic.

I have a ROMS file that I am trying to use, following the usefull scripts provided by Plagioclase and Delandmeter in issues #502 and #620. First I had to create a variable that contains depth of sigma levels (varying in time) which I did using (same nomenclature than code of plagioclase) :

depth_w = np.zeros((tdim, zdim, lon_rho.shape[0], lon_rho.shape[1])) for j in range(tdim): for k in range(zdim): S_ = (h_c*s_w_v[k] + h*Cs_w[k])/(h_c+h) depth_w[j,k,:,:] = zeta[j,:,:]+(zeta[j,:,:]+h)*S_

First thing I am not sure of is that in the Parcels documentation we can find : "To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes, which are located on the corners of the cells. (for indexing details: https://www.nemo-ocean.eu/doc/img360.png ) In 3D, the depth is the one corresponding to W nodes"

If it is right, why is it needed to interpolate depth at psi points (like plagioclase did) ? isn't the depth at w points sufficient for parcels to read the grid ? I tried boths ways but none of them gave satisfactory responses.

Then I had index issue due to the fact that u, v, w field had the following dimensions : current shape = (8, 42, 176, 325) but because dimension of 's_w' was (43,) this was leading to index issues. I "solved" them supposing that velocities were considered 0 at the bottom layer and deleting bottom depth of my depth_w grid to match sizes. Again I am not sure this is the correct way of doing it.

Finally I am experiencing problems, for instance when trying to show fields

File "<ipython-input-50-73c453234e2e>", line 1, in <module> fieldset.U.show()

File "C:\Users\Hugo\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\field.py", line 1154, in show projection=projection, land=land, vmin=vmin, vmax=vmax, savefile=savefile, **kwargs)

File "C:\Users\Hugo\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\plotting.py", line 147, in plotfield latN, latS, lonE, lonW = parsedomain(domain, fld)

File "C:\Users\Hugo\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\plotting.py", line 310, in parsedomain return field.grid.lon.shape[0], 0, field.grid.lon.shape[1], 0

IndexError: tuple index out of range

or when I try to execute the run. I always have out of bond issues whatever the position of particles (I checked and depth/lat lon is correct) but given the error above I suppose it has more to do with horizontal positioning of particles in the grid. I don't know if this issue is correlated to my definition of the depth grid but any help would be appreciated.

I can provide my python script, drifters file and roms data (data + grid files) if needed.

Hugo

[nvogtvincent]

Hi Hugo,

Just wanted to say that this is a funny coincidence because I have also just started trying to do this again and I've run into the same issues as you, and was going to make a post on this later today. I will edit this post later with more information (need to gather my thoughts together) but in short, I have some thoughts on the velocity problem with a possible solution but I'm not sure why it's working.

Secondly, I have the exact same problem with showing fields and I have no idea why. You can solve the "tuple index out of range" problem by changing the last part of that code to field.grid.lon.shape[0], 0, field.grid.lat.shape[0], 0 but this introduces further problems later on in the code including some bizarre errors where the code refuses to evaluate an 'if' statement which the console itself says is correct.

Noam

[HugoDENISFR]

Hi, good to know ! I managed to get rid of this issue too by doing the little change you are suggesting but as you say this gives birth to further problems and I didn't want to mess up more with the source code !

[nvogtvincent]

Edit: It appears that the issue with omega and w having different vertical dimensions is specific to Croco, not ROMS.

Firstly, regarding psi, my understanding is that the grid used by OceanParcels is psi in the horizontal and w in the vertical, so basically the unit cell used in OceanParcels has u, v, and w at the centre of each face of the cell. In the horizontal, the w-points are offset from the psi points so you have to carry out that interpolation to generate the correct grid.

The second point I wanted to make is that Croco (and presumably ROMS) outputs two variables associated with the vertical velocity - omega and w. Omega (W) is described as the terrain-following vertical velocity component whereas w (wvel) is described as the true vertical velocity component 'computed only for output purposes' (see ROMS wiki). I believe that the difference is that omega is perpendicular to the s-coordinate layers whereas w is the component in the true vertical. At the surface at least, the difference between the two is generally small except for tides, which seem to introduce quite a significant difference. I am not sure which is supposed to be used in OceanParcels, however there is another important difference between the two variables namely that w corresponds to the rho-levels (e.g. N vertical points) whereas omega is staggered at the w-points (e.g. N+1 vertical points). If you refer to the C-grid diagram referenced in the 3D C-grid tutorial, you'll see that the variable that's in the correct position for a C-grid is omega, not w. So I think that omega may be the variable we need to be using instead of w. Since the first and last layers of omega correspond to the sea surface and sea-floor (contrary to w), omega=0 in these layers. This may be why I was seeing particles leave the domain when I was trying to use w last year, since OceanParcels may have been extrapolating non-zero vertical velocities to the sea surface.

If you use omega instead of w, you get an error saying index 50 is out of bounds for axis 0 with size 50 (or however many levels you have). I do not understand why this is, because in a c-grid (as the documentation says) the w-points are staggered - for some reason Parcels seems to need all velocities to have the same number of vertical points? Anyway though, if you change line 2139 of field.py to include cgrid_velocity (by default it only applies to b-grids) then it works and runs.

However, this introduces yet another problem namely that the horizontal velocities are too small close to the ocean surface. This is because the above routine concatenates a layer of 0s to the 'bottom' layer of the horizontal velocities to give them the same number of vertical layers as the vertical velocity. However, in Croco (and presumably ROMS) the first vertical index is the sea-floor, not the sea surface. So it's actually adding 0s to the sea-surface velocities, not the sea-floor velocities. To solve this, you have to flip the s-coordinate dimension of your input files before running them in OceanParcels (very easy with cdo invertlev [file_in] [file_out]). For some reason, if you modify that concatenation step in field.py to add the 0s to the top layer rather than the bottom, it breaks everything.

From the testing that I've done today, this seems to work quite well - through inspection, trajectories seem to correspond well with the true depth structure of the flow. I am still concerned about the difference between w and omega. For the reasons explained above I think that omega is the variable we need, but I am not sure. If it isn't then this is problematic because whereas omega is on a c-grid, w is not.

Another concern I have is with SSH, because the ROMS grid is actually 4D (time-varying) since it depends on SSH. Although OceanParcels can apparently take 4D grids, this would increase the data size by at least 33% and this would be unfeasible for the very large domain I am working with. I think it would only be feasible if there is some way to continuously recalculate the grid online, but I don't think this capability currently exists. At the moment, I am using the mean SST to calculate the grid and I am hoping that the errors this introduces are acceptably small.

Finally, on plotting fields - yep, I can't get it to work either. plotTrajectoriesFile is working, but .show isn't. I am not sure how to fix this - I did my best but I also ended up at a point where I was tinkering with code I didn't understand.

For reference, this is the code I'm currently using (assuming that you've already flipped the vertical layers):

Code

``` #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Jun 15 11:44:14 2020 Test script to import the CROCO grid @author: noam """ from parcels import FieldSet, Field, grid, ParticleSet, JITParticle, AdvectionRK4_3D, plotting, ParticleFile, plotTrajectoriesFile from glob import glob from netCDF4 import Dataset from scipy.interpolate import RectBivariateSpline import numpy as np from datetime import timedelta as delta import numpy as np from datetime import datetime import time import os import sys from os import path def make_parcels_grid_file(lon, lat, depth): # Generates a netcdf file for the parcels grid # Create file nc = Dataset(pgrd_file, mode='w') # Global attributes nc.created = datetime.now().isoformat() nc.type = 'CROCO grid file for use in Ocean Parcels' # Dimensions nc.createDimension('x', lon_psi.shape[1]) nc.createDimension('y', lon_psi.shape[0]) nc.createDimension('z', zn) # nc.createDimension('time', t.shape[0]) - TIME is temporarily deactivated # Variables nc.createVariable('lon', 'f8', ('y', 'x'), zlib=True) nc.variables['lon'].long_name = 'longitude at psi points' nc.variables['lon'].units = 'degree_east' nc.variables['lon'][:] = lon nc.createVariable('lat', 'f8', ('y', 'x'), zlib=True) nc.variables['lat'].long_name = 'latitude at psi points' nc.variables['lat'].units = 'degree_north' nc.variables['lat'][:] = lat nc.createVariable('depth', 'f8', ('z', 'y', 'x'), zlib=True) nc.variables['depth'].long_name = 'depth at psi points' nc.variables['depth'].units = 'meter' nc.variables['depth'][:] = depth nc.close() # Data paths data_path = path.expanduser('~/croco_sim/Run_WIOMonClim_Coburg10/CROCO_FILES/') grd_file = data_path + 'croco_grd.nc' pgrd_file = data_path + 'parcels_grd.nc' avg_file = data_path + 'croco_avg.nc.24h.flip' aux_file = data_path + 'croco_avg.nc' # Temporary fix as hc was removed / avg # Load grid file fh = Dataset(grd_file, mode='r') fh.set_auto_mask(False) lon_psi = fh.variables['lon_psi'][:] lat_psi = fh.variables['lat_psi'][:] lon_rho = fh.variables['lon_rho'][:] lat_rho = fh.variables['lat_rho'][:] fh.close() # Load auxiliary file (temporary fix) fh = Dataset(aux_file, mode='r') hc = fh.variables['hc'][:] # Vertical stretching parameter fh.close() # Load averages file fh = Dataset(avg_file, mode='r') fh.set_auto_mask(False) Cs_w = fh.variables['Cs_w'][:] # Stretching function for s_w (C(s_w)) s_w = fh.variables['s_w'][:] # S coordinates (w) t = fh.variables['time'][:] h = fh.variables['h'][:] # Depth zeta = fh.variables['zeta'][:] # SSH. zn = s_w.size # Number of vertical points fh.close() # CAUTION: Replace zeta with 0 zeta = np.zeros_like(zeta[0,:,:]) # Construct psi grid (lon_psi, lat_psi, z_w) z_psi = np.zeros((zn, lon_psi.shape[0], lon_psi.shape[1])) # Interpolate h and zeta to the psi grid h_interp = RectBivariateSpline(lat_rho[:, 0], lon_rho[0, :], h) h_psi = h_interp(lat_psi[:, 0], lon_psi[0, :]) zeta_interp = RectBivariateSpline(lat_rho[:, 0], lon_rho[0, :], zeta) zeta_psi = zeta_interp(lat_psi[:, 0], lon_psi[0, :]) for k in range(zn): # Transformation: # http://croco-ocean.gitlabpages.inria.fr/croco_doc/model/model.grid.html # This generates the depth at each psi point # CAUTION: this is currently IGNORING the time-dependence on zeta S_ = (hc*s_w[k] + h_psi*Cs_w[k])/(hc + h_psi) z_psi[k, :, :] = -(zeta_psi + (zeta_psi + h_psi)*S_) # [Making depth positive] # z_psi = depth at psi points # Write coordinates to grid file if path.isfile(pgrd_file): print('Warning!') print('A parcels grid file already exists in the destination specified.') # overwrite_response = input('Overwrite? (y/n)') overwrite_response = 'y' if overwrite_response == 'y': print('Overwriting parcels grid file!') os.remove(pgrd_file) make_parcels_grid_file(lon_psi, lat_psi, z_psi) print('Writing complete.') elif overwrite_response == 'n': print('Using existing parcels grid file!') else: print('Response not understood. Aborting.') sys.exit() else: print('Writing parcels grid file!') make_parcels_grid_file(lon_psi, lat_psi, z_psi) print('Writing complete.') filenames = {'U': {'lon': pgrd_file, 'lat': pgrd_file, 'depth': pgrd_file, 'data': avg_file}, 'V': {'lon': pgrd_file, 'lat': pgrd_file, 'depth': pgrd_file, 'data': avg_file}, 'W': {'lon': pgrd_file, 'lat': pgrd_file, 'depth': pgrd_file, 'data': avg_file}} variables = {'U': 'u', 'V': 'v', 'W': 'omega'} dimensions = {'U': {'lon': 'lon', 'lat': 'lat', 'depth': 'depth', 'time': 'time'}, 'V': {'lon': 'lon', 'lat': 'lat', 'depth': 'depth', 'time': 'time'}, 'W': {'lon': 'lon', 'lat': 'lat', 'depth': 'depth', 'time': 'time'}} fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions, allow_time_extrapolation=True) pnum = 20 lons = np.linspace(52,52, num=pnum) lats = np.linspace(-14,-4, num=pnum) depths = np.linspace(10,10, num=pnum) # times = (1 + np.arange(0,pnum)) * delta(days=1).total_seconds() times = np.linspace(0,0, num=pnum) pset= ParticleSet(fieldset, JITParticle, lon=lons, lat=lats, depth=depths, time=times) pfile = ParticleFile("particles", pset, outputdt=delta(hours=6)) pset.execute(AdvectionRK4_3D, runtime=delta(days=15), dt=delta(seconds=1200), output_file=pfile) pfile.export() plotTrajectoriesFile("particles.nc"); ```

I hope this is somewhat helpful - as I say, I'm not sure if the approach I've taken is the correct one but it does seem to be broadly working. The questions I would have for those who understand OceanParcels better than me are:

1. Why do I have to add an extra layer of 0s to my horizontal velocities to give them the same number of vertical layers as my vertical velocities when vertical velocities are staggered on a c-grid? And what is this actually doing in the code, since this extra layer does actually affect the velocities?
2. How does OceanParcels know I'm using an s-grid?
3. Is the vertical velocity used by OceanParcels on an s-grid the true vertical component or the component perpendicular to the unit cell face?

[HugoDENISFR]

In my case, my ROMS output is pretty basic I think, and I don't have this variable omega so I am pretty sure there is no other way than using w.

Regarding w position, I agree there is an offset from psi points in the horizontal but I got confused because of what is written in parcels documentation at https://oceanparcels.org/gh-pages/html/ As I commented before : In 3D, the depth is the one corresponding to W nodes" so if we take it litterraly, it means the depth of the vertical grid we have to provide is the depth of w position and not psi positions.

In this case if there is no interpolation the top deptht correspond to zeta at this position and the bottom one to h at this position which I tought was the right mesh to obtain

In my case execute is not working either and I get this kind of error ; Particle P[0](lon=25.100000, lat=-34.000000, depth=5.000000, time=0.000000) Time: 347160608.0, timestep dt: 1800.000000 Out-of-bounds sampling by particle at (25.100000, -34.000000, 5.000000) whatever the position of the particle and deptht...

It thought it had to do with the error above but if execution is working for you I might have to search for another cause.

[JamiePringle]

I just want to jump in with some corrections for Plagioclase.

Plagioclase says:

"however there is another important difference between the two variables namely that w corresponds to the rho-levels (e.g. N vertical points) whereas omega is staggered at the w-points (e.g. N+1 vertical points). If you refer to the C-grid diagram referenced in the 3D C-grid tutorial, you'll see that the variable that's in the correct position for a C-grid is omega, not w."

This is not correct. Both w and omega are at the same horizontal positions as rho, but are staggered in the vertical. An easy way to confirm this is to look at the grid information in the ROMS history file, where omega and w are defined as:

`asd

float w(ocean_time, s_w, eta_rho, xi_rho) ;
    w:long_name = "vertical momentum component" ;
    w:units = "meter second-1" ;
    w:time = "ocean_time" ;
    w:standard_name = "upward_sea_water_velocity" ;
    w:grid = "grid" ;
    w:location = "face" ;
    w:coordinates = "x_rho y_rho s_w ocean_time" ;
    w:field = "w-velocity, scalar, series" ;
float omega(ocean_time, s_w, eta_rho, xi_rho) ;
    omega:long_name = "S-coordinate vertical momentum component" ;
    omega:units = "meter second-1" ;
    omega:time = "ocean_time" ;
    omega:grid = "grid" ;
    omega:location = "face" ;
    omega:coordinates = "x_rho y_rho s_w ocean_time" ;
    omega:field = "omega, scalar, series" ;

`

Where you can see both are on the vertical s-coordinate grid (the s_w in the definition). See also the figure in the Numerical Solution section of the ROMS wiki, including the figure showing the vertical grid.

as is discussed in the ROMS wiki, the boundary condition on the omega is zero at the top and bottom. w is not zero at the top because the free surface can change with time (and kinematically, d(zeta)/dt=w) and there can be a w at the bottom if the bottom is sloping, to conserve volume (w=u\dot\nabla H).

NOTE -- even if w is not zero at the surface, the particle will not leave the surface if, as is true in the model, the free surface moves up with w. Likewise, the particle should not go through the bottom even if w is non-zero at the bottom, if the code handles the divergence of the horizontal velocities correctly. The latter is tricky to get write in particle tracking code because of the interpolation of velocity from the grid to the position of the particle.

The conversion between omega and w is given by equations 12 and 13 in the web page above. It is standard s-coordinate stuff. Note that the units of omega in the history file appear to be wrong -- the Hz metric in the page above has units of Length, and so omega has units of 1/time, not length/time.

Also, for HugoDENISFR, in a hydrostatic model like ROMS, w and omega are purely diagnostic quantities and can be obtained from u, v and the time derivative of the free surface. This is what ROMS does, and I have a python code that can do it (though with some numerical inaccuracy if you save the history file in single precision). W is diagnosed from omega, u and v in the wvelocity.F file in ROMS.

I am interested in using oceanParcels for ROMS, so I am interested in what you find.

Jamie Pringle UNH

[nvogtvincent]

Hi Jamie,

Apologies if I was being unclear - I did not mean to suggest that omega and w are staggered horizontally, only that they are staggered vertically.

Although I'm surprised about the history file output you've given there. In that output, w has the same vertical dimension as omega which would suggest that they're not staggered? In my Croco output, w has the vertical dimension of s_rho, not s_w. Perhaps this is a difference between Croco and ROMS although I don't understand why that change would have been implemented.

float omega(time, s_w, eta_rho, xi_rho) ;
        omega:long_name = "averaged S-coordinate vertical momentum component" ;
        omega:units = "meter second-1" ;
        omega:field = "omega, scalar, series" ;
        omega:coordinates = "lat_rho lon_rho" ;
float w(time, s_rho, eta_rho, xi_rho) ;
        w:long_name = "averaged vertical momentum component" ;
        w:units = "meter second-1" ;
        w:field = "w-velocity, scalar, series" ;
        w:standard_name = "upward_sea_water_velocity" ;
        w:coordinates = "lat_rho lon_rho" ;

Also to clarify, when you say the units of omega are wrong, are you saying that the output of omega in the his/avg files is m/s (and that this is therefore not omega but omega*Hz) or that the output of omega in the his/avg files is not m/s?

[JamiePringle]

Dear Plagioclase --

In the ROMS version maintained by Rutgers, so the central ROMS version, w and omega are on the same horizontal and vertical grids. But CROCO could be doing something different -- but then perhaps it would be helpful to clarify what it is doing.

It is worth remembering that the models do not use w directly -- the use omega, the velocity parallel to the s-coordinate. At least in Rutgers ROMS, w (called wvel in the model code) is diagnosed for output by wvelocity.F, and you can see there how it is calculated.

On the units of omega, I think from the ROMS wiki that the units of omega are 1/s, and the text in the netcdf file is inaccurate. But I will check the code and get back to you.

Jamie

[JamiePringle]

Oh goodness, I just checked varinfo.dat, and it says that the units of omega are m^3/s. I am going to get to the bottom of this, and post a bug report to the ROMS mailing list.

From varinfo.dat: `.

'omega' ! Output 'S-coordinate vertical momentum component' 'meter3 second-1' ! [m3/s] 'omega, scalar, series' 'ocean_time' 'idOvel' 'w3dvar' 1.0d0 `

Examining the code, I think omega has units of 1/s, but have posted this question to the ROMS mailing list. For your purposes in CROCO I would look at Nonlinear/wvelocity.F as see how w is calculated from omega to answer it for your fork of the ROMS code.

For those doing this on there own, since W and omega are not always written to output, perhaps they should be diagnosed from U and V. I can share some example code, and some comments on its issues, if people desire.

Jamie

[nvogtvincent]

Thanks for investigating this. On the basis of this forum post I had assumed that (regardless of what the true units of omega are) that the variable called omega in the output did have the units of m/s but it will be interesting to see what the response is.

[JamiePringle]

I should have seen that post. I think Krystin Thyng is right, and omega is multiplied by Hz to give it units of m/s. I will try to bug Hernan to update varinfo.dat.

On which oceanParcels should use? I would have to dig into there numerics and perhaps they can jump in. But omega and w are NOT interchangeable. An easy way to see that is to imagine a steady flow with a horizontally uniform horizontal velocity which does not extend to the bottom but is over a sloping bottom. If w=0, then a parcel will not change its depth as it is advected horizontally, and will stay at a constant z. If omega=0, the particle will change its depth to partially follow the bathymetry as it is advected horizontally; it will follow the s (or sigma) coordinate that it was originally on.

It is my impression that oceanParcels can diagnose the vertical velocity from the horizontal velocity, but I am going to avoid digging into that for a bit to see if one of the developers jumps in.

Jamie

EDITED: yup, the value of omega that is written out is m/s, the netcdf file is correct, varinfo.dat is wrong.

[erikvansebille]

Thanks @JamiePringle, @Plagioclase and @HugoDENISFR. Really great to see this discussion unfolding!

I am not a ROMS/CROCO user myself, so not entirely sure of the details of these models. But I guess perhaps the most important point is this: Particles in Parcels live in 'physical' (lon, lat, depth) space. This means that (in e.g. the RK4 Kernel) non-zero velocities will change particles position in physical space too. The only time that the grid is used, is for the interpolation of the velocity to the particle position.

In other words, a non-zero W will cause the particle to move up and down the water column. And contrary, a zero W will cause the particle to stay at the same depth. @JamiePringle was hence right in that Parcels expects the w velocity, not the omega velocity.

And no, Parcels cannot diagnose W itself. So ideally you'd do that in pre-processing

Please free free to suggest changes/improvements to make Parcels better work with ROMS/CROCO. Perhaps even a Tutorial, in the style of the NEMO tutorial. Again, since we are not ROMS users I'm not sure if we can provide this. But if you write it, I'd be happy to host it on the oceanparcels tutorial list!

[HugoDENISFR]

Hi @erikvansebille, thanks for your answer, could you provide a feedback about parcels documentation regarding from_c_grid method ?

"To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes, which are located on the corners of the cells. (for indexing details: https://www.nemo-ocean.eu/doc/img360.png ) In 3D, the depth is the one corresponding to W nodes"

The depth mesh expected by parcels does correspond to the depth of W nodes or f/psi nodes (corner of the cell) as suggested by @Plagioclase ?

[erikvansebille]

Hi @HugoDENISFR. The reason is that the interpolation on C-grids (where velocities should be interpolated as fluxes across faces) is linear in only one of the dimensions: U is linearly interpolated in the lon-direction, V is linearly interpolated in the lat-direction and W is linearly interpolated in the depth direction. So Parcels needs to know the corner of the grid cell, see also Figure 4 of this GMD paper.

Perhaps it's easiest to refer to the code, so that you can see what's going on. For Scipy interpolation (JIT is identical, but in C), it's at https://github.com/OceanParcels/parcels/blob/master/parcels/field.py#L878-L882:

        elif self.interp_method == 'cgrid_velocity':
            # evaluating W velocity in c_grid
            f0 = self.data[ti, zi, yi+1, xi+1]
            f1 = self.data[ti, zi+1, yi+1, xi+1]
            return (1-zeta) * f0 + zeta * f1

Does this help?

[nvogtvincent]

Thanks for your answer, @erikvansebille. It looks to me like there are two challenges facing implementing Parcels for very large datasets in CROCO. The first issue is in common with ROMS, namely the fact that s-grids are time-varying. Calculating the time-varying 4D grids required for Parcels will be possible for short/coarse runs but highly problematic for large runs. These 4D grids would have similar dimensions to the U/V files which could be on the order of TBs, which would be computationally expensive and I'm not even sure if plain python can do that. Two ways I can think of getting around this is (1) if the errors introduced by neglecting the time-varying grid are small or (2) if it's possible to recalculate these grids online. I'm not sure how difficult the latter would be.

The second issue appears to be unique to CROCO (inherited from ROMS AGRIF) - apologies for the confusion I caused above, I didn't realise that the w output is different between CROCO and ROMS. For some reason, CROCO outputs the true vertical velocity at rho vertical levels (not w vertical levels), so the U/V/W output from CROCO isn't staggered in the vertical and is therefore not supplied at the sea surface or sea floor. @JamiePringle pointed out that w can be calculated from u/v when CROCO is run in hydrostatic mode but like the first problem, for very large datasets this would be very expensive. Can Parcels deal with a non-standard c(ish)-grid where the vertical velocity isn't staggered?

[JamiePringle]

Several comments on the above.

The fundamental issue is that the ROMS grid goes from the bottom to z=zeta, the free surface. Also, in a hydrostatic model, w and omega are purely diagnostic quantities.

So first, if CROCO is not too different from ROMS, there is no reason to store the vertical grids each output (though that option does exist in the ROMS code). The change in the vertical grid can be entirely diagnosed from the free surface elevation zeta. Please see ROMS wiki on the vertical coordinate equations 1 and 2 for the equations for the grid. The standard pyroms toolbox has these routines encoded, and for a ROMS grid object romsGrid, you can get the the z coordinate of the w points from romsGrid.vgrid.z_w(), and this routine takes zeta as an argument to give you the instantaneous z grid. Or you can write your own code to give z -- it is easy and the transformations are well documented.

Second, W can be calculated from u and v. I have a code that does so in python, as does Rob Hetland. The calculation is not very expensive. Mine follows wvelocity.F, his is pythonic and much faster. Both to be precise need to have the time derivative of the free surface zeta to the the velocity right near the free surface. You have three choices for that 1) calculate the time derivative zeta_t from the history file. This will be bad unless the history is saved very frequently. 2) Save zeta_t from the model run. This does not help you if using someone else's precomputed run. 3) ignore zeta_t in the calculation of w. This might be your best option.

If you ignore zeta_t in the calculation of w, then you MUST ignore the time variation of the z coordinate. If you think about this at the surface, if zeta_t=0, i.e. the free-surface is not changing height, then w=0 at the surface AND the location of the surface point is not changing in time. If you look at the code in wvelocity.F, you can see that this effect extends in the vertical linearly through the water column, becoming 0 at the bottom.

So it might be that your best bet, if you are working with pre-existing model files, is to compute w assuming that zeta_t=0, and then calculate z_w and z_r, the vertical grids on w and rho points, assuming zeta=0 and is unchanging.

Conversely, if you use the model's computation of w, then you MUST include the effect of the time changing vertical coordinate, for otherwise your particles will leave through the surface, and similar (though less evident) errors will occur elsewhere in the water column, most strongly near the surface. To see why, think of a particle at the free surface. If the free surface is going up, there will be a positive w at the surface. If you neglect the change in the z coordinate, the particle will go sailing off into the sky...

If you save omega, you can easily compute w either assuming a constant zeta=0 or a non-zero zeta_t.

If you calculate W or Omega from U and V, then there will be some error if the history files are saved with 32 bits instead of 64. I don't find it too bad in my tests.

Jamie

[nvogtvincent]

Jamie thank you for these suggestions. If I understand you correctly, a workflow could be as follows:

1. Run model, output U/V (+ mean ζ?)
2. Calculate W as u dz/dx + v dz/dy + omega (since I suspect the omega in the output is omega*Hz although we've not confirmed this yet, and I'm assuming that if we're approximating dζ/dt ~ 0 the dz/dx and dz/dy terms would act on the time-mean grid rather than the time-varying grid)
3. Particle tracking in Parcels using stationary grid using mean ζ

I'm not used to using python with extremely large datasets, do you know if the scripts that you/Rob Hetland have written to calculate W can deal with datasets on the order of TBs (or otherwise if this is something that is physically possible)? This is the reason why I've been hesitant about steps that require significant preprocessing of data (and why online processing is so helpful) since file size is already the limiting step in terms of what I can do, and I assume that this is a general problem.

[JamiePringle]

Roughly, but

1) write out u, v and zeta (no need to use mean zeta, just instantaneous). You could write out w or omeg as well 2) If you want, you can also calculate omega from u and v. That is what I do. In the model, omega is diagnosed from U and V, it is not a prognostic variable. My workflow is to write out U and V, and then calculate omega and w. 3) particle track with zeta=0, OR with time varying zeta and z grids.

If you have computers that can read in a single time-step of the model fields, calculating W and omega is quick, and usually much less time consuming then reading in the fields. W and omega depend on a single time slice of U and V. Certainly, any computer that can run the model can compute W and Omega. Remember it is not the size of the total data set, just the size of single time record.

Since it seems that CROCO has diverged substantially from ROMS, I STRONGLY suggest you confirm that there calculation of omega is the same as in ROMS. This should be easy to see in the source code for CROCO. In roms, omega is diagnosed in omega.F, is calculated to have units of m^3/s, and is called W, where W=[Hz/(pmpn)]omega. pm=1/dx, pn=1/dy.

I think with ocean parcels you can store the W field in a different file from the other fields.

Jamie

Caveat: I have not yet done this with oceanParcels, I have been using tracmass and ariane. I hope to use oceanParcels with ROMS sometime later this summer. I have been using tracmass with Mercator.

[HugoDENISFR]

Hi everyone, hope you had a great week-end !

First of all, Thanks a lot @Plagioclase for sharing your code with us, it helped me a lot running parcels with my ROMS dataset which gave pretty good results ! I am comparing the output with another Lagrangian model (Ichthyop) and both gave alike results (even though differencies tend to accumulate over time) which is reassuring regarding the implementation.

However I still have a small issue that and I can't figure out what I am doing wrong. I am pretty sure it has to do with the fieldset.W.set_scaling_factor(-1.) parameter. In my case I have to set it to -1 because it seems the convention used in my dataset differs from the one used in parcels (w increasing means decreasing z ) .

I think this parameter is not taken into account most of the time when I run my simulation because it doesn't change the output of the run without it and I obtain a dispersion pattern that is opposite to what I expect (for instance particles going west whereas I expect them to go east with a very similar trajectory). It may have to do with the order in which I pass the different commands (might be erased by another one) but I can't figure out what is the problem. I was once able to get the expected dispersal (particles going in the right direction according to the prediction of Ichthyop) without really knowing what I had done to make it work and since then I wasn't able to reproduce this output despite the fact that my code hasn't changed.

Here is the order in which I pass the different commands :

fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions,allow_time_extrapolation=True) fieldset.W.set_scaling_factor(-1.) pset = ParticleSet.from_list(fieldset=fieldset, pclass=JITParticle,lon=init_lon,lat=init_lat,depth=init_depth) kernels = pset.Kernel(AdvectionRK4_3D) output= pset.ParticleFile(name=out + 'test_parcels_dylan_w_adv.nc', outputdt=delta(seconds=5400)) pset.execute(kernels + k_sample, runtime=delta(hours=19.5), dt=delta(seconds=1800),output_file=output) output.export() output.close()

If any of you have an idea of what I have to change in this code, I would be glad to hear about it.

Hugo

[erikvansebille]

Dear @HugoDENISFR. Your code seems to be ok. Do you mean that the fieldset.W.set_scaling_factor(-1.) does not have any effect? Or that it's effect is the wrong way around?

In other words: do the results 'swap' if you change to fieldset.W.set_scaling_factor(1.)? Or do they stay the same?

[HugoDENISFR]

No you are right, I was wrong the outputs are slightly different. Still, I don't understand the radical change in outputs when changing AdvectionRK4 for AdvectionRK4_3D. I would expect that taking into account vertical advection would change a bit the trajectories and depth of particles (like I obtain in other models such as Ichthyop or Opendrift) and not reverse the whole trajectory (see below).

So there is clearly something I am missing and maybe not due to vertical motion as I thought before.

With AdvectionRK4 ->

With Advection RK4_3D (same inititial positions, fields, etc...) ->

[erikvansebille]

Sorry, I am not sure what I'm seeing in these figures..

Could you provide a minimal working example? That is a much more tractable way for us to explore/fix a potential bug

[nvogtvincent]

@HugoDENISFR, I had a similar issue with the horizontal current components being the wrong way round. The solution (if I recall correctly) was to change the sign of the depths you write to the depths file - see if that works.

[HugoDENISFR]

Thanks again @Plagioclase , it worked ! Seems like you are ready to make this ROMS tutorial ! ;)

Cheers,

Hugo

[HugoDENISFR]

Hello everyone,

Hope you are all well since the last discussion. I used for several simulations the code version for ROMS dataset and it worked well. However I came across several errors where particles reached the bottom of the water column which make me think again about the implementation and there is something I am a bit uncertain about :

In my ROMS dataset there are n+1 s_w levels and n s_rho levels. Thus u,v and w have n values in the vertical axis. Using your code version @Plagioclase gives me n+1 sigma level depths which is perfectly normal.

However Parcels doesn't work if the vertical dimensions of my sigma_levels and u,v,w variables don't match. I am not sure why but from what I have seen, it is due to the way NEMO grid is designed (same number of t-points and w-points in NEMO).

The work around I had use was to delete the bottom sigma levels. Therefore the sigma grid and fields had the same vertical dimension. Because of these errors I am a bit uncertain about if this is the good way of doing this or not. And because of what @JamiePringle said, adding a layer of u,v,w with null velocities to my dataset at the top or bottom doesn't seem a good idea in a physical point of view.

If you have any thoughts about this, I would be glad to hear them !

[LaurentDebreu]

Dear all, I am one of the developers of CROCO. Thanks for this very interesting discussion. I should give a try soon on the use of oceanparcels with CROCO outputs. Just a comment with respect to the location of w in the output files. As you probably now, there is not just one ROMS model but a family of ROMS models (even if you may consider the Rutgers version as the official version). CROCO is the following of ROMS_AGRIF and ROMS_AGRIF started from the UCLA branch of ROMS (and not the Rutgers branch). The UCLA ROMS outputs w at s_rho and not at s_w. So we didn't change anything. I don't know which is the best with respect to lagrangian trajectories. (both choices have to performed a vertical interpolation, either of Omega or of the horizontal divergence). And these vertical interpolations are hardcoded in the models (spline interpolation). So that if you think, that it is not the right way, as mentioned in a previous post, the best is probably to recompute
w with your preferred target vertical positioning / vertical interpolation scheme. Laurent