diff --git a/doc/source/documentation/cookbook/beaching.ipynb b/doc/source/documentation/cookbook/beaching.ipynb new file mode 100644 index 0000000..ec488e2 --- /dev/null +++ b/doc/source/documentation/cookbook/beaching.ipynb @@ -0,0 +1,55 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a4abeb9c-4943-4b50-8b16-3a6a824b0d2f", + "metadata": {}, + "source": [ + "# Beaching\n", + "\n", + "In the ocean, objects that are transported by ocean currents may become beached or stranded on land. If an object becomes beached in an intertidal zone, the beaching event may be short-lived, with the object being resuspended on the next flood tide. In contrast, if the object is deposited on the strand line during a spring tide, it may be several days or weeks before it has an opportunity to be resuspended. In addition, beaching will also be influenced by waves, the wind and storm surges. Objects may even become permanently or semi-permanently beached if deposited far inland during a storm surge. Beaching and resuspension dynamics are also dependent on the type of object under consideration.\n", + "\n", + "Particle tracking models handle beaching in different ways. The approach taken is strongly influenced by the type of driving data that is being used. Driving data from a global or broad scale ocean circulation model will not, in general, include wetting and drying with changes in the tide - the spatial resolution is too poor - meaning grid cells remain wet at all times. However, beaching may still be parameterised; for example, by imposing a probability that a particle becomes beached if it resides within a given distance of the model's land boundary, or if it crosses a land boundary during the simulation. Resuspension may be handled in a similar way. Alternatively, beaching may be ignored completely. In this scenario, if a simulated particle does cross a land boundary - which can happen for various reasons - it may be restored to its previous position or reflected back into the model domain. When using the output of fine scale models of coastal or estuarine systems which include wetting and drying, one is also forced to consider what happens to particles that become trapped in intertidal zones.\n", + "\n", + "## Beaching in PyLag\n", + "\n", + "Presently, PyLag does not include a paramaterisation of beaching in the absence of wetting and drying; if a land boundary is crossed, particles are either restored to their previous position or reflected back into the model domain depending on the type of boundary condition calculator used. However, it does include a parameterisation for wetting and drying. There are two options which are controlled through the binary flag `allow_beaching`, which should be set in the `[SIMULATION]` section of the PyLag run configuration file.\n", + "\n", + "If `allow_beaching=True`, then particles can become beached in intertidal zones. If this occurs, they are flagged as having beached using the particle attribute `is_beached`. Particles that are beached are passed over until the host cell they are in is flooded with water once more. If `allow_beaching=False`, particles are prevented from beaching by moving them to the centroid of the nearest wet element.\n", + "\n", + "Although in reality intertidal zones are free of water, this is not the case in models. Typically, models will impose a threshold water depth at which a cell is deemed \"dry\". If a particle is in an element when this threshold is met, it will become beached. The situation is complicated slightly when doing offline particle tracking simulations using discrete model outputs which are (necessarily!) saved infrequently relative to the model time step. In this event, the tidal cycle may be poorly resolved - or even aliased and the period for which the particle remains beached may become artificially long.\n", + "\n", + "When beaching is excluded, and the run is fully deterministic (i.e., diffusion is neglected), care should also be taken when particles beach in the same location at the same time. If this occurs, the selected particles will be moved to the same point in space - the centroid of the nearest wet element - and continue to move as one." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecbed4b4-e69b-4f38-ba7e-06a3aa138fef", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/doc/source/documentation/cookbook/index.rst b/doc/source/documentation/cookbook/index.rst index 353765e..7377827 100644 --- a/doc/source/documentation/cookbook/index.rst +++ b/doc/source/documentation/cookbook/index.rst @@ -11,6 +11,7 @@ In the cookbook, you can find hints and tips for working with PyLag; explanation initial_positions boundary_conditions + beaching fixed_vertical_positioning boundary_layer_dynamics settling diff --git a/pylag/data_reader.pxd b/pylag/data_reader.pxd index 97932e2..d5731b4 100644 --- a/pylag/data_reader.pxd +++ b/pylag/data_reader.pxd @@ -33,6 +33,8 @@ cdef class DataReader: DTYPE_FLOAT_t particle_pathline[2], DTYPE_FLOAT_t intersection[2]) + cdef find_nearest_wet_host(self, DTYPE_FLOAT_t time, Particle *particle) + cdef set_default_location(self, Particle *particle) cdef set_local_coordinates(self, Particle *particle) diff --git a/pylag/data_reader.pyx b/pylag/data_reader.pyx index c9f8f1c..42d11c7 100644 --- a/pylag/data_reader.pyx +++ b/pylag/data_reader.pyx @@ -178,6 +178,29 @@ cdef class DataReader: DTYPE_FLOAT_t intersection[2]): raise NotImplementedError + def find_nearest_wet_host_wrapper(self, DTYPE_FLOAT_t time, + ParticleSmartPtr particle): + """ Python wrapper for finding the nearest wet host element + + This function is used to find the nearest wet host element to the + particle's current position. + + Parameters + ---------- + time : float + The current time. + + particle : pylag.particle_cpp_wrapper.ParticleSmartPtr + The particle at its current position. + + Returns + ------- + """ + return self.find_nearest_wet_host(time, particle.get_ptr()) + + cdef find_nearest_wet_host(self, DTYPE_FLOAT_t time, Particle *particle): + raise NotImplementedError + def set_default_location_wrapper(self, ParticleSmartPtr particle): """ Python wrapper for setting the default location of a particle diff --git a/pylag/fvcom_data_reader.pyx b/pylag/fvcom_data_reader.pyx index 46c02ee..179b361 100644 --- a/pylag/fvcom_data_reader.pyx +++ b/pylag/fvcom_data_reader.pyx @@ -28,7 +28,7 @@ from pylag.data_types_cython cimport DTYPE_INT_t, DTYPE_FLOAT_t from libcpp.string cimport string from libcpp.vector cimport vector -from libc.math cimport log +from libc.math cimport log, sqrt # PyLag cython imports from pylag.parameters cimport cartesian, geographic, deg_to_radians @@ -463,6 +463,98 @@ cdef class FVCOMDataReader(DataReader): return + cdef find_nearest_wet_host(self, DTYPE_FLOAT_t time, Particle *particle): + """ Find the nearest wet host element + + Find the nearest wet host element to the particle's position and relocate + the particle to the centroid of this element. Begin by searching the + the particle's host element and then its immediate neighbours. If no wet + elements are found, the neighbours of neighbours are searched and so on. + If multiple wet elements are found, the one closest to the particle is + selected. + """ + cdef vector[DTYPE_INT_t] current_nbe_checks + cdef vector[DTYPE_INT_t] next_nbe_checks + cdef vector[DTYPE_INT_t] elems_checked + cdef vector[DTYPE_INT_t] wet_elems + + cdef DTYPE_FLOAT_t dist, min_dist + cdef DTYPE_INT_t current_host, nbe, n_wet_elems, new_host + cdef DTYPE_INT_t i, j + cdef bint has_been_checked + + # Confirm the particle currently resides in a dry element + if self.is_wet(time, particle) == 1: + return + + # Save the current host + current_host = particle.get_host_horizontal_elem(self._name) + + # Add the current host to the list of elements whose neighbours we will check next + next_nbe_checks.push_back(current_host) + + # Flag that we have already checked the element itself + elems_checked.push_back(current_host) + + # Perform an ever widening search to find the nearest wet element + while True: + current_nbe_checks = next_nbe_checks + + # Clear the list of elements to check next - we will repopulate this as we go + next_nbe_checks.clear() + + # Loop over the current list of elements whose neighbours we will check + for elem in current_nbe_checks: + + # Search the element's immediate neighbours + for i in range(3): + nbe = self._nbe[i, elem] + + # Check if sea point + if nbe == LAND_BDY_CROSSED or nbe == OPEN_BDY_CROSSED: + continue + + # Have we cheked this element before? + has_been_checked = False + for j in range(elems_checked.size()): + if nbe == elems_checked[j]: + has_been_checked = True + break + + if has_been_checked: + continue + else: + # Check if the neighbour is wet + particle.set_host_horizontal_elem(self._name, nbe) + + if self.is_wet(time, particle) == 1: + wet_elems.push_back(nbe) + + next_nbe_checks.push_back(nbe) + elems_checked.push_back(nbe) + + n_wet_elems = wet_elems.size() + if n_wet_elems == 0: + continue + else: + if n_wet_elems == 1: + new_host = wet_elems[0] + else: + # If multiple wet elements are found, select the one closest to the particle + min_dist = 1.e20 + for i in range(wet_elems.size()): + # TODO - Will work with geographic coordinates? + dist = sqrt((particle.get_x1() - self._xc[wet_elems[i]])**2 + + (particle.get_x2() - self._yc[wet_elems[i]])**2) + if dist < min_dist: + min_dist = dist + new_host = wet_elems[i] + + particle.set_host_horizontal_elem(self._name, new_host) + self.set_default_location(particle) + + return + cdef set_local_coordinates(self, Particle *particle): """ Set local coordinates diff --git a/pylag/numerics.pyx b/pylag/numerics.pyx index 8a830ff..5eaa7c3 100644 --- a/pylag/numerics.pyx +++ b/pylag/numerics.pyx @@ -107,6 +107,8 @@ cdef class StdNumMethod(NumMethod): """ cdef DTYPE_FLOAT_t _time_step + cdef bint _allow_beaching + cdef ItMethod _iterative_method cdef HorizBoundaryConditionCalculator _horiz_bc_calculator cdef VertBoundaryConditionCalculator _vert_bc_calculator @@ -123,6 +125,11 @@ cdef class StdNumMethod(NumMethod): self._time_step = self._iterative_method.get_time_step() + try: + self._allow_beaching = config.getboolean('SIMULATION', 'allow_beaching') + except (configparser.NoOptionError) as e: + self._allow_beaching = True + cdef DTYPE_INT_t step(self, DataReader data_reader, DTYPE_FLOAT_t time, Particle *particle) except INT_ERR: """ Perform one iteration of the numerical method @@ -176,6 +183,15 @@ cdef class StdNumMethod(NumMethod): # If the cell is still dry, pass over return IN_DOMAIN else: + # Apply depth/height restoring (the position of the free surface may have changed) + if _depth_restoring is True: + zmax = data_reader.get_zmax(time, &_particle_copy) + _particle_copy.set_x3(_fixed_depth_below_surface + zmax) + + elif _height_restoring is True: + zmin = data_reader.get_zmin(time, &_particle_copy) + _particle_copy.set_x3(_fixed_height_above_bed + zmin) + # Set vertical grid vars, which may have changed while the particle was beached flag = data_reader.set_vertical_grid_vars(time, &_particle_copy) @@ -210,69 +226,46 @@ cdef class StdNumMethod(NumMethod): if flag != IN_DOMAIN: return flag - # Restore to a fixed depth or height? - if _depth_restoring is True: - zmax = data_reader.get_zmax(time+self._time_step, &_particle_copy) - _particle_copy.set_x3(_fixed_depth_below_surface + zmax) + # Check for beaching + if data_reader.is_wet(time+self._time_step, &_particle_copy) == 0: - # Only try to set vertical grid vars if the particle is not beached - if data_reader.is_wet(time+self._time_step, &_particle_copy) == 1: + if self._allow_beaching: + _particle_copy.set_is_beached(1) - # Determine the new host zlayer - flag = data_reader.set_vertical_grid_vars(time+self._time_step, &_particle_copy) + # Copy back particle properties + particle[0] = _particle_copy - # Return if failure recorded - if flag != IN_DOMAIN: - return flag + return flag else: - _particle_copy.set_is_beached(1) + # Relocate to the nearest wet host. NB used t + dt! + data_reader.find_nearest_wet_host(time+self._time_step, &_particle_copy) - # Copy back particle properties - particle[0] = _particle_copy + # Restore to a fixed depth or height? + if _depth_restoring is True: + zmax = data_reader.get_zmax(time+self._time_step, &_particle_copy) + _particle_copy.set_x3(_fixed_depth_below_surface + zmax) - return flag - elif _height_restoring is True: zmin = data_reader.get_zmin(time+self._time_step, &_particle_copy) _particle_copy.set_x3(_fixed_height_above_bed + zmin) - # Only try to set vertical grid vars if the particle is not beached - if data_reader.is_wet(time+self._time_step, &_particle_copy) == 1: + # Determine the new host zlayer + flag = data_reader.set_vertical_grid_vars(time+self._time_step, &_particle_copy) - # Determine the new host zlayer - flag = data_reader.set_vertical_grid_vars(time+self._time_step, &_particle_copy) - - # Return if failure recorded - if flag != IN_DOMAIN: - return flag - else: - _particle_copy.set_is_beached(1) - - # Copy back particle properties - particle[0] = _particle_copy - - return flag - - # Check to see if the particle has beached. Only set vertical grid vars if it is in water. - if data_reader.is_wet(time+self._time_step, &_particle_copy) == 1: - # Set vertical grid vars. NB use t + dt! - flag = data_reader.set_vertical_grid_vars(time+self._time_step, &_particle_copy) - - # Apply surface/bottom boundary conditions if required + if flag != IN_DOMAIN: + # TODO - Depth/height restoring may fail in shallow water. Apply BC or return an error? + flag = self._vert_bc_calculator.apply(data_reader, time+self._time_step, &_particle_copy) + + # Return without updating the particle's position if failure recorded if flag != IN_DOMAIN: - flag = self._vert_bc_calculator.apply(data_reader, time+self._time_step, &_particle_copy) - - # Return if failure recorded - if flag != IN_DOMAIN: - return flag - else: - _particle_copy.set_is_beached(1) - + return flag + # Copy back particle properties particle[0] = _particle_copy return flag + cdef class OS0NumMethod(NumMethod): """ Numerical method that employs operator splitting @@ -317,6 +310,8 @@ cdef class OS0NumMethod(NumMethod): cdef DTYPE_FLOAT_t _diff_time_step cdef DTYPE_INT_t _n_sub_time_steps + cdef bint _allow_beaching + cdef ItMethod _adv_iterative_method cdef ItMethod _diff_iterative_method cdef HorizBoundaryConditionCalculator _horiz_bc_calculator @@ -351,6 +346,11 @@ cdef class OS0NumMethod(NumMethod): self._n_sub_time_steps = int(self._adv_time_step / self._diff_time_step) + try: + self._allow_beaching = config.getboolean('SIMULATION', 'allow_beaching') + except (configparser.NoOptionError) as e: + self._allow_beaching = True + cdef DTYPE_INT_t step(self, DataReader data_reader, DTYPE_FLOAT_t time, Particle *particle) except INT_ERR: """ Perform one iteration of the numerical integration @@ -410,6 +410,15 @@ cdef class OS0NumMethod(NumMethod): # If the cell is still dry, pass over return IN_DOMAIN else: + # Apply depth/height restoring (the position of the free surface may have changed) + if _depth_restoring is True: + zmax = data_reader.get_zmax(time, &_particle_copy_a) + _particle_copy_a.set_x3(_fixed_depth_below_surface + zmax) + + elif _height_restoring is True: + zmin = data_reader.get_zmin(time, &_particle_copy_a) + _particle_copy_a.set_x3(_fixed_height_above_bed + zmin) + # Set vertical grid vars, which may have changed while the particle was beached flag = data_reader.set_vertical_grid_vars(time, &_particle_copy_a) @@ -448,11 +457,15 @@ cdef class OS0NumMethod(NumMethod): # Check to see if the particle has beached. NB use time `time', # since this is when the diffusion loop starts. if data_reader.is_wet(time, &_particle_copy_a) == 0: - _particle_copy_a.set_is_beached(1) + if self._allow_beaching: + _particle_copy_a.set_is_beached(1) - particle[0] = _particle_copy_a + particle[0] = _particle_copy_a - return flag + return flag + else: + # Relocate to the nearest wet host + data_reader.find_nearest_wet_host(time, &_particle_copy_a) # Set vertical grid vars flag = data_reader.set_vertical_grid_vars(time, &_particle_copy_a) @@ -497,11 +510,15 @@ cdef class OS0NumMethod(NumMethod): # Check to see if the particle has beached if data_reader.is_wet(t+self._diff_time_step, &_particle_copy_b) == 0: - _particle_copy_b.set_is_beached(1) + if self._allow_beaching: + _particle_copy_b.set_is_beached(1) - particle[0] = _particle_copy_b + particle[0] = _particle_copy_b + + return flag + else: + data_reader.find_nearest_wet_host(t+self._diff_time_step, &_particle_copy_b) - return flag flag = data_reader.set_vertical_grid_vars(t+self._diff_time_step, &_particle_copy_b) @@ -516,57 +533,42 @@ cdef class OS0NumMethod(NumMethod): # Save the particle's last position to help with host element searching _particle_copy_a = _particle_copy_b - # Restore to a fixed depth or height? - if _depth_restoring is True: - zmax = data_reader.get_zmax(time+self._adv_time_step, &_particle_copy_b) - _particle_copy_b.set_x3(_fixed_depth_below_surface + zmax) - - # Only try to set vertical grid vars if the particle is not beached - if data_reader.is_wet(time+self._adv_time_step, &_particle_copy_b) == 1: - - # Determine the new host zlayer - flag = data_reader.set_vertical_grid_vars(time+self._adv_time_step, &_particle_copy_b) + # Check to see if the particle has beached + if data_reader.is_wet(time+self._adv_time_step, &_particle_copy_b) == 0: + if self._allow_beaching: + _particle_copy_b.set_is_beached(1) - # Return if failure recorded - if flag != IN_DOMAIN: - return flag + particle[0] = _particle_copy_b + return flag else: - _particle_copy_b.set_is_beached(1) - - particle[0] = _particle_copy_b + data_reader.find_nearest_wet_host(time+self._adv_time_step, &_particle_copy_b) - return flag + # Restore to a fixed depth or height? + if _depth_restoring is True: + zmax = data_reader.get_zmax(time+self._adv_time_step, &_particle_copy_b) + _particle_copy_b.set_x3(_fixed_depth_below_surface + zmax) elif _height_restoring is True: zmin = data_reader.get_zmin(time+self._adv_time_step, &_particle_copy_b) _particle_copy_b.set_x3(_fixed_height_above_bed + zmin) - # Only try to set vertical grid vars if the particle is not beached - if data_reader.is_wet(time+self._adv_time_step, &_particle_copy_b) == 1: - - # Determine the new host zlayer - flag = data_reader.set_vertical_grid_vars(time+self._adv_time_step, &_particle_copy_b) + # Determine the new host zlayer + flag = data_reader.set_vertical_grid_vars(time+self._adv_time_step, &_particle_copy_b) - # Return if failure recorded - if flag != IN_DOMAIN: - return flag - - else: - _particle_copy_b.set_is_beached(1) - - particle[0] = _particle_copy_b - - return flag + if flag != IN_DOMAIN: + # TODO - Depth/height restoring may fail in shallow water. Apply BC or return an error? + flag = self._vert_bc_calculator.apply(data_reader, t+self._diff_time_step, &_particle_copy_b) - # Check to see if the particle has beached - if data_reader.is_wet(time+self._adv_time_step, &_particle_copy_b) == 0: - _particle_copy_b.set_is_beached(1) + # Return if failure recorded + if flag != IN_DOMAIN: + return flag particle[0] = _particle_copy_b return flag + cdef class OS1NumMethod(NumMethod): """ Numerical method that employs strang splitting @@ -603,6 +605,8 @@ cdef class OS1NumMethod(NumMethod): cdef DTYPE_FLOAT_t _adv_time_step cdef DTYPE_FLOAT_t _diff_time_step + cdef bint _allow_beaching + cdef ItMethod _adv_iterative_method cdef ItMethod _diff_iterative_method cdef HorizBoundaryConditionCalculator _horiz_bc_calculator @@ -629,6 +633,11 @@ cdef class OS1NumMethod(NumMethod): "the numerical integration scheme OS1NumMethod."\ "".format(self._adv_time_step, self._diff_time_step)) + try: + self._allow_beaching = config.getboolean('SIMULATION', 'allow_beaching') + except (configparser.NoOptionError) as e: + self._allow_beaching = True + cdef DTYPE_INT_t step(self, DataReader data_reader, DTYPE_FLOAT_t time, Particle *particle) except INT_ERR: """ Perform one iteration of the numerical integration @@ -687,6 +696,15 @@ cdef class OS1NumMethod(NumMethod): # If the cell is still dry, pass over return IN_DOMAIN else: + # Apply depth/height restoring (the position of the free surface may have changed) + if _depth_restoring is True: + zmax = data_reader.get_zmax(time, &_particle_copy_a) + _particle_copy_a.set_x3(_fixed_depth_below_surface + zmax) + + elif _height_restoring is True: + zmin = data_reader.get_zmin(time, &_particle_copy_a) + _particle_copy_a.set_x3(_fixed_height_above_bed + zmin) + # Set vertical grid vars, which may have changed while the particle was beached flag = data_reader.set_vertical_grid_vars(time, &_particle_copy_a) @@ -723,11 +741,14 @@ cdef class OS1NumMethod(NumMethod): # Check is beached status. NB evaluated at time `time', since this is when the # advection update starts if data_reader.is_wet(time, &_particle_copy_a) == 0: - _particle_copy_a.set_is_beached(1) + if self._allow_beaching: + _particle_copy_a.set_is_beached(1) - particle[0] = _particle_copy_a + particle[0] = _particle_copy_a - return flag + return flag + else: + data_reader.find_nearest_wet_host(time, &_particle_copy_a) flag = data_reader.set_vertical_grid_vars(time, &_particle_copy_a) @@ -776,11 +797,14 @@ cdef class OS1NumMethod(NumMethod): # Check is beached status if data_reader.is_wet(t, &_particle_copy_b) == 0: - _particle_copy_a.set_is_beached(1) + if self._allow_beaching: + _particle_copy_a.set_is_beached(1) - particle[0] = _particle_copy_b + particle[0] = _particle_copy_b - return flag + return flag + else: + data_reader.find_nearest_wet_host(t, &_particle_copy_b) flag = data_reader.set_vertical_grid_vars(t, &_particle_copy_b) @@ -820,71 +844,36 @@ cdef class OS1NumMethod(NumMethod): t = time + self._adv_time_step - # Restore to a fixed depth or height? - if _depth_restoring is True: - zmax = data_reader.get_zmax(t, &_particle_copy_b) - _particle_copy_b.set_x3(_fixed_depth_below_surface + zmax) - - if data_reader.is_wet(t, &_particle_copy_b) == 1: - - # Determine the new host zlayer - flag = data_reader.set_vertical_grid_vars(t, &_particle_copy_b) - - # Apply surface/bottom boundary conditions if required - if flag != IN_DOMAIN: - flag = self._vert_bc_calculator.apply(data_reader, t, &_particle_copy_b) + # Has the particle beached? + if data_reader.is_wet(t, &_particle_copy_b) == 0: + if self._allow_beaching: + _particle_copy_b.set_is_beached(1) - # Return if failure recorded - if flag != IN_DOMAIN: - return flag + particle[0] = _particle_copy_b + return flag else: - _particle_copy_b.set_is_beached(1) - - # Copy back particle properties - particle[0] = _particle_copy_b + data_reader.find_nearest_wet_host(t, &_particle_copy_b) - return flag + # Restore to a fixed depth or height? + if _depth_restoring is True: + zmax = data_reader.get_zmax(t, &_particle_copy_b) + _particle_copy_b.set_x3(_fixed_depth_below_surface + zmax) elif _height_restoring is True: zmin = data_reader.get_zmin(t, &_particle_copy_b) _particle_copy_b.set_x3(_fixed_height_above_bed + zmin) - if data_reader.is_wet(t, &_particle_copy_b) == 1: - - # Determine the new host zlayer - flag = data_reader.set_vertical_grid_vars(t, &_particle_copy_b) - - # Apply surface/bottom boundary conditions if required - if flag != IN_DOMAIN: - flag = self._vert_bc_calculator.apply(data_reader, t, &_particle_copy_b) - - # Return if failure recorded - if flag != IN_DOMAIN: - return flag - - else: - _particle_copy_b.set_is_beached(1) - - # Copy back particle properties - particle[0] = _particle_copy_b - - return flag - - # Set vertical grid vars if the particle is in water - if data_reader.is_wet(t, &_particle_copy_b) == 1: + # Determine the new host zlayer + flag = data_reader.set_vertical_grid_vars(t, &_particle_copy_b) - flag = data_reader.set_vertical_grid_vars(t, &_particle_copy_b) + if flag != IN_DOMAIN: + # TODO - Depth/height restoring may fail in shallow water. Apply BC or return an error? + flag = self._vert_bc_calculator.apply(data_reader, t, &_particle_copy_b) - # Apply surface/bottom boundary conditions if required + # Return if failure recorded if flag != IN_DOMAIN: - flag = self._vert_bc_calculator.apply(data_reader, t, &_particle_copy_b) - - # Return if failure recorded - if flag != IN_DOMAIN: - return flag - else: - _particle_copy_b.set_is_beached(1) + return flag # Copy back particle properties particle[0] = _particle_copy_b