Fixing unit conversion for C-grid velocity interpolation

[erikvansebille]

This PR reverts back the change in c311fba, to avoid having to set fieldset.V.units = GeographicPolar()

Note that this aligns with the implementation of unit conversion in v3, where u and v are converted by the same meshJac

Parcels/parcels/field.py

Lines 1615 to 1635 in 1ce5908

	rad = np.pi / 180.0
	deg2m = 1852 * 60.0
	if applyConversion:
	meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1
	else:
	meshJac = deg2m if grid.mesh == "spherical" else 1
	jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
	u = (
	(-(1 - eta) * U - (1 - xsi) * V) * px[0]
	+ ((1 - eta) * U - xsi * V) * px[1]
	+ (eta * U + xsi * V) * px[2]
	+ (-eta * U + (1 - xsi) * V) * px[3]
	) / jac
	v = (
	(-(1 - eta) * U - (1 - xsi) * V) * py[0]
	+ ((1 - eta) * U - xsi * V) * py[1]
	+ (eta * U + xsi * V) * py[2]
	+ (-eta * U + (1 - xsi) * V) * py[3]
	) / jac

This has been the case since the very first implementation of C-grid interpolation in 19ef089

The question, however, is whether this is indeed correct to do. I never noticed it, but it may feel a bit strange to apply the same unit converter to lon and lat.

From discussions today, I now realise that our 'gold standard test', the zonal flow in nemo_curvilinear, in not a good test for this case, because v==np.close(0.) in that case, so it doesn't matter if it's decided by 1852*60 or by 1852*60*cos(lat)

I will go back to Philippe's original paper to see why we originally implemented the same conversion for u and v in Parcels v2.0. Hopefully that clarifies this....

[erikvansebille]

OK, I dug a bit deeper and now think that the original interpolation in v3 (and v4) is correct. See the new unit test I created in 5436136.

Parcels/tests/test_advection.py

Lines 459 to 481 in 5436136

	def test_nemo_curvilinear_with_vvel():
	"""Testing that a constant meridional velocity field in a NEMO curvilinear grid
	results in correct particle movement.
	Only works in Southern hemisphere (away from the tripolar rotated grid).
	"""
	data_folder = parcels.download_example_dataset("NemoCurvilinear_data")
	ds_fields = xr.open_mfdataset(
	data_folder.glob("*.nc4"),
	data_vars="minimal",
	coords="minimal",
	compat="override",
	)
	ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True)
	vvel = 0.5 # m/s
	ds_fields["V"].data[:] = vvel # set a constant meridional velocity
	fieldset = parcels.FieldSet.from_nemo(ds_fields)
	lonp, latp = 30, -70
	pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)
	pset.execute(AdvectionEE, runtime=np.timedelta64(10, "D"), dt=np.timedelta64(1, "D"))
	v = (pset.lat - latp) / pset.time * 1852 * 60
	np.testing.assert_allclose(v, vvel, atol=1e-5)

This sets a meridional vvel=0.5 on all grid points of the fieldset.V of NemoCurvilinear_data, and then tests whether the average velocity of particles advected indeed is vvel (accounting for the fact that one degree of latitude is 1852*60 meters).

The test passes with the current v4 (and also v3; I checked manually), but not with the change below in /src/parcels/_core/field.py

@@ -312,8 +312,8 @@ class VectorField:
 
         if applyConversion:
             if self.U.grid._mesh == "spherical":
-                meshJac = 1852 * 60.0 * np.cos(np.deg2rad(y))
-                u = u / meshJac
+                meshJac = 1852 * 60.0
+                u = u / (meshJac * np.cos(np.deg2rad(y)))
                 v = v / meshJac
             if "3D" in self.vector_type:
                 w = self.W.units.to_target(w, z, y, x)

So this means that u and v need to be converted with latitude-dependence. I find it counter-intuitive (anyone who can explain me why; please be my guest!), but the good news is that this does not seem to be a bug ☺️

[MeikeBos]

I think the same unit converter meshJac for both u and v arises from the fact that we use the fluxes to calculate the velocities for the c-grid, see equation (7) in Philippe's original paper. For the U-flux, the velocity conversion gives the deg2m * math.cos(rad * y) factor and the edge gives the deg2m factor, for the V-flux, the velocity conversion gives the deg2m factor and the edge gives the deg2m * math.cos(rad * y). Thus the total conversion factor for both the U and V flux is meshJac = deg2m * deg2m * math.cos(rad * y)

[erikvansebille]

Thanks for this explanation, @MeikeBos. Makes sense, yes!

It's a relief that we now also understand why the cosine-factor is also needed for the meridional velocity 😃