Harmonise Field dimensions order

parcels/application_kernels/advection.py

@@ -184,8 +184,8 @@
         time_i = np.linspace(0, fieldset.U.grid.time[ti + 1] - fieldset.U.grid.time[ti], I_s)
         ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])
 
-    xsi, eta, zeta, xi, yi, zi = fieldset.U._search_indices(
-        particle.lon, particle.lat, particle.depth, particle=particle
+    zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices(
+        -1, particle.depth, particle.lat, particle.lon, particle=particle
     )
     if withW:
         if abs(xsi - 1) < tol:
@@ -232,14 +232,14 @@
     else:
         dz = 1.0
 
-    c1 = fieldset.UV.dist(px[0], px[1], py[0], py[1], grid.mesh, np.dot(i_u.phi2D_lin(xsi, 0.0), py))
-    c2 = fieldset.UV.dist(px[1], px[2], py[1], py[2], grid.mesh, np.dot(i_u.phi2D_lin(1.0, eta), py))
-    c3 = fieldset.UV.dist(px[2], px[3], py[2], py[3], grid.mesh, np.dot(i_u.phi2D_lin(xsi, 1.0), py))
-    c4 = fieldset.UV.dist(px[3], px[0], py[3], py[0], grid.mesh, np.dot(i_u.phi2D_lin(0.0, eta), py))
+    c1 = fieldset.UV.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
+    c2 = fieldset.UV.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
+    c3 = fieldset.UV.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
+    c4 = fieldset.UV.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
     rad = np.pi / 180.0
     deg2m = 1852 * 60.0
     meshJac = (deg2m * deg2m * math.cos(rad * particle.lat)) if grid.mesh == "spherical" else 1
-    dxdy = fieldset.UV.jacobian(xsi, eta, px, py) * meshJac
+    dxdy = fieldset.UV.jacobian(py, px, eta, xsi) * meshJac
 
     if withW:
         U0 = direction * fieldset.U.data[ti, zi + 1, yi + 1, xi] * c4 * dz

parcels/compilation/codegenerator.py

@@ -834,7 +834,7 @@
             statements_croco = [
                 c.Statement(f"float cs_w[] = {*Cs_w, }".replace("(", "{").replace(")", "}")),
                 c.Statement(
-                    f"{node.var} = croco_from_z_to_sigma(U, H, Zeta, {args[3]}, {args[2]}, {args[1]}, time, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], hc, &cs_w)"
+                    f"{node.var} = croco_from_z_to_sigma(time, {args[1]}, {args[2]}, {args[3]}, U, H, Zeta, &particles->ti[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->xi[pnum*ngrid], hc, &cs_w)"
                 ),
             ]
             args = (args[0], node.var, args[2], args[3])
@@ -863,7 +863,7 @@
             statements_croco = [
                 c.Statement(f"float cs_w[] = {*Cs_w, }".replace("(", "{").replace(")", "}")),
                 c.Statement(
-                    f"{node.var4} = croco_from_z_to_sigma(U, H, Zeta, {args[3]}, {args[2]}, {args[1]}, time, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], hc, &cs_w)"
+                    f"{node.var4} = croco_from_z_to_sigma(time, {args[1]}, {args[2]}, {args[3]}, U, H, Zeta, &particles->ti[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->xi[pnum*ngrid], hc, &cs_w)"
                 ),
             ]
             args = (args[0], node.var4, args[2], args[3])

parcels/fieldfilebuffer.py

@@ -92,7 +92,7 @@ def parse_name(self, name):
         return name
 
     @property
-    def lonlat(self):
+    def latlon(self):
         lon = self.dataset[self.dimensions["lon"]]
         lat = self.dataset[self.dimensions["lat"]]
         if self.nolonlatindices and self.gridindexingtype not in ["croco"]:
@@ -141,7 +141,7 @@ def lonlat(self):
             if rectilinear:
                 lon_subset = lon_subset[0, :]
                 lat_subset = lat_subset[:, 0]
-        return lon_subset, lat_subset
+        return lat_subset, lon_subset
 
     @property
     def depth(self):

parcels/include/index_search.h

@@ -127,9 +127,10 @@ static inline StatusCode search_indices_vertical_z(type_coord z, int zdim, float
   return SUCCESS;
 }
 
-static inline StatusCode search_indices_vertical_s(type_coord z, int xdim, int ydim, int zdim, float *zvals,
-                                    int xi, int yi, int *zi, double xsi, double eta, double *zeta,
-                                    int z4d, int ti, int tdim, double time, double t0, double t1, int interp_method)
+static inline StatusCode search_indices_vertical_s(double time, type_coord z,
+                                                   int tdim, int zdim, int ydim, int xdim, float *zvals,
+                                                   int ti, int *zi, int yi, int xi, double *zeta, double eta, double xsi,
+                                                   int z4d, double t0, double t1, int interp_method)
 {
   if (interp_method == BGRID_VELOCITY || interp_method == BGRID_W_VELOCITY || interp_method == BGRID_TRACER){
     xsi = 1;
@@ -184,7 +185,7 @@ static inline StatusCode search_indices_vertical_s(type_coord z, int xdim, int y
   return SUCCESS;
 }
 
-static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, int onlyX, int sphere_mesh)
+static inline void reconnect_bnd_indices(int *yi, int *xi, int ydim, int xdim, int onlyX, int sphere_mesh)
 {
   if (*xi < 0){
     if (sphere_mesh)
@@ -211,10 +212,11 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i
 }
 
 
-static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype,
-                                                   int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
-                                                   int ti, double time, double t0, double t1, int interp_method,
-                                                   int gridindexingtype)
+static inline StatusCode search_indices_rectilinear(double time, type_coord z, type_coord y, type_coord x,
+                                                    CStructuredGrid *grid, GridType gtype,
+                                                    int ti, int *zi, int *yi, int *xi, double *zeta, double *eta, double *xsi,
+                                                    double t0, double t1, int interp_method,
+                                                    int gridindexingtype)
 {
   int xdim = grid->xdim;
   int ydim = grid->ydim;
@@ -261,7 +263,7 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
         ++(*xi);
       else if (xvalsi > x)
         --(*xi);
-      reconnect_bnd_indices(xi, yi, xdim, ydim, 1, 1);
+      reconnect_bnd_indices(yi, xi, ydim, xdim, 1, 1);
       xvalsi = xvals[*xi];
       if (xvalsi < x - 225) xvalsi += 360;
       if (xvalsi > x + 225) xvalsi -= 360;
@@ -294,9 +296,9 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
         status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
         break;
       case RECTILINEAR_S_GRID:
-        status = search_indices_vertical_s(z, xdim, ydim, zdim, zvals,
-                                        *xi, *yi, zi, *xsi, *eta, zeta,
-                                        z4d, ti, tdim, time, t0, t1, interp_method);
+        status = search_indices_vertical_s(time, z, tdim, zdim, ydim, xdim, zvals,
+                                           ti, zi, *yi, *xi, zeta, *eta, *xsi,
+                                           z4d, t0, t1, interp_method);
         break;
       default:
         status = ERRORINTERPOLATION;
@@ -321,10 +323,12 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
 }
 
 
-static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype,
-                                                   int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
-                                                   int ti, double time, double t0, double t1, int interp_method,
-                                                   int gridindexingtype)
+static inline StatusCode search_indices_curvilinear(double time, type_coord z, type_coord y, type_coord x,
+                                                    CStructuredGrid *grid, GridType gtype,
+                                                    int ti, int *zi, int *yi, int *xi,
+                                                    double *zeta, double *eta, double *xsi,
+                                                    double t0, double t1, int interp_method,
+                                                    int gridindexingtype)
 {
   int xi_old = *xi;
   int yi_old = *yi;
@@ -407,23 +411,23 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
       (*yi)--;
     if (*eta > 1+tol)
       (*yi)++;
-    reconnect_bnd_indices(xi, yi, xdim, ydim, 0, sphere_mesh);
+    reconnect_bnd_indices(yi, xi, ydim, xdim, 0, sphere_mesh);
     it++;
     if ( it > maxIterSearch){
-      printf("Correct cell not found for (%f, %f) after %d iterations\n", x, y, maxIterSearch);
+      printf("Correct cell not found for (lat, lon) = (%f, %f) after %d iterations\n", y, x, maxIterSearch);
       printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old);
       printf("            new particle indices: (yi, xi) %d %d\n", *yi, *xi);
       printf("            Mesh 2d shape:  %d %d\n", ydim, xdim);
-      printf("            Relative particle position:  (xsi, eta) %1.16e %1.16e\n", *xsi, *eta);
+      printf("            Relative particle position:  (eta, xsi) %1.16e %1.16e\n", *eta, *xsi);
       return ERROROUTOFBOUNDS;
     }
   }
   if ( (*xsi != *xsi) || (*eta != *eta) ){  // check if nan
-      printf("Correct cell not found for (%f, %f))\n", x, y);
+      printf("Correct cell not found for (lat, lon) = (%f, %f))\n", y, x);
       printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old);
       printf("            new particle indices: (yi, xi) %d %d\n", *yi, *xi);
       printf("            Mesh 2d shape:  %d %d\n", ydim, xdim);
-      printf("            Relative particle position:  (xsi, eta) %1.16e %1.16e\n", *xsi, *eta);
+      printf("            Relative particle position:  (eta, xsi) %1.16e %1.16e\n", *eta, *xsi);
       return ERROROUTOFBOUNDS;
   }
   if (*xsi < 0) *xsi = 0;
@@ -438,9 +442,9 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
         status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
         break;
       case CURVILINEAR_S_GRID:
-        status = search_indices_vertical_s(z, xdim, ydim, zdim, zvals,
-                                        *xi, *yi, zi, *xsi, *eta, zeta,
-                                        z4d, ti, tdim, time, t0, t1, interp_method);
+        status = search_indices_vertical_s(time, z, tdim, ydim, xdim, zdim, zvals,
+                                           ti, zi, *yi, *xi, zeta, *eta, *xsi,
+                                           z4d, t0, t1, interp_method);
         break;
       default:
         status = ERRORINTERPOLATION;
@@ -460,21 +464,23 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
 /* Local linear search to update grid index
  * params ti, sizeT, time. t0, t1 are only used for 4D S grids
  * */
-static inline StatusCode search_indices(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid,
-                                       int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
-                                       GridType gtype, int ti, double time, double t0, double t1, int interp_method,
-                                       int gridindexingtype)
+static inline StatusCode search_indices(double time, type_coord z, type_coord y, type_coord x,
+                                        CStructuredGrid *grid,
+                                        int ti, int *zi, int *yi, int *xi,
+                                        double *zeta, double *eta, double *xsi,
+                                        GridType gtype, double t0, double t1, int interp_method,
+                                        int gridindexingtype)
 {
   switch(gtype){
     case RECTILINEAR_Z_GRID:
     case RECTILINEAR_S_GRID:
-      return search_indices_rectilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta,
-                                   ti, time, t0, t1, interp_method, gridindexingtype);
+      return search_indices_rectilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi,
+                                        t0, t1, interp_method, gridindexingtype);
       break;
     case CURVILINEAR_Z_GRID:
     case CURVILINEAR_S_GRID:
-      return search_indices_curvilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta,
-                                   ti, time, t0, t1, interp_method, gridindexingtype);
+      return search_indices_curvilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi,
+                                        t0, t1, interp_method, gridindexingtype);
       break;
     default:
       printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n");

parcels/include/interpolation_utils.h

@@ -8,7 +8,7 @@ typedef enum
   } Orientation;
 
 
-static inline void phi2D_lin(double xsi, double eta, double *phi)
+static inline void phi2D_lin(double eta, double xsi, double *phi)
 {
     phi[0] = (1-xsi) * (1-eta);
     phi[1] =    xsi  * (1-eta);
@@ -25,40 +25,40 @@ static inline void phi1D_quad(double xsi, double *phi)
 }
 
 
-static inline void dphidxsi3D_lin(double xsi, double eta, double zet, double *dphidxsi, double *dphideta, double *dphidzet)
+static inline void dphidxsi3D_lin(double zeta, double eta, double xsi, double *dphidzeta, double *dphideta, double *dphidxsi)
 {
-  dphidxsi[0] = - (1-eta) * (1-zet);
-  dphidxsi[1] =   (1-eta) * (1-zet);
-  dphidxsi[2] =   (  eta) * (1-zet);
-  dphidxsi[3] = - (  eta) * (1-zet);
-  dphidxsi[4] = - (1-eta) * (  zet);
-  dphidxsi[5] =   (1-eta) * (  zet);
-  dphidxsi[6] =   (  eta) * (  zet);
-  dphidxsi[7] = - (  eta) * (  zet);
-
-  dphideta[0] = - (1-xsi) * (1-zet);
-  dphideta[1] = - (  xsi) * (1-zet);
-  dphideta[2] =   (  xsi) * (1-zet);
-  dphideta[3] =   (1-xsi) * (1-zet);
-  dphideta[4] = - (1-xsi) * (  zet);
-  dphideta[5] = - (  xsi) * (  zet);
-  dphideta[6] =   (  xsi) * (  zet);
-  dphideta[7] =   (1-xsi) * (  zet);
-
-  dphidzet[0] = - (1-xsi) * (1-eta);
-  dphidzet[1] = - (  xsi) * (1-eta);
-  dphidzet[2] = - (  xsi) * (  eta);
-  dphidzet[3] = - (1-xsi) * (  eta);
-  dphidzet[4] =   (1-xsi) * (1-eta);
-  dphidzet[5] =   (  xsi) * (1-eta);
-  dphidzet[6] =   (  xsi) * (  eta);
-  dphidzet[7] =   (1-xsi) * (  eta);
+  dphidxsi[0] = - (1-eta) * (1-zeta);
+  dphidxsi[1] =   (1-eta) * (1-zeta);
+  dphidxsi[2] =   (  eta) * (1-zeta);
+  dphidxsi[3] = - (  eta) * (1-zeta);
+  dphidxsi[4] = - (1-eta) * (  zeta);
+  dphidxsi[5] =   (1-eta) * (  zeta);
+  dphidxsi[6] =   (  eta) * (  zeta);
+  dphidxsi[7] = - (  eta) * (  zeta);
+
+  dphideta[0] = - (1-xsi) * (1-zeta);
+  dphideta[1] = - (  xsi) * (1-zeta);
+  dphideta[2] =   (  xsi) * (1-zeta);
+  dphideta[3] =   (1-xsi) * (1-zeta);
+  dphideta[4] = - (1-xsi) * (  zeta);
+  dphideta[5] = - (  xsi) * (  zeta);
+  dphideta[6] =   (  xsi) * (  zeta);
+  dphideta[7] =   (1-xsi) * (  zeta);
+
+  dphidzeta[0] = - (1-xsi) * (1-eta);
+  dphidzeta[1] = - (  xsi) * (1-eta);
+  dphidzeta[2] = - (  xsi) * (  eta);
+  dphidzeta[3] = - (1-xsi) * (  eta);
+  dphidzeta[4] =   (1-xsi) * (1-eta);
+  dphidzeta[5] =   (  xsi) * (1-eta);
+  dphidzeta[6] =   (  xsi) * (  eta);
+  dphidzeta[7] =   (1-xsi) * (  eta);
 }
 
-static inline void dxdxsi3D_lin(double *px, double *py, double *pz, double xsi, double eta, double zet, double *jacM, int sphere_mesh)
+static inline void dxdxsi3D_lin(double *pz, double *py, double *px, double zeta, double eta, double xsi, double *jacM, int sphere_mesh)
 {
-  double dphidxsi[8], dphideta[8], dphidzet[8];
-  dphidxsi3D_lin(xsi, eta, zet, dphidxsi, dphideta, dphidzet);
+  double dphidxsi[8], dphideta[8], dphidzeta[8];
+  dphidxsi3D_lin(zeta, eta, xsi, dphidzeta, dphideta, dphidxsi);
 
   int i;
   for(i=0; i<9; ++i)
@@ -76,20 +76,22 @@ static inline void dxdxsi3D_lin(double *px, double *py, double *pz, double xsi,
   for(i=0; i<8; ++i){
     jacM[3*0+0] += px[i] * dphidxsi[i] * jac_lon; // dxdxsi
     jacM[3*0+1] += px[i] * dphideta[i] * jac_lon; // dxdeta
-    jacM[3*0+2] += px[i] * dphidzet[i] * jac_lon; // dxdzet
+    jacM[3*0+2] += px[i] * dphidzeta[i] * jac_lon; // dxdzeta
     jacM[3*1+0] += py[i] * dphidxsi[i] * jac_lat; // dydxsi
     jacM[3*1+1] += py[i] * dphideta[i] * jac_lat; // dydeta
-    jacM[3*1+2] += py[i] * dphidzet[i] * jac_lat; // dydzet
+    jacM[3*1+2] += py[i] * dphidzeta[i] * jac_lat; // dydzeta
     jacM[3*2+0] += pz[i] * dphidxsi[i];           // dzdxsi
     jacM[3*2+1] += pz[i] * dphideta[i];           // dzdeta
-    jacM[3*2+2] += pz[i] * dphidzet[i];           // dzdzet
+    jacM[3*2+2] += pz[i] * dphidzeta[i];           // dzdzeta
   }
 }
 
-static inline double jacobian3D_lin_face(double *px, double *py, double *pz, double xsi, double eta, double zet, Orientation orientation, int sphere_mesh)
+static inline double jacobian3D_lin_face(double *pz, double *py, double *px,
+                                         double zeta, double eta, double xsi,
+                                         Orientation orientation, int sphere_mesh)
 {
   double jacM[9];
-  dxdxsi3D_lin(px, py, pz, xsi, eta, zet, jacM, sphere_mesh);
+  dxdxsi3D_lin(pz, py, px, zeta, eta, xsi, jacM, sphere_mesh);
 
   double j[3];
 
@@ -112,10 +114,12 @@ static inline double jacobian3D_lin_face(double *px, double *py, double *pz, dou
   return sqrt(j[0]*j[0]+j[1]*j[1]+j[2]*j[2]);
 }
 
-static inline double jacobian3D_lin(double *px, double *py, double *pz, double xsi, double eta, double zet, int sphere_mesh)
+static inline double jacobian3D_lin(double *pz, double *py, double *px,
+                                    double zeta, double eta, double xsi,
+                                    int sphere_mesh)
 {
   double jacM[9];
-  dxdxsi3D_lin(px, py, pz, xsi, eta, zet, jacM, sphere_mesh);
+  dxdxsi3D_lin(pz, py, px, zeta, eta, xsi, jacM, sphere_mesh);
 
   double jac = jacM[3*0+0] * (jacM[3*1+1]*jacM[3*2+2] - jacM[3*2+1]*jacM[3*1+2])
              - jacM[3*0+1] * (jacM[3*1+0]*jacM[3*2+2] - jacM[3*2+0]*jacM[3*1+2])