FieldSet.from_mom does not consider grid rotation

[erikvansebille]

As discovered by @eavellashaw, FieldSet.from_mom() does not consider grid rotation and thus gives erroneous results in regions where the grid is not rectilinear in the lon/lat direction. This needs to be fixed

[erikvansebille]

A quick fix is to add rotation fields from e.g. the ice.static.nc file, as use a custom RotAdvection kernel (here in Euler Forward for simplicity):

fset = fields.from_mom(...)
fset.add_field(Field.from_netcdf(rotfile, 'COSROT', {'lon': 'GEOLON', 'lat': 'GEOLAT'}))
fset.add_field(Field.from_netcdf(rotfile, 'SINROT', {'lon': 'GEOLON', 'lat': 'GEOLAT'}))

ROTparticle = JITParticle.add_variables({Variable('cosrot'), Variable('sinrot')})

def RotAdvection(particle, fieldset, time):
    particle.cosrot = fieldset.COSROT.eval(time, particle.depth, particle.lat, particle.lon, particle)
    particle.sinrot = fieldset.SINROT.eval(time, particle.depth, particle.lat, particle.lon, particle)
    u, v = fieldset.UV.eval(time, particle.depth, particle.lat, particle.lon, particle, applyConversion=False)
    ur = (u * particle.cosrot + v * particle.sinrot)
    vr = (- u * particle.sinrot + v * particle.cosrot)

    # Euler Forward step and convert from m/s to degrees/s
    particle_dlon += ur * particle.dt / (1852. * 60.0 * math.cos(particle.lat * math.pi / 180.))
    particle_dlat += vr * particle.dt / (1852. * 60.0)

[eavellashaw]

After a lot of discussions (and a bit of debugging 😉), @eavellashaw, @noemieplanat, @michaeldenes, and @erikvansebille found a way to implement and account for the grid rotation in the AdvectionRK4_3D Kernel:

fset = FieldSet.from_mom5(...)

fset.add_field(Field.from_netcdf(rotfile, 'COSROT', {'lon': 'GEOLON', 'lat': 'GEOLAT'}))
fset.add_field(Field.from_netcdf(rotfile, 'SINROT', {'lon': 'GEOLON', 'lat': 'GEOLAT'}))

ROTparticle = JITParticle.add_variables({Variable('cosrot'), Variable('sinrot')})

def Rot_AdvectionRK4_3D(particle, fieldset, time): # Modified AdvectionRK4_3D Kernel to account for grid rotation
    particle.cosrot = fieldset.COSROT.eval(time, particle.depth, particle.lat, particle.lon, particle)
    particle.sinrot = fieldset.SINROT.eval(time, particle.depth, particle.lat, particle.lon, particle)
    
    (u1, v1, w1) = fieldset.UVW.eval(time, particle.depth, particle.lat, particle.lon, particle, applyConversion=False)
    
    ur1 = u1 * particle.cosrot + v1 * particle.sinrot
    vr1 = - u1 * particle.sinrot + v1 * particle.cosrot

    lon1 = particle.lon + ur1 * .5 * particle.dt / (1852. * 60. * math.cos(particle.lat * math.pi / 180.))
    lat1 = particle.lat + vr1 * .5 * particle.dt / (1852. * 60.)
    dep1 = particle.depth + w1 * .5 * particle.dt
    
    
    particle.cosrot = fieldset.COSROT.eval(time, dep1, lat1, lon1, particle)
    particle.sinrot = fieldset.SINROT.eval(time, dep1, lat1, lon1, particle)
    
    (u2, v2, w2) = fieldset.UVW.eval(time + .5 * particle.dt, dep1, lat1, lon1, particle, applyConversion=False)

    ur2 = u2 * particle.cosrot + v2 * particle.sinrot
    vr2 = - u2 * particle.sinrot + v2 * particle.cosrot
    
    lon2 = particle.lon + ur2 * .5 * particle.dt / (1852. * 60. * math.cos(particle.lat * math.pi / 180.))
    lat2 = particle.lat + vr2 * .5 * particle.dt / (1852. * 60.)
    dep2 = particle.depth + w2 * .5 * particle.dt
    
    
    particle.cosrot = fieldset.COSROT.eval(time, dep2, lat2, lon2, particle)
    particle.sinrot = fieldset.SINROT.eval(time, dep2, lat2, lon2, particle)
    
    (u3, v3, w3) = fieldset.UVW.eval(time + .5 * particle.dt, dep2, lat2, lon2, particle, applyConversion=False)
    
    ur3 = u3 * particle.cosrot + v3 * particle.sinrot
    vr3 = - u3 * particle.sinrot + v3 * particle.cosrot
    
    lon3 = particle.lon + ur3 * particle.dt / (1852. * 60. * math.cos(particle.lat * math.pi / 180.))
    lat3 = particle.lat + vr3 * particle.dt / (1852. * 60.)
    dep3 = particle.depth + w3 * particle.dt
    
    
    particle.cosrot = fieldset.COSROT.eval(time, dep3, lat3, lon3, particle)
    particle.sinrot = fieldset.SINROT.eval(time, dep3, lat3, lon3, particle)

    (u4, v4, w4) = fieldset.UVW.eval(time + particle.dt, dep3, lat3, lon3, particle, applyConversion=False)
    
    ur4 = u4 * particle.cosrot + v4 * particle.sinrot
    vr4 = - u4 * particle.sinrot + v4 * particle.cosrot
    
    particle_dlon += (ur1 + 2*ur2 + 2*ur3 + ur4) / 6. * particle.dt / (1852. * 60. * math.cos(particle.lat * math.pi / 180.))  # noqa
    particle_dlat += (vr1 + 2*vr2 + 2*vr3 + vr4) / 6. * particle.dt / (1852. * 60.)  # noqa
    particle_ddepth += (w1 + 2*w2 + 2*w3 + w4) / 6. * particle.dt  # noqa

When comparing to CM2.5 model outputs, this agrees very well with the pathways the particles should follow, here are two examples when using Rot_AdvectionRK4_3D:

This solves the problem! The next step would be to implement it in the Parcels code directly so the user does not have to deal with this issue. Unfortunately, I don't have time to delve into this at the moment, but this fix allows MOM5 users to track particles in the Arctic, finally!