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 [email protected]

[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

[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/(pm*pn)]*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

[claudiofgcardoso]

Hello all,

Thank you for this discussion! It really helped me understand which vertical velocities from ROMS to use.
I undertand that the vertical location of depth should be at the top of the cell (s_w) but I am still strugling in understanding where the horizontal location should be loaded to parcels: if in the original coordinates at the center of the w face (rho) or if interpolated to psi points as suggested by @nvogtvincent above: "...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."

This makes perfect sense. However, when loading the 4D depth variable to the fieldset with "fieldset.from_c_grid_dataset" to then attribute depth to U, V and W using the recently-implemented method "set_depth_from_field", this field has "interp_method='cgrid_tracer'". Does this mean that during the computation of particle depth in the simulation, depth will also be interpolated based on c-grid interpolation scheme? If that's the case, the loading of depth variable should be done in the original horizontal coordinates (without previous interpolation to psi points as suggested above, right?

I tried running parcels using depth at w points (vertically) and: a) at the original horizontal rho-points and b) horizontally interpolated to psi. No error was raised with a) or b), but the results were slighly different (colors in the trajectories indicate depth):

a)

b)

Looking forward for getting some feedback on this!
Cheers,
Cláudio