Merge pull request #2455 from Parcels-code/unitconversion_for_cgridvelocity

Fixing unit conversion for C-grid velocity interpolation

Author: erikvansebille

src/parcels/_core/field.py

@@ -311,8 +311,10 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
             (u, v, w) = self._vector_interp_method(particle_positions, grid_positions, self)
 
         if applyConversion:
-            u = self.U.units.to_target(u, z, y, x)
-            v = self.V.units.to_target(v, z, y, x)
+            if self.U.grid._mesh == "spherical":
+                meshJac = 1852 * 60.0 * np.cos(np.deg2rad(y))
+                u = u / meshJac
+                v = v / meshJac
             if "3D" in self.vector_type:
                 w = self.W.units.to_target(w, z, y, x)
 

src/parcels/_core/fieldset.py

@@ -10,7 +10,6 @@
 import xarray as xr
 import xgcm
 
-from parcels._core.converters import GeographicPolar
 from parcels._core.field import Field, VectorField
 from parcels._core.utils import sgrid
 from parcels._core.utils.string import _assert_str_and_python_varname
@@ -309,7 +308,6 @@ def from_nemo(ds: xr.Dataset):
             ).to_attrs(),
         )
         fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical")
-        fieldset.V.units = GeographicPolar()
         if "UV" in fieldset.fields:
             fieldset.UV.vector_interp_method = CGrid_Velocity
         if "UVW" in fieldset.fields:

tests/test_advection.py

@@ -456,6 +456,31 @@ def test_nemo_curvilinear_fieldset():
     np.testing.assert_allclose(pset.lat, latp, atol=1e-1)
 
 
+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)
+
+
 @pytest.mark.parametrize("kernel", [AdvectionRK4, AdvectionRK4_3D])
 def test_nemo_3D_curvilinear_fieldset(kernel):
     data_folder = parcels.download_example_dataset("NemoNorthSeaORCA025-N006_data")