Tripolar B-grid issue in the North Pole (MOM5 velocity fields)

[eavellashaw]

Tripolar B-grid issue in the North Pole (MOM5 velocity fields)

[eavellashaw]

This issue has been solved in #1509 !

[erikvansebille]

Thanks for raising this @eavellashaw. I'm surprised that this doesn't work; so far it has always worked on C-grids. See e.g. the homepage animation tutorial, where particles move happily through the Arctic in curvilinear NEMO simulation (see also the Curvilinear Grid tutorial.

Can you show what is wrong with your particles?

[eavellashaw]

Thanks so much for your fast reply @erikvansebille! A fellow student used Parcels in a NEMO configuration and it worked, that is why I am confused too ...

Here are two figures I made showing the advection of the particles, colors corresponding to the depth the particles were at at time = 0. The first one is the result when run with xu_ocean and yu_ocean, and the second one when run with geolon_c and geolat_c.

The particles are generally supposed to circulate along the slope (black arrow), not go straight toward the North Pole (1st figure). In the second one, it is also acting weirdly as the particles seem to be "blocked", but they are not that deep and the bathymetry should not be a problem at their depths.

I was thinking that when interpolating, Parcels may get the wrong indices due to the fact that the dimensions of the velocity fields are x and y, not geolon and geolat. I am fairly new to all of this so I might be completely wrong...

If there is any kind of unclarity don't hesitate to let me know!

[erikvansebille]

Dear @eavellashaw, I must say that the stop at exactly 80N in the right figure looks conspicuous. But I have no idea why this would go wrong. As you say, there is no problem with NEMO's C-Grid, so it might have to do something with MOM5's B-grid?

I'm sorry, but I can't help you much more. I don't have access to B-grid models (we don't use them in my team) and thus also no time to dig much deeper into this. Perhaps other MOM5/Parcels users might have solutions?

[eavellashaw]

I think it definitely has something to do with MOM5's B-grid, I am now reaching out to people having worked with MOM5 B-grid velocities (typically the ones mentioned in #860), I'll let you know if a solution pops out!

Thanks for your help!

[eavellashaw]

Hi @erikvansebille, as I have a bit more insight on the situation, I have a follow-up question.

So, I ran two same simulations but creating the fieldset with from_mom5() in the first one and with from_nemo() in the second one. I chose to do so to be sure that the error lied elsewhere than in the curvilinearity and tripolarity of my grid.

I don't know if you are familiar with the surface circulation in the Arctic, but the results in the second figure are way more convincing than in the first one.

My question is the following: Do you know what is the main difference when indexing the latitude and longitude coordinates of the particles between from_nemo() and from_mom5() when performing the interpolation? I think the error might come from parcels looking at the wrong coordinates.

As said earlier, here are the dimensions of my velocity fields:

u (st_ocean, yu_ocean, xu_ocean)
v (st_ocean, yu_ocean, xu_ocean)
w (sw_ocean, yt_ocean, xt_ocean)

These are 1-D arrays and when indexing, parcels might be searching for indices in these arrays and not in geolat_c and geolon_c.

This is of course an hypothesis, tell me if it makes sense and if you know where I should look into the parcels code for testing it! I already looked into the code but I am not sure where to look specifically and what to change.

Cheers,

[erikvansebille]

Thanks for diving into this, @eavellashaw. I agree with you that the second plot (with from_nemo()) looks more realistic.

The biggest difference between from_mom5 and from_nemo is that the first expects an Arakawa B-grid, whereas from_nemo interprets the data as on a C-grid.

So if your data is indeed on a C-grid, it makes much more sense to use from_nemo (or from_mitgcm, which has slightly different conventions on where the grid corners are)

Anyways, the specific calls for the MOM5 interpolation are at every instance where MOM5 is used in the code here

parcels/parcels/include/parcels.h

Lines 437 to 551 in dad04e3

	static inline StatusCode temporal_interpolation_structured_grid(type_coord x, type_coord y, type_coord z, double time, CField *f,
	GridCode gcode, int *xi, int *yi, int *zi, int *ti,
	float *value, int interp_method, int gridindexingtype)
	{
	StatusCode status;
	CStructuredGrid *grid = f->grid->grid;
	int igrid = f->igrid;
	/* Find time index for temporal interpolation */
	if (f->time_periodic == 0 && f->allow_time_extrapolation == 0 && (time < grid->time[0] || time > grid->time[grid->tdim-1])){
	return ERROR_TIME_EXTRAPOLATION;
	}
	status = search_time_index(&time, grid->tdim, grid->time, &ti[igrid], f->time_periodic, grid->tfull_min, grid->tfull_max, grid->periods); CHECKSTATUS(status);
	double xsi, eta, zeta;
	float data2D[2][2][2];
	float data3D[2][2][2][2];
	// if we're in between time indices, and not at the end of the timeseries,
	// we'll make sure to interpolate data between the two time values
	// otherwise, we'll only use the data at the current time index
	int tii = (ti[igrid] < grid->tdim-1 && time > grid->time[ti[igrid]]) ? 2 : 1;
	float val[2] = {0.0f, 0.0f};
	double t0 = grid->time[ti[igrid]];
	// we set our second time bound and search time depending on the
	// index critereon above
	double t1 = (tii == 2) ? grid->time[ti[igrid]+1] : t0+1;
	double tsrch = (tii == 2) ? time : t0;
	status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid],
	&xsi, &eta, &zeta, gcode, ti[igrid],
	tsrch, t0, t1, interp_method, gridindexingtype);
	CHECKSTATUS(status);
	if (grid->zdim == 1) {
	// last param is a flag, which denotes that we only want the first timestep
	// (rather than both)
	status = getCell2D(f, xi[igrid], yi[igrid], ti[igrid], data2D, tii == 1); CHECKSTATUS(status);
	} else {
	if ((gridindexingtype == MOM5) && (zi[igrid] == -1)) {
	status = getCell3D(f, xi[igrid], yi[igrid], 0, ti[igrid], data3D, tii == 1); CHECKSTATUS(status);
	} else if ((gridindexingtype == POP) && (zi[igrid] == grid->zdim-2)) {
	status = getCell3D(f, xi[igrid], yi[igrid], zi[igrid]-1, ti[igrid], data3D, tii == 1); CHECKSTATUS(status);
	} else {
	status = getCell3D(f, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D, tii == 1); CHECKSTATUS(status);
	}
	}
	// define a helper macro that will select the appropriate interpolation method
	// depending on whether we need 2D or 3D
	#define INTERP(fn_2d, fn_3d) \
	do { \
	if (grid->zdim == 1) { \
	for (int i = 0; i < tii; i++) { \
	status = fn_2d(xsi, eta, data2D[i], &val[i]); \
	CHECKSTATUS(status); \
	} \
	} else { \
	for (int i = 0; i < tii; i++) { \
	status = fn_3d(xsi, eta, zeta, data3D[i], &val[i]); \
	CHECKSTATUS(status); \
	} \
	} \
	} while (0)
	if ((interp_method == LINEAR) || (interp_method == CGRID_VELOCITY) ||
	(interp_method == BGRID_VELOCITY) || (interp_method == BGRID_W_VELOCITY)) {
	// adjust the normalised coordinate for flux-based interpolation methods
	if ((interp_method == CGRID_VELOCITY) || (interp_method == BGRID_W_VELOCITY)) {
	if ((gridindexingtype == NEMO) || (gridindexingtype == MOM5) || (gridindexingtype == POP)) {
	// velocity is on the northeast of a tracer cell
	xsi = 1;
	eta = 1;
	} else if (gridindexingtype == MITGCM) {
	// velocity is on the southwest of a tracer cell
	xsi = 0;
	eta = 0;
	}
	} else if (interp_method == BGRID_VELOCITY) {
	if (gridindexingtype == MOM5) {
	zeta = 1;
	} else {
	zeta = 0;
	}
	}
	if ((gridindexingtype == MOM5) && (zi[igrid] == -1)) {
	INTERP(spatial_interpolation_bilinear, spatial_interpolation_trilinear_surface);
	} else if ((gridindexingtype == POP) && (zi[igrid] == grid->zdim-2)) {
	INTERP(spatial_interpolation_bilinear, spatial_interpolation_trilinear_bottom);
	} else {
	INTERP(spatial_interpolation_bilinear, spatial_interpolation_trilinear);
	}
	} else if (interp_method == NEAREST) {
	INTERP(spatial_interpolation_nearest2D, spatial_interpolation_nearest3D);
	} else if ((interp_method == CGRID_TRACER) || (interp_method == BGRID_TRACER)) {
	if ((gridindexingtype == POP) && (zi[igrid] == grid->zdim-2)) {
	INTERP(spatial_interpolation_tracer_bc_grid_2D, spatial_interpolation_tracer_bc_grid_bottom);
	} else {
	INTERP(spatial_interpolation_tracer_bc_grid_2D, spatial_interpolation_tracer_bc_grid_3D);
	}
	} else if (interp_method == LINEAR_INVDIST_LAND_TRACER) {
	INTERP(spatial_interpolation_bilinear_invdist_land, spatial_interpolation_trilinear_invdist_land);
	} else {
	return ERROR;
	}
	// tsrch = t0 in the case where val[1] isn't populated, so this
	// gives the right interpolation in either case
	*value = val[0] + (val[1] - val[0]) * (float)((tsrch - t0) / (t1 - t0));
	return SUCCESS;
	#undef INTERP
	}

and also

parcels/parcels/include/index_search.h

Lines 101 to 105 in dad04e3

	if (zvals[zdim-1] > zvals[0]){
	if ((z < zvals[0]) && (gridindexingtype == MOM5) && (z > 2 * zvals[0] - zvals[1])){
	*zi = -1;
	*zeta = z / zvals[0];
	return SUCCESS;

But note that the code above is the C/JIT version, which may be complicated to read. If you prefer to look at the scripy version, that is here

parcels/parcels/field.py

Lines 1076 to 1080 in dad04e3

	if self.interp_method == 'bgrid_velocity':
	if self.gridindexingtype == 'mom5':
	zeta = 1.
	else:
	zeta = 0.

and

parcels/parcels/field.py

Lines 698 to 703 in dad04e3

	if z < grid.depth[0]:
	# Since MOM5 is indexed at cell bottom, allow z at depth[0] - dz where dz = (depth[1] - depth[0])
	if self.gridindexingtype == "mom5" and z > 2*grid.depth[0] - grid.depth[1]:
	return (-1, z / grid.depth[0])
	else:
	raise FieldOutOfBoundSurfaceError(0, 0, z, field=self)

[eavellashaw]

Thank you for this fast reply, I'll dive into the code and let you know if I make an interesting discovery!

I am working on an Arakawa B-grid with mom5 velocity fields (CM2-O models suite), that is why I used from_mom5(). But I thought about using from_nemo() as the biggest difference between these two are the interpolation methods and the locations of the points in the vertical. This was to check if something was wrong with my datasets, or if the problem might lie in Parcels, handling in the wrong way the coordinates of my velocity fields.

Anyway, I'll let you know!