Particles crossing North Pole - follow up on #2598 and #2521

[noemieplanat]

Parcels version

v4

Description

Hi @fluidnumericsJoe,
Following up on your fix for #2598 and #2521, I have tried advecting particles through the North Pole. Advecting only one particle, the trajectory looks weird (see code and image below), and I have the impression the interpolation of velocities is not correctly done when crossing the Pole ... ? Any idea why this is the case ?
Let me know if I can test further this issue or simplify the code,
Noemie

Code sample

from parcels import FieldSet, ParticleSet, Variable, convert, ParticleFile, StatusCode # type: ignore
from parcels.kernels import AdvectionRK4_3D
import xarray as xr
from datetime import timedelta as delta
import numpy as np
import datetime
import warnings
import pylab as plt
from cartopy import crs as ccrs
import cftime as cf
import pandas as pd 

listeU = ['/storage/nplanat/CREG_REF12/1d/CREG12.L75-REF12_y2001m04.1d_gridU.nc']
listeV = ['/storage/nplanat/CREG_REF12/1d/CREG12.L75-REF12_y2001m04.1d_gridV.nc']
listeW = ['/storage/nplanat/CREG_REF12/1d/CREG12.L75-REF12_y2001m04.1d_gridW.nc']
mesh_file_z = '/storage/nplanat/Data_CREG_masks/CREG12.L75-REF09_domain_cfg_20230801_Z.nc'

ds_fields = xr.open_mfdataset(
    listeU +listeV+listeW,
    data_vars={"vozocrtx","vomecrty", "vovecrtz"},
    coords="minimal",
    compat="override",
)
ds_fields = ds_fields.drop({"depthu_bounds", \
                           "time_centered_bounds",\
                          "time_counter_bounds",\
                           "e3u","e3v", "sozotaux",\
                           "utau_atmoce", "utau_iceoce",\
                           "uwspd10", "ahtu_3d", "depthv_bounds",\
                           "sometauy", "vtau_atmoce",\
                           "vtau_iceoce","vwspd10",\
                           "ahtv_3d", "depthw_bounds",\
                           "votkeavt","e3w"}
                          )
ds_coords = xr.open_dataset(mesh_file_z, decode_times=False)
ds_coords = ds_coords.rename({'t':'time_counter'})
ds_coords = ds_coords.set_coords({'time_counter', 'nav_lev'})

ds_fset = convert.nemo_to_sgrid(fields=dict(U=ds_fields["vozocrtx"], V=ds_fields["vomecrty"], W = ds_fields["vovecrtz"]), coords=ds_coords)
fset = FieldSet.from_sgrid_conventions(ds_fset, "spherical")
pset = ParticleSet(fieldset=fset,   # the fields on which the particles are advected
               lon= [150],
               lat = [89.9],
               z = [50])
outfile = '/storage/nplanat/ADV_7/Debug_bi17.parquet'
output_file = ParticleFile(outfile, outputdt=delta(hours = 6))
pset.execute([AdvectionRK4_3D], runtime=delta(days=15), dt=delta(hours=6), output_file=output_file, verbose_progress=True)

dp = pd.read_parquet(outfile, engine = 'pyarrow')
f = plt.figure()
ax = f.add_subplot(1,1,1, projection = ccrs.NorthPolarStereo())
ax.plot(dp['lon'].values,dp['lat'].values,'.', ms = 4, transform = ccrs.PlateCarree())
ax.gridlines()

[fluidnumericsJoe]

Hey @noemieplanat , the plot thickens! I'll give this a test drive and take a look at the underlying velocity field in the data set as well. Do you happen to know if the velocity field components are lat-lon oriented, along curvilinear grid lines (covariant), or across curvilinear grid lines (contravariant) ?