Argo floats in shallow waters

src/virtualship/instruments/argo_float.py

@@ -53,6 +53,7 @@ class ArgoFloat:
         Variable("vertical_speed", dtype=np.float32),
         Variable("cycle_days", dtype=np.int32),
         Variable("drift_days", dtype=np.int32),
+        Variable("grounded", dtype=np.int32, initial=0),
     ]
 )
 
@@ -64,10 +65,23 @@ class ArgoFloat:
 def _argo_float_vertical_movement(particle, fieldset, time):
     if particle.cycle_phase == 0:
         # Phase 0: Sinking with vertical_speed until depth is drift_depth
-        particle_ddepth += (  # noqa Parcels defines particle_* variables, which code checkers cannot know.
+        particle_ddepth += (  # noqa
             particle.vertical_speed * particle.dt
         )
-        if particle.depth + particle_ddepth <= particle.drift_depth:
+
+        # bathymetry at particle location
+        loc_bathy = fieldset.bathymetry.eval(
+            time, particle.depth, particle.lat, particle.lon
+        )
+        if particle.depth + particle_ddepth <= loc_bathy:
+            particle_ddepth = loc_bathy - particle.depth + 50.0  # 50m above bathy
+            particle.cycle_phase = 1
+            particle.grounded = 1
+            print(
+                "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to drift depth. Raising by 50m above bathymetry and continuing cycle."
+            )
+
+        elif particle.depth + particle_ddepth <= particle.drift_depth:
             particle_ddepth = particle.drift_depth - particle.depth
             particle.cycle_phase = 1
 
@@ -81,7 +95,17 @@ def _argo_float_vertical_movement(particle, fieldset, time):
     elif particle.cycle_phase == 2:
         # Phase 2: Sinking further to max_depth
         particle_ddepth += particle.vertical_speed * particle.dt
-        if particle.depth + particle_ddepth <= particle.max_depth:
+        loc_bathy = fieldset.bathymetry.eval(
+            time, particle.depth, particle.lat, particle.lon
+        )
+        if particle.depth + particle_ddepth <= loc_bathy:
+            particle_ddepth = loc_bathy - particle.depth + 50.0  # 50m above bathy
+            particle.cycle_phase = 3
+            particle.grounded = 1
+            print(
+                "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to max depth. Raising by 50m above bathymetry and continuing cycle."
+            )
+        elif particle.depth + particle_ddepth <= particle.max_depth:
             particle_ddepth = particle.max_depth - particle.depth
             particle.cycle_phase = 3
 
@@ -91,6 +115,7 @@ def _argo_float_vertical_movement(particle, fieldset, time):
         particle.cycle_age += (
             particle.dt
         )  # solve issue of not updating cycle_age during ascent
+        particle.grounded = 0
         if particle.depth + particle_ddepth >= particle.min_depth:
             particle_ddepth = particle.min_depth - particle.depth
             particle.temperature = (
@@ -148,7 +173,7 @@ def __init__(self, expedition, from_data):
         super().__init__(
             expedition,
             variables,
-            add_bathymetry=False,
+            add_bathymetry=True,
             allow_time_extrapolation=False,
             verbose_progress=True,
             spacetime_buffer_size=spacetime_buffer_size,
@@ -171,6 +196,22 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
+        shallow_waypoints = {}
+        for i, m in enumerate(measurements):
+            loc_bathy = fieldset.bathymetry.eval(
+                time=0,
+                z=0,
+                y=m.spacetime.location.lat,
+                x=m.spacetime.location.lon,
+            )
+            if abs(loc_bathy) < 50.0:
+                shallow_waypoints[f"Waypoint {i + 1}"] = f"{abs(loc_bathy):.2f}m depth"
+
+        if len(shallow_waypoints) > 0:
+            raise ValueError(
+                f"{self.__class__.__name__} cannot be deployed in waters shallower than 50m. The following waypoints are too shallow: {shallow_waypoints}."
+            )
+
         # define parcel particles
         argo_float_particleset = ParticleSet(
             fieldset=fieldset,
@@ -191,7 +232,7 @@ def simulate(self, measurements, out_path) -> None:
         out_file = argo_float_particleset.ParticleFile(
             name=out_path,
             outputdt=OUTPUT_DT,
-            chunks=[len(argo_float_particleset), 100],
+            chunks=[len(argo_float_particleset), 100],  # TODO: is this necessary?
         )
 
         # get earliest between fieldset end time and provide end time

tests/instruments/test_argo_float.py

@@ -27,6 +27,7 @@ def test_simulate_argo_floats(tmpdir) -> None:
     u = np.full((2, 2, 2), 1.0)
     t = np.full((2, 2, 2), CONST_TEMPERATURE)
     s = np.full((2, 2, 2), CONST_SALINITY)
+    bathy = np.full((2, 2), -5000.0)
 
     fieldset = FieldSet.from_data(
         {"V": v, "U": u, "T": t, "S": s},
@@ -39,6 +40,15 @@ def test_simulate_argo_floats(tmpdir) -> None:
             ],
         },
     )
+    fieldset.add_field(
+        FieldSet.from_data(
+            {"bathymetry": bathy},
+            {
+                "lon": np.array([0.0, 10.0]),
+                "lat": np.array([0.0, 10.0]),
+            },
+        ).bathymetry
+    )
 
     # argo floats to deploy
     argo_floats = [