Using particles.dt directly in sigmagrids

Author: erikvansebille

src/parcels/kernels/_sigmagrids.py

@@ -39,47 +39,46 @@ def AdvectionRK4_3D_CROCO(particles, fieldset):  # pragma: no cover
     This kernel assumes the vertical velocity is the 'w' field from CROCO output and works on sigma-layers.
     It also uses linear interpolation of the W field, which gives much better results than the default C-grid interpolation.
     """
-    dt = particles.dt
     sigma = particles.z / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x]
 
     sig = convert_z_to_sigma_croco(fieldset, particles.t, particles.z, particles.y, particles.x, particles)
     (u1, v1) = fieldset.UV[particles.t, sig, particles.y, particles.x, particles]
     w1 = fieldset.W[particles.t, sig, particles.y, particles.x, particles]
     w1 *= sigma / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x]
-    x1 = particles.x + u1 * 0.5 * dt
-    y1 = particles.y + v1 * 0.5 * dt
-    sig_dep1 = sigma + w1 * 0.5 * dt
+    x1 = particles.x + u1 * 0.5 * particles.dt
+    y1 = particles.y + v1 * 0.5 * particles.dt
+    sig_dep1 = sigma + w1 * 0.5 * particles.dt
     dep1 = sig_dep1 * fieldset.h[particles.t, np.zeros_like(particles.z), y1, x1]
 
-    sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * dt, dep1, y1, x1, particles)
-    (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, sig1, y1, x1, particles]
-    w2 = fieldset.W[particles.t + 0.5 * dt, sig1, y1, x1, particles]
-    w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y1, x1]
-    x2 = particles.x + u2 * 0.5 * dt
-    y2 = particles.y + v2 * 0.5 * dt
-    sig_dep2 = sigma + w2 * 0.5 * dt
-    dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y2, x2]
+    sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * particles.dt, dep1, y1, x1, particles)
+    (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles]
+    w2 = fieldset.W[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles]
+    w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y1, x1]
+    x2 = particles.x + u2 * 0.5 * particles.dt
+    y2 = particles.y + v2 * 0.5 * particles.dt
+    sig_dep2 = sigma + w2 * 0.5 * particles.dt
+    dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y2, x2]
 
-    sig2 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * dt, dep2, y2, x2, particles)
-    (u3, v3) = fieldset.UV[particles.t + 0.5 * dt, sig2, y2, x2, particles]
-    w3 = fieldset.W[particles.t + 0.5 * dt, sig2, y2, x2, particles]
-    w3 *= sig_dep2 / fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y2, x2]
-    x3 = particles.x + u3 * dt
-    y3 = particles.y + v3 * dt
-    sig_dep3 = sigma + w3 * dt
-    dep3 = sig_dep3 * fieldset.h[particles.t + dt, np.zeros_like(particles.z), y3, x3]
+    sig2 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * particles.dt, dep2, y2, x2, particles)
+    (u3, v3) = fieldset.UV[particles.t + 0.5 * particles.dt, sig2, y2, x2, particles]
+    w3 = fieldset.W[particles.t + 0.5 * particles.dt, sig2, y2, x2, particles]
+    w3 *= sig_dep2 / fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y2, x2]
+    x3 = particles.x + u3 * particles.dt
+    y3 = particles.y + v3 * particles.dt
+    sig_dep3 = sigma + w3 * particles.dt
+    dep3 = sig_dep3 * fieldset.h[particles.t + particles.dt, np.zeros_like(particles.z), y3, x3]
 
-    sig3 = convert_z_to_sigma_croco(fieldset, particles.t + dt, dep3, y3, x3, particles)
-    (u4, v4) = fieldset.UV[particles.t + dt, sig3, y3, x3, particles]
-    w4 = fieldset.W[particles.t + dt, sig3, y3, x3, particles]
-    w4 *= sig_dep3 / fieldset.h[particles.t + dt, np.zeros_like(particles.z), y3, x3]
-    x4 = particles.x + u4 * dt
-    y4 = particles.y + v4 * dt
-    sig_dep4 = sigma + w4 * dt
+    sig3 = convert_z_to_sigma_croco(fieldset, particles.t + particles.dt, dep3, y3, x3, particles)
+    (u4, v4) = fieldset.UV[particles.t + particles.dt, sig3, y3, x3, particles]
+    w4 = fieldset.W[particles.t + particles.dt, sig3, y3, x3, particles]
+    w4 *= sig_dep3 / fieldset.h[particles.t + particles.dt, np.zeros_like(particles.z), y3, x3]
+    x4 = particles.x + u4 * particles.dt
+    y4 = particles.y + v4 * particles.dt
+    sig_dep4 = sigma + w4 * particles.dt
 
-    dep4 = sig_dep4 * fieldset.h[particles.t + dt, np.zeros_like(particles.z), y4, x4]
-    particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
-    particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
+    dep4 = sig_dep4 * fieldset.h[particles.t + particles.dt, np.zeros_like(particles.z), y4, x4]
+    particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * particles.dt
+    particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * particles.dt
     particles.dz += (
         (dep1 - particles.z) * 2 + 2 * (dep2 - particles.z) * 2 + 2 * (dep3 - particles.z) + dep4 - particles.z
     ) / 6