Skip to content

Commit ff44d72

Browse files
erikvansebilleerikvansebille
committed
Doing time interpolation inside c-grid interpolation function
1 parent 7ce1663 commit ff44d72

File tree

1 file changed

+61
-58
lines changed

1 file changed

+61
-58
lines changed

‎parcels/field.py‎

Lines changed: 61 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -955,9 +955,9 @@ def _check_grid_dimensions(grid1, grid2):
955955
and np.allclose(grid1.time_full, grid2.time_full)
956956
)
957957

958-
def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, applyConversion=True):
958+
def spatial_c_grid_interpolation2D(self, time, z, y, x, particle=None, applyConversion=True):
959959
grid = self.U.grid
960-
(_, _, eta, xsi, _, zi, yi, xi) = self.U._search_indices(time, z, y, x, particle=particle)
960+
(tau, _, eta, xsi, ti, zi, yi, xi) = self.U._search_indices(time, z, y, x, particle=particle)
961961

962962
if grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
963963
px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]])
@@ -977,54 +977,63 @@ def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, apply
977977
c2 = i_u._geodetic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
978978
c3 = i_u._geodetic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
979979
c4 = i_u._geodetic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
980-
if grid.zdim == 1:
981-
if self.gridindexingtype == "nemo":
982-
U0 = self.U.data[ti, yi + 1, xi] * c4
983-
U1 = self.U.data[ti, yi + 1, xi + 1] * c2
984-
V0 = self.V.data[ti, yi, xi + 1] * c1
985-
V1 = self.V.data[ti, yi + 1, xi + 1] * c3
986-
elif self.gridindexingtype in ["mitgcm", "croco"]:
987-
U0 = self.U.data[ti, yi, xi] * c4
988-
U1 = self.U.data[ti, yi, xi + 1] * c2
989-
V0 = self.V.data[ti, yi, xi] * c1
990-
V1 = self.V.data[ti, yi + 1, xi] * c3
991-
else:
992-
if self.gridindexingtype == "nemo":
993-
U0 = self.U.data[ti, zi, yi + 1, xi] * c4
994-
U1 = self.U.data[ti, zi, yi + 1, xi + 1] * c2
995-
V0 = self.V.data[ti, zi, yi, xi + 1] * c1
996-
V1 = self.V.data[ti, zi, yi + 1, xi + 1] * c3
997-
elif self.gridindexingtype in ["mitgcm", "croco"]:
998-
U0 = self.U.data[ti, zi, yi, xi] * c4
999-
U1 = self.U.data[ti, zi, yi, xi + 1] * c2
1000-
V0 = self.V.data[ti, zi, yi, xi] * c1
1001-
V1 = self.V.data[ti, zi, yi + 1, xi] * c3
1002-
U = (1 - xsi) * U0 + xsi * U1
1003-
V = (1 - eta) * V0 + eta * V1
1004-
rad = np.pi / 180.0
1005-
deg2m = 1852 * 60.0
1006-
if applyConversion:
1007-
meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1
1008-
else:
1009-
meshJac = deg2m if grid.mesh == "spherical" else 1
1010-
1011-
jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
1012-
1013-
u = (
1014-
(-(1 - eta) * U - (1 - xsi) * V) * px[0]
1015-
+ ((1 - eta) * U - xsi * V) * px[1]
1016-
+ (eta * U + xsi * V) * px[2]
1017-
+ (-eta * U + (1 - xsi) * V) * px[3]
1018-
) / jac
1019-
v = (
1020-
(-(1 - eta) * U - (1 - xsi) * V) * py[0]
1021-
+ ((1 - eta) * U - xsi * V) * py[1]
1022-
+ (eta * U + xsi * V) * py[2]
1023-
+ (-eta * U + (1 - xsi) * V) * py[3]
1024-
) / jac
1025-
if isinstance(u, da.core.Array):
1026-
u = u.compute()
1027-
v = v.compute()
980+
981+
def _calc_UV(ti, yi, xi):
982+
if grid.zdim == 1:
983+
if self.gridindexingtype == "nemo":
984+
U0 = self.U.data[ti, yi + 1, xi] * c4
985+
U1 = self.U.data[ti, yi + 1, xi + 1] * c2
986+
V0 = self.V.data[ti, yi, xi + 1] * c1
987+
V1 = self.V.data[ti, yi + 1, xi + 1] * c3
988+
elif self.gridindexingtype in ["mitgcm", "croco"]:
989+
U0 = self.U.data[ti, yi, xi] * c4
990+
U1 = self.U.data[ti, yi, xi + 1] * c2
991+
V0 = self.V.data[ti, yi, xi] * c1
992+
V1 = self.V.data[ti, yi + 1, xi] * c3
993+
else:
994+
if self.gridindexingtype == "nemo":
995+
U0 = self.U.data[ti, zi, yi + 1, xi] * c4
996+
U1 = self.U.data[ti, zi, yi + 1, xi + 1] * c2
997+
V0 = self.V.data[ti, zi, yi, xi + 1] * c1
998+
V1 = self.V.data[ti, zi, yi + 1, xi + 1] * c3
999+
elif self.gridindexingtype in ["mitgcm", "croco"]:
1000+
U0 = self.U.data[ti, zi, yi, xi] * c4
1001+
U1 = self.U.data[ti, zi, yi, xi + 1] * c2
1002+
V0 = self.V.data[ti, zi, yi, xi] * c1
1003+
V1 = self.V.data[ti, zi, yi + 1, xi] * c3
1004+
U = (1 - xsi) * U0 + xsi * U1
1005+
V = (1 - eta) * V0 + eta * V1
1006+
rad = np.pi / 180.0
1007+
deg2m = 1852 * 60.0
1008+
if applyConversion:
1009+
meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1
1010+
else:
1011+
meshJac = deg2m if grid.mesh == "spherical" else 1
1012+
1013+
jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
1014+
1015+
u = (
1016+
(-(1 - eta) * U - (1 - xsi) * V) * px[0]
1017+
+ ((1 - eta) * U - xsi * V) * px[1]
1018+
+ (eta * U + xsi * V) * px[2]
1019+
+ (-eta * U + (1 - xsi) * V) * px[3]
1020+
) / jac
1021+
v = (
1022+
(-(1 - eta) * U - (1 - xsi) * V) * py[0]
1023+
+ ((1 - eta) * U - xsi * V) * py[1]
1024+
+ (eta * U + xsi * V) * py[2]
1025+
+ (-eta * U + (1 - xsi) * V) * py[3]
1026+
) / jac
1027+
if isinstance(u, da.core.Array):
1028+
u = u.compute()
1029+
v = v.compute()
1030+
return (u, v)
1031+
1032+
u, v = _calc_UV(ti, yi, xi)
1033+
if ti < self.U.grid.tdim - 1 and time > self.U.grid.time[ti]:
1034+
ut1, vt1 = _calc_UV(ti + 1, yi, xi)
1035+
u = (1 - tau) * u + tau * ut1
1036+
v = (1 - tau) * v + tau * vt1
10281037
return (u, v)
10291038

10301039
def spatial_c_grid_interpolation3D_full(self, ti, z, y, x, time, particle=None):
@@ -1255,7 +1264,7 @@ def spatial_c_grid_interpolation3D(self, ti, z, y, x, time, particle=None, apply
12551264
else:
12561265
if self.gridindexingtype == "croco":
12571266
z = _croco_from_z_to_sigma_scipy(self.fieldset, time, z, y, x, particle=particle)
1258-
(u, v) = self.spatial_c_grid_interpolation2D(ti, z, y, x, time, particle=particle)
1267+
(u, v) = self.spatial_c_grid_interpolation2D(time, z, y, x, particle=particle)
12591268
w = self.W.eval(time, z, y, x, particle=particle, applyConversion=False)
12601269
if applyConversion:
12611270
w = self.W.units.to_target(w, z, y, x)
@@ -1387,14 +1396,8 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True):
13871396
elif self.U.interp_method == "cgrid_velocity":
13881397
tau, ti = self.U._time_index(time)
13891398
(u, v) = self.spatial_c_grid_interpolation2D(
1390-
ti, z, y, x, time, particle=particle, applyConversion=applyConversion
1399+
time, z, y, x, particle=particle, applyConversion=applyConversion
13911400
)
1392-
if ti < self.U.grid.tdim - 1 and time > self.U.grid.time[ti]:
1393-
(u1, v1) = self.spatial_c_grid_interpolation2D(
1394-
ti + 1, z, y, x, time, particle=particle, applyConversion=applyConversion
1395-
)
1396-
u = u * (1 - tau) + u1 * tau
1397-
v = v * (1 - tau) + v1 * tau
13981401
if "3D" in self.vector_type:
13991402
w = self.W.eval(time, z, y, x, particle=particle, applyConversion=applyConversion)
14001403
return (u, v, w)

0 commit comments

Comments
 (0)