-
Notifications
You must be signed in to change notification settings - Fork 133
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
FieldSet.from_mom does not consider grid rotation #1509
Comments
A quick fix is to add rotation fields from e.g. the 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) |
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 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 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! |
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 fixedThe text was updated successfully, but these errors were encountered: