kernels.py

from parcels import StatusCode
import math

# Copied kernels from plasticparcels
def PolyTEOS10_bsq(particle, fieldset, time):
    """A seawater density kernel.

    Description:
    ----------
    A kernel to calculate the seawater density surrounding a particle based on
    the polyTEOS10-bsq algorithm from Appendix A.1 and A.2 of
    https://doi.org/10.1016/j.ocemod.2015.04.002

    Parameter Requirements
    ----------
    particle :
        - seawater_density
    fieldset :
        - `fieldset.conservative_temperature` - The conservative temperature field. Units [C].
        - `fieldset.absolute_salinity` - The absolute salinity field. Units [g/kg].

    Kernel Requirements:
    ----------
    Order of Operations:
        This kernel must be run before any kernel that requires an updated `particle.seawater_density` value.

    References
    ----------
     Roquet, F., Madec, G., McDougall, T. J., Barker, P. M., 2014: Accurate
      polynomial expressions for the density and specific volume of
      seawater using the TEOS-10 standard. Ocean Modelling.
     McDougall, T. J., D. R. Jackett, D. G. Wright and R. Feistel, 2003:
      Accurate and computationally efficient algorithms for potential
      temperature and density of seawater.  Journal of Atmospheric and
      Oceanic Technology, 20, 730-741.

    """

    Z = - math.fabs(particle.depth)  # note: use negative depths!
    SA = fieldset.absolute_salinity[time, particle.depth, particle.lat, particle.lon]
    CT = fieldset.conservative_temperature[time, particle.depth, particle.lat, particle.lon]

    SAu = 40. * 35.16504 / 35.
    CTu = 40.
    Zu = 1.0e+04
    deltaS = 32.

    zz = -Z / Zu
    R00 = 4.6494977072e+01
    R01 = -5.2099962525e+00
    R02 = 2.2601900708e-01
    R03 = 6.4326772569e-02
    R04 = 1.5616995503e-02
    R05 = -1.7243708991e-03
    r0 = (((((R05 * zz + R04) * zz + R03) * zz + R02) * zz + R01) * zz + R00) * zz

    R000 = 8.0189615746e+02
    R100 = 8.6672408165e+02
    R200 = -1.7864682637e+03
    R300 = 2.0375295546e+03
    R400 = -1.2849161071e+03
    R500 = 4.3227585684e+02
    R600 = -6.0579916612e+01
    R010 = 2.6010145068e+01
    R110 = -6.5281885265e+01
    R210 = 8.1770425108e+01
    R310 = -5.6888046321e+01
    R410 = 1.7681814114e+01
    R510 = -1.9193502195e+00
    R020 = -3.7074170417e+01
    R120 = 6.1548258127e+01
    R220 = -6.0362551501e+01
    R320 = 2.9130021253e+01
    R420 = -5.4723692739e+00
    R030 = 2.1661789529e+01
    R130 = -3.3449108469e+01
    R230 = 1.9717078466e+01
    R330 = -3.1742946532e+00
    R040 = -8.3627885467e+00
    R140 = 1.1311538584e+01
    R240 = -5.3563304045e+00
    R050 = 5.4048723791e-01
    R150 = 4.8169980163e-01
    R060 = -1.9083568888e-01
    R001 = 1.9681925209e+01
    R101 = -4.2549998214e+01
    R201 = 5.0774768218e+01
    R301 = -3.0938076334e+01
    R401 = 6.6051753097e+00
    R011 = -1.3336301113e+01
    R111 = -4.4870114575e+00
    R211 = 5.0042598061e+00
    R311 = -6.5399043664e-01
    R021 = 6.7080479603e+00
    R121 = 3.5063081279e+00
    R221 = -1.8795372996e+00
    R031 = -2.4649669534e+00
    R131 = -5.5077101279e-01
    R041 = 5.5927935970e-01
    R002 = 2.0660924175e+00
    R102 = -4.9527603989e+00
    R202 = 2.5019633244e+00
    R012 = 2.0564311499e+00
    R112 = -2.1311365518e-01
    R022 = -1.2419983026e+00
    R003 = -2.3342758797e-02
    R103 = -1.8507636718e-02
    R013 = 3.7969820455e-01
    ss = math.sqrt((SA + deltaS) / SAu)
    tt = CT / CTu
    zz = -Z / Zu
    rz3 = R013 * tt + R103 * ss + R003
    rz2 = (R022 * tt + R112 * ss + R012) * tt + (R202 * ss + R102) * ss + R002
    rz1 = (((R041 * tt + R131 * ss + R031) * tt + (R221 * ss + R121) * ss + R021) * tt + ((R311 * ss + R211) * ss + R111) * ss + R011) * tt + (((R401 * ss + R301) * ss + R201) * ss + R101) * ss + R001
    rz0 = (((((R060 * tt + R150 * ss + R050) * tt + (R240 * ss + R140) * ss + R040) * tt + ((R330 * ss + R230) * ss + R130) * ss + R030) * tt + (((R420 * ss + R320) * ss + R220) * ss + R120) * ss + R020) * tt + ((((R510 * ss + R410) * ss + R310) * ss + R210) * ss + R110) * ss + R010) * tt + (((((R600 * ss + R500) * ss + R400) * ss + R300) * ss + R200) * ss + R100) * ss + R000
    r = ((rz3 * zz + rz2) * zz + rz1) * zz + rz0

    particle.seawater_density = r0 + r



def PolyTEOS10_bsq_stif(particle, fieldset, time):
    """A seawater density kernel.

    Description:
    ----------
    A kernel to calculate the seawater density surrounding a particle based on
    the polyTEOS10-bsq algorithm from Appendix B.2 of
    https://doi.org/10.1016/j.ocemod.2015.04.002

    Parameter Requirements
    ----------
    particle :
        - seawater_density
    fieldset :
        - `fieldset.conservative_temperature` - The conservative temperature field. Units [C].
        - `fieldset.absolute_salinity` - The absolute salinity field. Units [g/kg].

    Kernel Requirements:
    ----------
    Order of Operations:
        This kernel must be run before any kernel that requires an updated `particle.seawater_density` value.

    References
    ----------
     Roquet, F., Madec, G., McDougall, T. J., Barker, P. M., 2014: Accurate
      polynomial expressions for the density and specific volume of
      seawater using the TEOS-10 standard. Ocean Modelling.
     McDougall, T. J., D. R. Jackett, D. G. Wright and R. Feistel, 2003:
      Accurate and computationally efficient algorithms for potential
      temperature and density of seawater.  Journal of Atmospheric and
      Oceanic Technology, 20, 730-741.

    """

    SAu = 40. * 35.16504 / 35.
    CTu = 40.
    Zu = 1.0e+04
    deltaS = 32.

    Z = - particle.depth  # note: use negative depths!
    SA = fieldset.absolute_salinity[time, particle.depth, particle.lat, particle.lon]
    CT = fieldset.conservative_temperature[time, particle.depth, particle.lat, particle.lon]

    zz = - Z / Zu
    R10 = 4.5238001132e-02
    R11 = -5.0691457704e-03
    R12 = 2.1990865986e-04
    R13 = 6.2587720090e-05
    R14 = 1.5194795322e-05
    R15 = -1.6777531159e-06
    r1 = (((((R15 * zz+R14) * zz+R13 ) * zz+R12 ) * zz+R11) * zz+R10) * zz + 1.0

    R000 = 8.0185969865e+02
    R100 = 8.6694400009e+02
    R200 = -1.7869886804e+03
    R300 = 2.0381548497e+03
    R400 = -1.2853207957e+03
    R500 = 4.3240996601e+02
    R600 = -6.0597694893e+01
    R010 = 2.6018938602e+01
    R110 = -6.5349778881e+01
    R210 = 8.1938301114e+01
    R310 = -5.7075044019e+01
    R410 = 1.7778972887e+01
    R510 = -1.9385277061e+00
    R020 = -3.7047586671e+01
    R120 = 6.1469677240e+01
    R220 = -6.0273564965e+01
    R320 = 2.9086148686e+01
    R420 = -5.4641152062e+00
    R030 = 2.1645370911e+01
    R130 = -3.3415216177e+01
    R230 = 1.9694119502e+01
    R330 = -3.1710487657e+00
    R040 = -8.3587252144e+00
    R140 = 1.1301873449e+01
    R240 = -5.3494910385e+00
    R050 = 5.4258439083e-01
    R150 = 4.7964142211e-01
    R060 = -1.9098973897e-01
    R001 = 2.1989266503e+01
    R101 = -4.2043785572e+01
    R201 = 4.8565183139e+01
    R301 = -3.0473875201e+01
    R401 = 6.5025798027e+00
    R011 = -1.3731591552e+01
    R111 = -4.3667272405e+00
    R211 = 5.2899251242e+00
    R311 = -7.1323430544e-01
    R021 = 7.4843328725e+00
    R121 = 3.1443007263e+00
    R221 = -1.8141790613e+00
    R031 = -2.6010169530e+00
    R131 = -4.9866824866e-01
    R041 = 5.5882379134e-01
    R002 = 1.1144133231e+00
    R102 = -4.5413503279e+00
    R202 = 2.7242118606e+00
    R012 = 2.8508453933e+00
    R112 = -4.4471541371e-01
    R022 = -1.5059281279e+00
    R003 = 1.9817129968e-01
    R103 = -1.7905418130e-01
    R013 = 2.5254195922e-01

    ss = math.sqrt((SA + deltaS) / SAu)
    tt = CT / CTu
    zz = -Z / Zu
    rz3 = R013 * tt + R103 * ss + R003
    rz2 = (R022 * tt + R112 * ss + R012) * tt + (R202 * ss + R102) * ss + R002
    rz1 = (((R041 * tt + R131 * ss + R031) * tt + (R221 * ss + R121) * ss + R021) * tt + ((R311 * ss + R211) * ss + R111) * ss + R011) * tt + (((R401 * ss + R301) * ss + R201) * ss + R101) * ss + R001
    rz0 = (((((R060 * tt + R150 * ss + R050) * tt + (R240 * ss + R140) * ss + R040) * tt + ((R330 * ss + R230) * ss + R130) * ss + R030) * tt + (((R420 * ss + R320) * ss + R220) * ss + R120) * ss + R020) * tt + ((((R510 * ss + R410) * ss + R310) * ss + R210) * ss + R110) * ss + R010) * tt + (((((R600 * ss + R500) * ss + R400) * ss + R300) * ss + R200) * ss + R100) * ss + R000
    rdot = ((rz3 * zz + rz2) * zz + rz1) * zz + rz0

    particle.seawater_density = r1 * rdot


def StokesDrift(particle, fieldset, time):
    """Stokes drift kernel.

    Description
    ----------
    Using the approach in [1] assuming a Phillips wave spectrum to determine
    the depth dependent Stokes drift. Specifically, the 'Stokes drift velocity'
    :math:`u_s` is computed as per Eq. (19) in [1].

    We treat the Stokes drift as a linear addition to the velocity field
        :math:`u(x,t) = u_c(x,t) + C_s * u_s(x,t)`
    where :math:`u_c` is the current velocity, :math:`u_s` is the Stokes drift velocity,
    and :math:`C_s` is the depth-varying decay factor.

    For further description, see https://plastic.oceanparcels.org/en/latest/physicskernels.html#stokes-drift

    Parameter Requirements
    ----------
    fieldset :
        - `fieldset.Stokes_U` and `fieldset.Stokes_V`, the Stokes drift velocity fields. Units [m s-1]
        - `fieldset.wave_Tp`, the peak wave period field (:math:`T_p`). Units [s].

    Kernel Requirements
    ----------
    Order of Operations:
        Must be run after the advection kernels

    References
    ----------
    [1] Breivik (2016) - https://doi.org/10.1016/j.ocemod.2016.01.005

    """
    # Sample the peak wave period
    T_p = fieldset.wave_Tp[time, particle.depth, particle.lat, particle.lon]

    # Compute the local bathymetry / water depth with a margin of error
    local_bathymetry = fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon]

    # Since AdvectionRK2 has already been run, let's piggy-back off the sampled U/V velocities and not resample
    # This is a bit of a trick, but avoids redundant field sampling
    #(U_ocean, V_ocean) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
    U_ocean = u1
    V_ocean = v1

    # Only compute displacements if the peak wave period is large enough and the particle is in the water
    if T_p > 1E-14 and particle.depth < 0.99*local_bathymetry and math.fabs(U_ocean**2 + V_ocean**2) > 1E-14:
        # Sample the U / V components of Stokes drift
        stokes_U = fieldset.Stokes_U[time, particle.depth, particle.lat, particle.lon]
        stokes_V = fieldset.Stokes_V[time, particle.depth, particle.lat, particle.lon]

        # Peak wave frequency
        omega_p = 2. * math.pi / T_p

        # Peak wave number
        k_p = (omega_p ** 2) / fieldset.G

        # Repeated inner term of Eq. (19) - note depth is negative in this formulation, but model depths are positive by convention
        kp_z_2 = 2. * k_p * particle.depth

        # Decay factor in Eq. (19) -- Where beta=1 for the Phillips spectrum
        decay = math.exp(-kp_z_2) - math.sqrt(math.pi * kp_z_2) * math.erfc(math.sqrt(kp_z_2))

        # Apply Eq. (19) and compute particle displacement
        particle_dlon += stokes_U * decay * particle.dt  # noqa
        particle_dlat += stokes_V * decay * particle.dt  # noqa

def StokesDriftCopernicus(particle, fieldset, time):
    """Stokes drift kernel.

    Description
    ----------
    Using the approach in [1] assuming a Phillips wave spectrum to determine
    the depth dependent Stokes drift. Specifically, the 'Stokes drift velocity'
    :math:`u_s` is computed as per Eq. (19) in [1].

    We treat the Stokes drift as a linear addition to the velocity field
        :math:`u(x,t) = u_c(x,t) + C_s * u_s(x,t)`
    where :math:`u_c` is the current velocity, :math:`u_s` is the Stokes drift velocity,
    and :math:`C_s` is the depth-varying decay factor.

    For further description, see https://plastic.oceanparcels.org/en/latest/physicskernels.html#stokes-drift

    Parameter Requirements
    ----------
    fieldset :
        - `fieldset.Stokes_U` and `fieldset.Stokes_V`, the Stokes drift velocity fields. Units [m s-1]
        - `fieldset.wave_Tp`, the peak wave period field (:math:`T_p`). Units [s].

    Kernel Requirements
    ----------
    Order of Operations:
        Must be run after the advection kernels

    References
    ----------
    [1] Breivik (2016) - https://doi.org/10.1016/j.ocemod.2016.01.005

    """
    # Sample the peak wave period
    T_p = fieldset.wave_Tp[time, particle.depth, particle.lat, particle.lon]

    # Compute the local bathymetry / water depth with a margin of error
    local_bathymetry = fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon]

    # Since AdvectionRK2 has already been run, let's piggy-back off the sampled U/V velocities and not resample
    # This is a bit of a trick, but avoids redundant field sampling
    #(U_ocean, V_ocean) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
    U_ocean = u1
    V_ocean = v1

    # Only compute displacements if the peak wave period is large enough and the particle is in the water
    if T_p > 1E-14 and particle.depth < 0.99*local_bathymetry and math.fabs(U_ocean**2 + V_ocean**2) > 1E-14 and particle.lat<89.80: # max lat of Copernicus wave data
        # Sample the U / V components of Stokes drift
        stokes_U = fieldset.Stokes_U[time, particle.depth, particle.lat, particle.lon]
        stokes_V = fieldset.Stokes_V[time, particle.depth, particle.lat, particle.lon]

        # Peak wave frequency
        omega_p = 2. * math.pi / T_p

        # Peak wave number
        k_p = (omega_p ** 2) / fieldset.G

        # Repeated inner term of Eq. (19) - note depth is negative in this formulation, but model depths are positive by convention
        kp_z_2 = 2. * k_p * particle.depth

        # Decay factor in Eq. (19) -- Where beta=1 for the Phillips spectrum
        decay = math.exp(-kp_z_2) - math.sqrt(math.pi * kp_z_2) * math.erfc(math.sqrt(kp_z_2))

        # Apply Eq. (19) and compute particle displacement
        particle_dlon += stokes_U * decay * particle.dt  # noqa
        particle_dlat += stokes_V * decay * particle.dt  # noqa

def Biofouling(particle, fieldset, time):
    r"""Settling velocity due to biofouling kernel.

    Description
    ----------
    Kernel to compute the settling velocity of particles due to changes in ambient algal
    concentrations, growth and death of attached algae based on [2]. The settling velocity
    of the particle is computed as per [1], however the particle size and density is
    affected by a biofouling process. The algae attached to the particle :math:`A` has a growth
    rate described by
        :math:`\frac{dA}{dt} = C + G - M - R`

    Here, :math:`C` models fouling of the plastic through collision with algae
        :math:`C = \beta_A \cdot A_A / \theta_{pl}`,
    where :math:`\beta_A` is the encounter rate, :math:`A_A` is the ambient algal concentration, and
    :math:`\theta_{pl}` is the surface area of plastic particle.

    :math:`G` models the growth of the algae attached to the surface of the particle
        :math:`G = \mu_A \cdot A`,
    where :math:`\mu_A` is the algal growth rate.

    :math:`M` models the grazing mortality of the algae attached to the surface of the particle
        :math:`M = m_A \cdot A`,
    where :math:`m_A` is the algal mortality rate.

    :math:`R` models the respiration of the algae attached to the surface of the particle
        :math:`R = Q_{10}^{(T-20)/10}R_{20}A`,
    where :math:`Q_{10}` is the temperature coefficient, which indicates how much the respiration
    increases when the temperature increases by 10C. :math:`T` is the temperature, and :math:`R_{20}`
    is the respiration rate.

    Calculation steps:
        1. Compute the seawater dynamic viscosity from Eq. (27) in [2]
        2. Compute the kinematic viscosity from Eq. (25) in [2]
        3. Compute the dimensionless particle diameter from Eq. (4) in [2], equivalent to Eq. (6) in [1]
        4. Compute the dimensionless settling velocity from Eq. (3) in [2], equivalent to Eq. (8) in [1]
        5. Compute the settling velocity of the particle from Eq. (2) in [2], equivalent to Eq. (9) in [1]


    Parameter Requirements
    ----------
    particle :
        - diameter
        - density
        - seawater_density
    fieldset :
        - `fieldset.G` - Gravity constant. Units [m s-2].
        - `fieldset.conservative_temperature` - The conservative temperature field. Units [C].
        - `fieldset.absolute_salinity` - The absolute salinity field. Units [g/kg].
        - `fieldset.algae_cell_volume` - The volume of 1 algal cell [m-3]
        - `fieldset.biofilm_density` - The density of the biofilm [kg m-3]


    Kernel Requirements
    ----------
    Order of Operations:
        This kernel must run after the PolyTEOS10_bsq kernel, which sets the particle.seawater_density variable, relied on by this.

    References
    ----------
    [1] Dietrich (1982) - https://doi.org/10.1029/WR018i006p01615
    [2] Kooi et al. (2017) - https://doi.org/10.1021/acs.est.6b04702
    [3] Menden-Deuer and Lessard (2000) - https://doi.org/10.4319/lo.2000.45.3.0569
    [4] Bernard and Remond (2012) - https://doi.org/10.1016/j.biortech.2012.07.022
    """
    
    # Piggy-back of advection kernel for u/v sampling
    ocean_U = u1
    ocean_V = v1

    # Only apply the biofouling kernel if the particle is in the ocean
    # Otherwise, skip this timestep (even if the biofilm should remineralise, there will be issues with the temp/salt values on land!)
    if math.sqrt(ocean_U**2 + ocean_V**2) > 1e-14:
        # seawater_density = particle.seawater_density  # [kg m-3]
        # Piggy-back off the PolyTEOS10_bsq kernel for the temperature and salinity sampling
        temperature = CT #fieldset.conservative_temperature[time, particle.depth, particle.lat, particle.lon]
        seawater_salinity = SA/1000. #fieldset.absolute_salinity[time, particle.depth, particle.lat, particle.lon] / 1000. # Convert from [g/kg] to [kg/kg]
        
        particle_radius = 0.5 * particle.plastic_diameter
        # particle_density = particle.plastic_density
        initial_settling_velocity = particle.settling_velocity  # settling velocity [m s-1]

        # Compute the seawater dynamic viscosity and kinematic viscosity
        water_dynamic_viscosity = 4.2844E-5 + (1. / ((0.156 * (temperature + 64.993) ** 2) - 91.296))  # Eq. (26) from [2]
        A = 1.541 + 1.998E-2 * temperature - 9.52E-5 * temperature ** 2  # Eq. (28) from [2]
        B = 7.974 - 7.561E-2 * temperature + 4.724E-4 * temperature ** 2  # Eq. (29) from [2]
        seawater_dynamic_viscosity = water_dynamic_viscosity * (1. + A * seawater_salinity + B * seawater_salinity ** 2)  # Eq. (27) from [2]
        seawater_kinematic_viscosity = seawater_dynamic_viscosity / particle.seawater_density  # Eq. (25) from [2]

        # Compute the algal growth component
        # Sample fields
        mol_concentration_diatoms = fieldset.bio_diatom[time, particle.depth, particle.lat, particle.lon]  # Mole concentration of Diatoms expressed as carbon in sea water
        mol_concentration_nanophytoplankton = fieldset.bio_nanophy[time, particle.depth, particle.lat, particle.lon]  # Mole concentration of Nanophytoplankton expressed as carbon in seawater
        total_primary_production_of_phyto = fieldset.pp_phyto[time, particle.depth, particle.lat, particle.lon]  # mg C /m3/day #pp_phyto_
        median_mg_carbon_per_cell = 2726.0e-9  # Median mg of Carbon per cell selected from [3]. From [2] pg7966 - "The conversion from carbon to algae cells is highly variable, ranging between 35339 to 47.76 pg per carbon cell [3]. We choose the median value, 2726 x 10^-9 mg per carbon cell."
        # carbon_molecular_weight = fieldset.carbon_molecular_weight # grams C per mol of C #Wt_C

        # Compute concentration numbers
        number_concentration_diatoms = mol_concentration_diatoms * (fieldset.carbon_molecular_weight / median_mg_carbon_per_cell)  # conversion from [mmol C m-3] to [mg C m-3] to [no. m-3]
        number_concentration_diatoms = max(number_concentration_diatoms, 0.)  # Ensure non-negative
        number_concentration_nanophytoplankton = mol_concentration_nanophytoplankton * (fieldset.carbon_molecular_weight / median_mg_carbon_per_cell)
        number_concentration_nanophytoplankton = max(number_concentration_nanophytoplankton, 0.)  # Ensure non-negative
        number_concentration_total = number_concentration_diatoms + number_concentration_nanophytoplankton  # conversion from [mmol C m-3] to [mg C m-3] to [no. m-3]

        # Compute primary production
        primary_production_per_cell = total_primary_production_of_phyto / number_concentration_total  # primary productivity per cell, in mg C / cell / day
        primary_production_numcell_per_cell = primary_production_per_cell / median_mg_carbon_per_cell  # primary productivity in terms of amount of cells per cell, in cells / cell / day
        primary_production_numcell_per_cell = max(primary_production_numcell_per_cell, 0.)  # Ensure non-negative

        # Compute growth rates
        max_growth_rate = 1.85  # Maximum growth rate (per day), defined in Table S1 of [2], based on the nominal value estimated in Table 4 of [4]
        mu_a = min(primary_production_numcell_per_cell, max_growth_rate)/86400.  # Set growth rate per second

        # Compute the algal growth of the algae already on the particle
        algae_growth = mu_a * particle.algae_amount  # productivity in amount of cells/m2/s

        # Compute the radius, surface area, volume and thickness of the particle including potential biofilm
        particle_volume = (4. / 3.) * math.pi * particle_radius ** 3.  # volume of plastic [m3]
        particle_surface_area = 4. * math.pi * particle_radius ** 2.  # surface area of plastic particle [m2]
        algal_cell_radius = ((3. / 4.) * (fieldset.algae_cell_volume / math.pi)) ** (1. / 3.)  # radius of an algal cell [m]
        biofilm_volume = fieldset.algae_cell_volume * particle.algae_amount * particle_surface_area  # volume of living biofilm [m3]

        total_volume = biofilm_volume + particle_volume  # volume of total (biofilm + plastic) [m3]
        total_radius = ((total_volume * (3. / (4. * math.pi))) ** (1. / 3.))  # total radius [m]

        # Compute diffusivities
        total_density = (particle_volume * particle.plastic_density + biofilm_volume * fieldset.biofilm_density) / total_volume  # total density [kg m-3]
        plastic_diffusivity = fieldset.K * (temperature + 273.16) / (6. * math.pi * seawater_dynamic_viscosity * total_radius)  # diffusivity of plastic particle [m2 s-1]
        algae_diffusivity = fieldset.K * (temperature + 273.16) / (6. * math.pi * seawater_dynamic_viscosity * algal_cell_radius)  # diffusivity of algal cells [m2 s-1]

        # Compute the encounter rates
        beta_abrown = 4. * math.pi * (plastic_diffusivity + algae_diffusivity) * (total_radius + algal_cell_radius)  # Brownian motion [m3 s-1]
        beta_ashear = 1.3 * fieldset.Gamma * ((total_radius + algal_cell_radius) ** 3.)  # advective shear [m3 s-1]
        beta_aset = (1. / 2.) * math.pi * total_radius ** 2. * math.fabs(initial_settling_velocity)  # differential settling [m3 s-1]
        beta_a = beta_abrown + beta_ashear + beta_aset  # collision rate [m3 s-1]

        # Compute the algal growth via collision
        a_collision = fieldset.collision_probability * (beta_a * number_concentration_diatoms) / particle_surface_area  # [no. m-2 s-1] collisions with diatoms (Eq. 11 in [2])

        # Compute the algal decay due to respiration
        a_respiration = fieldset.algae_respiration_f * (fieldset.Q10 ** ((temperature - 20.) / 10.)) * fieldset.R20 * particle.algae_amount  # [no. m-2 s-1] respiration

        # Compute the algal decay due to grazing
        a_grazing = fieldset.algae_mortality_rate * particle.algae_amount

        # Compute the final algal amount
        algae_amount_change = (a_collision + algae_growth - a_grazing - a_respiration) * particle.dt
        if particle.algae_amount + algae_amount_change < 0.:
            particle.algae_amount = 0.
        else:
            particle.algae_amount += algae_amount_change

        # Compute the new settling velocity
        particle_diameter = 2. * (total_radius)  # equivalent spherical diameter [m], calculated from Dietrich (1982) from A = pi/4 * dn**2

        # Compute the density difference of the particle
        normalised_density_difference = (total_density - particle.seawater_density) / particle.seawater_density  # normalised difference in density between the plastic particle and seawater [-]

        # Compute the dimensionless particle diameter D_* using Eq. (4) from [2]
        dimensionless_diameter = (math.fabs(total_density - particle.seawater_density) * fieldset.G * particle_diameter ** 3.) / (particle.seawater_density * seawater_kinematic_viscosity ** 2.)  # [-]

        # Compute the dimensionless settling velocity w_*
        if dimensionless_diameter > 5E9:  # "The boundary layer around the sphere becomes fully turbulent, causing a reduction in drag and an increase in settling velocity" - [1]
            dimensionless_velocity = 265000.  # Set a maximum dimensionless settling velocity
        elif dimensionless_diameter < 0.05:  # "At values of D_* less than 0.05, (9) deviates signficantly ... from Stokes' law and (8) should be used." - [1]
            dimensionless_velocity = (dimensionless_diameter ** 2.) / 5832.  # Using Eq. (8) in [1]
        else:
            dimensionless_velocity = 10. ** (-3.76715 + (1.92944 * math.log10(dimensionless_diameter)) - (0.09815 * math.log10(dimensionless_diameter) ** 2.) - (
                                            0.00575 * math.log10(dimensionless_diameter) ** 3.) + (0.00056 * math.log10(dimensionless_diameter) ** 4.))  # Using Eq. (9) in [1]

        # Compute the settling velocity of the particle using Eq. (5) from [1] (solving for the settling velocity)
        sign_of_density_difference = math.copysign(1., normalised_density_difference)
        settling_velocity = sign_of_density_difference * (fieldset.G * seawater_kinematic_viscosity * dimensionless_velocity * math.fabs(normalised_density_difference)) ** (1. / 3.)  # m s-1

        # Update the settling velocity
        particle.settling_velocity = settling_velocity

        # Update particle depth
        particle_ddepth += particle.settling_velocity * particle.dt  # noqa


def BiofoulingCopernicus(particle, fieldset, time):
    r"""Settling velocity due to biofouling kernel using Copernicus Hindcast data (no diatoms)

    Description
    ----------
    Kernel to compute the settling velocity of particles due to changes in ambient algal
    concentrations, growth and death of attached algae based on [2]. The settling velocity
    of the particle is computed as per [1], however the particle size and density is
    affected by a biofouling process. The algae attached to the particle :math:`A` has a growth
    rate described by
        :math:`\frac{dA}{dt} = C + G - M - R`

    Here, :math:`C` models fouling of the plastic through collision with algae
        :math:`C = \beta_A \cdot A_A / \theta_{pl}`,
    where :math:`\beta_A` is the encounter rate, :math:`A_A` is the ambient algal concentration, and
    :math:`\theta_{pl}` is the surface area of plastic particle.

    :math:`G` models the growth of the algae attached to the surface of the particle
        :math:`G = \mu_A \cdot A`,
    where :math:`\mu_A` is the algal growth rate.

    :math:`M` models the grazing mortality of the algae attached to the surface of the particle
        :math:`M = m_A \cdot A`,
    where :math:`m_A` is the algal mortality rate.

    :math:`R` models the respiration of the algae attached to the surface of the particle
        :math:`R = Q_{10}^{(T-20)/10}R_{20}A`,
    where :math:`Q_{10}` is the temperature coefficient, which indicates how much the respiration
    increases when the temperature increases by 10C. :math:`T` is the temperature, and :math:`R_{20}`
    is the respiration rate.

    Calculation steps:
        1. Compute the seawater dynamic viscosity from Eq. (27) in [2]
        2. Compute the kinematic viscosity from Eq. (25) in [2]
        3. Compute the dimensionless particle diameter from Eq. (4) in [2], equivalent to Eq. (6) in [1]
        4. Compute the dimensionless settling velocity from Eq. (3) in [2], equivalent to Eq. (8) in [1]
        5. Compute the settling velocity of the particle from Eq. (2) in [2], equivalent to Eq. (9) in [1]


    Parameter Requirements
    ----------
    particle :
        - diameter
        - density
        - seawater_density
    fieldset :
        - `fieldset.G` - Gravity constant. Units [m s-2].
        - `fieldset.conservative_temperature` - The conservative temperature field. Units [C].
        - `fieldset.absolute_salinity` - The absolute salinity field. Units [g/kg].
        - `fieldset.algae_cell_volume` - The volume of 1 algal cell [m-3]
        - `fieldset.biofilm_density` - The density of the biofilm [kg m-3]


    Kernel Requirements
    ----------
    Order of Operations:
        This kernel must run after the PolyTEOS10_bsq kernel, which sets the particle.seawater_density variable, relied on by this.

    References
    ----------
    [1] Dietrich (1982) - https://doi.org/10.1029/WR018i006p01615
    [2] Kooi et al. (2017) - https://doi.org/10.1021/acs.est.6b04702
    [3] Menden-Deuer and Lessard (2000) - https://doi.org/10.4319/lo.2000.45.3.0569
    [4] Bernard and Remond (2012) - https://doi.org/10.1016/j.biortech.2012.07.022
    """
    
    # Piggy-back of advection kernel for u/v sampling
    ocean_U = u1
    ocean_V = v1

    # Only apply the biofouling kernel if the particle is in the ocean
    # Otherwise, skip this timestep (even if the biofilm should remineralise, there will be issues with the temp/salt values on land!)
    if math.sqrt(ocean_U**2 + ocean_V**2) > 1e-14:
        # seawater_density = particle.seawater_density  # [kg m-3]
        # Piggy-back off the PolyTEOS10_bsq kernel for the temperature and salinity sampling
        temperature = CT #fieldset.conservative_temperature[time, particle.depth, particle.lat, particle.lon]
        seawater_salinity = SA/1000. #fieldset.absolute_salinity[time, particle.depth, particle.lat, particle.lon] / 1000. # Convert from [g/kg] to [kg/kg]
        
        particle_radius = 0.5 * particle.plastic_diameter
        # particle_density = particle.plastic_density
        initial_settling_velocity = particle.settling_velocity  # settling velocity [m s-1]

        # Compute the seawater dynamic viscosity and kinematic viscosity
        water_dynamic_viscosity = 4.2844E-5 + (1. / ((0.156 * (temperature + 64.993) ** 2) - 91.296))  # Eq. (26) from [2]
        A = 1.541 + 1.998E-2 * temperature - 9.52E-5 * temperature ** 2  # Eq. (28) from [2]
        B = 7.974 - 7.561E-2 * temperature + 4.724E-4 * temperature ** 2  # Eq. (29) from [2]
        seawater_dynamic_viscosity = water_dynamic_viscosity * (1. + A * seawater_salinity + B * seawater_salinity ** 2)  # Eq. (27) from [2]
        seawater_kinematic_viscosity = seawater_dynamic_viscosity / particle.seawater_density  # Eq. (25) from [2]

        # Compute the algal growth component
        # Sample fields
        mol_concentration_diatoms = 0. #There is no diatom field in Copernicus hindcast data. #fieldset.bio_diatom[time, particle.depth, particle.lat, particle.lon]  # Mole concentration of Diatoms expressed as carbon in sea water
        mol_concentration_nanophytoplankton = fieldset.bio_nanophy[time, particle.depth, particle.lat, particle.lon]  # Mole concentration of Nanophytoplankton expressed as carbon in seawater
        total_primary_production_of_phyto = fieldset.pp_phyto[time, particle.depth, particle.lat, particle.lon]  # mg C /m3/day #pp_phyto_
        median_mg_carbon_per_cell = 2726.0e-9  # Median mg of Carbon per cell selected from [3]. From [2] pg7966 - "The conversion from carbon to algae cells is highly variable, ranging between 35339 to 47.76 pg per carbon cell [3]. We choose the median value, 2726 x 10^-9 mg per carbon cell."
        # carbon_molecular_weight = fieldset.carbon_molecular_weight # grams C per mol of C #Wt_C

        # Compute concentration numbers
        number_concentration_diatoms = mol_concentration_diatoms * (fieldset.carbon_molecular_weight / median_mg_carbon_per_cell)  # conversion from [mmol C m-3] to [mg C m-3] to [no. m-3]
        number_concentration_diatoms = max(number_concentration_diatoms, 0.)  # Ensure non-negative
        number_concentration_nanophytoplankton = mol_concentration_nanophytoplankton * (fieldset.carbon_molecular_weight / median_mg_carbon_per_cell)
        number_concentration_nanophytoplankton = max(number_concentration_nanophytoplankton, 0.)  # Ensure non-negative
        number_concentration_total = number_concentration_diatoms + number_concentration_nanophytoplankton  # conversion from [mmol C m-3] to [mg C m-3] to [no. m-3]

        # Compute primary production
        primary_production_per_cell = total_primary_production_of_phyto / number_concentration_total  # primary productivity per cell, in mg C / cell / day
        primary_production_numcell_per_cell = primary_production_per_cell / median_mg_carbon_per_cell  # primary productivity in terms of amount of cells per cell, in cells / cell / day
        primary_production_numcell_per_cell = max(primary_production_numcell_per_cell, 0.)  # Ensure non-negative

        # Compute growth rates
        max_growth_rate = 1.85  # Maximum growth rate (per day), defined in Table S1 of [2], based on the nominal value estimated in Table 4 of [4]
        mu_a = min(primary_production_numcell_per_cell, max_growth_rate)/86400.  # Set growth rate per second

        # Compute the algal growth of the algae already on the particle
        algae_growth = mu_a * particle.algae_amount  # productivity in amount of cells/m2/s

        # Compute the radius, surface area, volume and thickness of the particle including potential biofilm
        particle_volume = (4. / 3.) * math.pi * particle_radius ** 3.  # volume of plastic [m3]
        particle_surface_area = 4. * math.pi * particle_radius ** 2.  # surface area of plastic particle [m2]
        algal_cell_radius = ((3. / 4.) * (fieldset.algae_cell_volume / math.pi)) ** (1. / 3.)  # radius of an algal cell [m]
        biofilm_volume = fieldset.algae_cell_volume * particle.algae_amount * particle_surface_area  # volume of living biofilm [m3]

        total_volume = biofilm_volume + particle_volume  # volume of total (biofilm + plastic) [m3]
        total_radius = ((total_volume * (3. / (4. * math.pi))) ** (1. / 3.))  # total radius [m]

        # Compute diffusivities
        total_density = (particle_volume * particle.plastic_density + biofilm_volume * fieldset.biofilm_density) / total_volume  # total density [kg m-3]
        plastic_diffusivity = fieldset.K * (temperature + 273.16) / (6. * math.pi * seawater_dynamic_viscosity * total_radius)  # diffusivity of plastic particle [m2 s-1]
        algae_diffusivity = fieldset.K * (temperature + 273.16) / (6. * math.pi * seawater_dynamic_viscosity * algal_cell_radius)  # diffusivity of algal cells [m2 s-1]

        # Compute the encounter rates
        beta_abrown = 4. * math.pi * (plastic_diffusivity + algae_diffusivity) * (total_radius + algal_cell_radius)  # Brownian motion [m3 s-1]
        beta_ashear = 1.3 * fieldset.Gamma * ((total_radius + algal_cell_radius) ** 3.)  # advective shear [m3 s-1]
        beta_aset = (1. / 2.) * math.pi * total_radius ** 2. * math.fabs(initial_settling_velocity)  # differential settling [m3 s-1]
        beta_a = beta_abrown + beta_ashear + beta_aset  # collision rate [m3 s-1]

        # Compute the algal growth via collision
        a_collision = fieldset.collision_probability * (beta_a * number_concentration_diatoms) / particle_surface_area  # [no. m-2 s-1] collisions with diatoms (Eq. 11 in [2])

        # Compute the algal decay due to respiration
        a_respiration = fieldset.algae_respiration_f * (fieldset.Q10 ** ((temperature - 20.) / 10.)) * fieldset.R20 * particle.algae_amount  # [no. m-2 s-1] respiration

        # Compute the algal decay due to grazing
        a_grazing = fieldset.algae_mortality_rate * particle.algae_amount

        # Compute the final algal amount
        algae_amount_change = (a_collision + algae_growth - a_grazing - a_respiration) * particle.dt
        if particle.algae_amount + algae_amount_change < 0.:
            particle.algae_amount = 0.
        else:
            particle.algae_amount += algae_amount_change

        # Compute the new settling velocity
        particle_diameter = 2. * (total_radius)  # equivalent spherical diameter [m], calculated from Dietrich (1982) from A = pi/4 * dn**2

        # Compute the density difference of the particle
        normalised_density_difference = (total_density - particle.seawater_density) / particle.seawater_density  # normalised difference in density between the plastic particle and seawater [-]

        # Compute the dimensionless particle diameter D_* using Eq. (4) from [2]
        dimensionless_diameter = (math.fabs(total_density - particle.seawater_density) * fieldset.G * particle_diameter ** 3.) / (particle.seawater_density * seawater_kinematic_viscosity ** 2.)  # [-]

        # Compute the dimensionless settling velocity w_*
        if dimensionless_diameter > 5E9:  # "The boundary layer around the sphere becomes fully turbulent, causing a reduction in drag and an increase in settling velocity" - [1]
            dimensionless_velocity = 265000.  # Set a maximum dimensionless settling velocity
        elif dimensionless_diameter < 0.05:  # "At values of D_* less than 0.05, (9) deviates signficantly ... from Stokes' law and (8) should be used." - [1]
            dimensionless_velocity = (dimensionless_diameter ** 2.) / 5832.  # Using Eq. (8) in [1]
        else:
            dimensionless_velocity = 10. ** (-3.76715 + (1.92944 * math.log10(dimensionless_diameter)) - (0.09815 * math.log10(dimensionless_diameter) ** 2.) - (
                                            0.00575 * math.log10(dimensionless_diameter) ** 3.) + (0.00056 * math.log10(dimensionless_diameter) ** 4.))  # Using Eq. (9) in [1]

        # Compute the settling velocity of the particle using Eq. (5) from [1] (solving for the settling velocity)
        sign_of_density_difference = math.copysign(1., normalised_density_difference)
        settling_velocity = sign_of_density_difference * (fieldset.G * seawater_kinematic_viscosity * dimensionless_velocity * math.fabs(normalised_density_difference)) ** (1. / 3.)  # m s-1

        # Update the settling velocity
        particle.settling_velocity = settling_velocity

        # Update particle depth
        particle_ddepth += particle.settling_velocity * particle.dt  # noqa


# Because every displacement kernel is an Euler-forward scheme, except advection, the only kernel to throw OOB in the current timestep is the advection kernel.
def checkErrorThroughSurface_2DAdvectionRK4(particle, fieldset, time):
    # This is a kernel to handle 3D advection leading to ErrorThroughSurface
    if particle.state == StatusCode.ErrorThroughSurface:
        # Perform 2D horizontal advection only!
        """Advection of particles using fourth-order Runge-Kutta integration."""
        (u1, v1) = fieldset.UV[particle]
        lon1, lat1 = (particle.lon + u1 * 0.5 * particle.dt, particle.lat + v1 * 0.5 * particle.dt)
        (u2, v2) = fieldset.UV[time + 0.5 * particle.dt, particle.depth, lat1, lon1, particle]
        lon2, lat2 = (particle.lon + u2 * 0.5 * particle.dt, particle.lat + v2 * 0.5 * particle.dt)
        (u3, v3) = fieldset.UV[time + 0.5 * particle.dt, particle.depth, lat2, lon2, particle]
        lon3, lat3 = (particle.lon + u3 * particle.dt, particle.lat + v3 * particle.dt)
        (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
        particle_dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * particle.dt  # noqa
        particle_dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particle.dt  # noqa

        particle.state = StatusCode.Success



def AdvectionRK2_3D(particle, fieldset, time):  # pragma: no cover
    """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
    (u1, v1, w1) = fieldset.UVW[particle]
    lon1 = particle.lon + u1 * 0.5 * particle.dt
    lat1 = particle.lat + v1 * 0.5 * particle.dt
    dep1 = particle.depth + w1 * 0.5 * particle.dt
    (u2, v2, w2) = fieldset.UVW[time + 0.5 * particle.dt, dep1, lat1, lon1, particle]
    particle_dlon += u2 * particle.dt  # noqa
    particle_dlat += v2 * particle.dt  # noqa
    particle_ddepth += w2 * particle.dt  # noqa


def AdvectionRK2(particle, fieldset, time):  # pragma: no cover
    """Advection of particles using fourth-order Runge-Kutta integration."""
    (u1, v1) = fieldset.UV[particle]
    lon1 = particle.lon + u1 * 0.5 * particle.dt
    lat1 = particle.lat + v1 * 0.5 * particle.dt
    (u2, v2) = fieldset.UV[time + 0.5 * particle.dt, particle.depth, lat1, lon1, particle]
    particle_dlon += u2 * particle.dt  # noqa
    particle_dlat += v2 * particle.dt  # noqa


# Because every displacement kernel is an Euler-forward scheme, except advection, the only kernel to throw OOB in the current timestep is the advection kernel.
def checkErrorThroughSurface_2DAdvectionRK2(particle, fieldset, time):
    # This is a kernel to handle 3D advection leading to ErrorThroughSurface
    if particle.state == StatusCode.ErrorThroughSurface:
        # Perform 2D horizontal advection only!
        """Advection of particles using second-order Runge-Kutta integration."""
        (u1, v1) = fieldset.UV[particle]
        lon1, lat1 = (particle.lon + u1 * 0.5 * particle.dt, particle.lat + v1 * 0.5 * particle.dt)
        (u2, v2) = fieldset.UV[time + 0.5 * particle.dt, particle.depth, lat1, lon1, particle]
        particle_dlon += u2 * particle.dt  # noqa
        particle_dlat += v2 * particle.dt  # noqa

        particle.state = StatusCode.Success



def checkThroughSurface(particle, fieldset, time):
    """ Kernel to check if the new particle position will be through the surface. """
    updated_depth = particle.depth + particle_ddepth
    if updated_depth <= 0.:
        particle_ddepth = 0. # This will set the particle to the surface
        particle.depth = 0.

def checkThroughSurfaceCopernicus(particle, fieldset, time):
    """ Kernel to check if the new particle position will be through the surface. """
    updated_depth = particle.depth + particle_ddepth
    if updated_depth <= fieldset.z_start:
        particle_ddepth = 0. # This will set the particle to the surface
        particle.depth = fieldset.z_start

def deleteParticle(particle, fieldset, time):
    """ Kernel to delete a particle. """
    if particle.state >= 50:
        particle.delete()


def stopExecution(particle, fieldset, time):
    """ Kernel to stop execution for a particle. """
    if particle.state >= 50:
        particle.state = StatusCode.StopExecution

def belowLatitude(particle, fieldset, time):
    """ Kernel to stop execution for a particle below a certain latitude. """
    if particle.lat < 26. and particle.state == StatusCode.Success:
        particle.state = StatusCode.StopExecution

def periodicBC(particle, fieldset, time):
    if particle.lon < 0.:
        particle.lon = particle.lon + 360.
    elif particle.lon > 360.:
        particle.lon = particle.lon - 360.

# def reflectAtSurface(particle, fieldset, time):
#     """ Kernel to reflect the particle at the ocean surface if it goes through the surface.
#     """
#     potential_depth = particle.depth + particle_ddepth
#     if potential_depth < 0.:# Particle is above the surface
#         particle_ddepth = - potential_depth - particle.depth
#         # End of kernel loop, particle.depth = particle.depth - potential_depth - particle.depth = -potential_depth, so the reflection has been applied

def reflectAtSurface(particle, fieldset, time):
    """ Kernel to reflect the particle at the ocean surface if it goes through the surface.
    For simplicity, we will just set the depth, and ignore the particle_ddepth variable - this will lead to a warning
    """
    potential_depth = particle.depth + particle_ddepth
    if potential_depth < 0.:# Particle is above the surface
        particle.depth = -potential_depth
        particle_ddepth = 0. # Set particle_ddepth to 0, as we have already updated the depth


def reflectAtSurfaceCopernicus(particle, fieldset, time):
    """ Kernel to reflect the particle at the ocean surface if it goes through the surface.
    For simplicity, we will just set the depth, and ignore the particle_ddepth variable - this will lead to a warning
    """
    potential_depth = particle.depth + particle_ddepth
    if potential_depth < fieldset.z_start: # Particle is above the surface
        particle.depth = fieldset.z_start - potential_depth
        particle_ddepth = 0. # Set particle_ddepth to 0, as we have already updated the depth



def reflectAtBathymetry(particle, fieldset, time):
    """ Kernel to reflect the particle at the ocean bathymetry if it goes through the bathymetry.
    For simplicity, we will just set the depth, and ignore the particle_ddepth variable - this will lead to a warning
    """
    #local_bathymetry = fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon]
    # local_bathymetry is already computed in the StokesDrift kernel, so we can piggy-back off that value
    potential_depth = particle.depth + particle_ddepth

    if potential_depth > 100:
        local_bathymetry = 0.99*local_bathymetry # Handle linear interpolation issues for deep particles
        # NOTE: This sets the local_bathymetry variable to something different to what was sampled, be careful using the variable in later kernels!
    
    if potential_depth > local_bathymetry: # If the potential particle position is below the bathymetry
        beyond_depth = potential_depth - local_bathymetry # The extent to which the particle goes beyond the bathymetry
        particle.depth = local_bathymetry - beyond_depth # Reflect the particle back above the bathymetry
        particle_ddepth = 0. # Set particle_ddepth to 0, as we have already updated the depth



def unbeachingBySamplingAfterwards(particle, fieldset, time):
    """Unbeaching particles by sampling the velocity field and moving in the direction of velocities!
    Note, a particle can spend at least one full timestep on land before being unbeached, as we only
    unbeach when the dlon and dlat are both (near) zero.
    """
    # Use velocities from the advection kernel
    vel_u = u1
    vel_v = v1

    # in the case of being beached, these displacement will be zero
    new_lon = particle.lon + particle_dlon
    new_lat = particle.lat + particle_dlat
    new_depth = particle.depth + particle_ddepth

    #if math.sqrt(vel_u**2 + vel_v**2) < 1e-14: 
    if math.fabs(particle_dlon) + math.fabs(particle_dlat) < 1e-14: # The particle hasn't moved more than horizontally
        # Case 1: we move onto land -> most likely case because particles should be initialised in water!
        #if math.fabs(particle_dlon) > 1e-14 or math.fabs(particle_dlat) > 1e-14:
        #    # Case 2: We have moved onto land, so rather, don't move, but still allow for the vertical movement.
        #    particle_dlon = 0.
        #    particle_dlat = 0.

        ##Case 2: In the rare case that we start on land -> We haven't moved and our velocity is zero - we need to find which direction to move in.
        #else: # if math.fabs(particle_dlon) == 0. and math.fabs(particle_dlat) == 0.:
        displacement = 1./8. # Degree displacement to sample the velocity field

        unbeach_U = 1. / (1852. * 60. * math.cos(particle.lat * math.pi / 180.)) # Convert 1m/2s to degrees/s at the particle latitude in zonal direction
        unbeach_V = 1. / (1852. * 60.) # Convert 1m/2s to degrees/s in meridional direction

        (U_left, V_left) = fieldset.UV[time, new_depth, new_lat, new_lon - displacement]
        (U_right, V_right) = fieldset.UV[time, new_depth, new_lat, new_lon + displacement]
        (U_up, V_up) = fieldset.UV[time, new_depth, new_lat + displacement, new_lon]
        (U_down, V_down) = fieldset.UV[time, new_depth, new_lat - displacement, new_lon]

        # Find the direction of the highest velocity
        left = math.sqrt(U_left**2 + V_left**2)
        right = math.sqrt(U_right**2 + V_right**2)
        up = math.sqrt(U_up**2 + V_up**2)
        down = math.sqrt(U_down**2 + V_down**2)
        
        max_vel = 0.
        U_dir = 0.
        V_dir = 0.
        
        if left + right + up + down > 1e-14:
            # Annoyingly in C, we can't do python like comparisons, so be basic.
            # Let's just iterate over the four possibilities
            max_vel = left
            U_dir = -1.
            V_dir = 0.
            if max_vel < right:
                max_vel = right
                U_dir = 1.
                V_dir = 0.
            if max_vel < up:
                max_vel = up
                U_dir = 0.
                V_dir = 1.
            if max_vel < down:
                max_vel = down
                U_dir = 0.
                V_dir = -1.

        else: # In the case that all four directions are zero - i.e., we are stuck in a corner
            (U_left_up, V_left_up) = fieldset.UV[time, new_depth, new_lat + displacement, new_lon - displacement]
            (U_left_down, V_left_down) = fieldset.UV[time, new_depth, new_lat - displacement, new_lon - displacement]
            (U_right_up, V_right_up) = fieldset.UV[time, new_depth, new_lat + displacement, new_lon + displacement]
            (U_right_down, V_right_down) = fieldset.UV[time, new_depth, new_lat - displacement, new_lon + displacement]

            left_up = math.sqrt(U_left_up**2 + V_left_up**2)
            left_down = math.sqrt(U_left_down**2 + V_left_down**2)
            right_up = math.sqrt(U_right_up**2 + V_right_up**2)
            right_down = math.sqrt(U_right_down**2 + V_right_down**2)

            max_vel = left_up
            U_dir = -1.
            V_dir = 1.
            if max_vel < left_down:
                max_vel = left_down
                U_dir = 1.
                V_dir = -1.
            if max_vel < right_up:
                max_vel = right_up
                U_dir = 1.
                V_dir = 1.
            if max_vel < right_down:
                max_vel = right_down
                U_dir = 1.
                V_dir = -1.

        # Move the particle in the direction of the highest velocity
        # First checking immediate neighbours, then checking diagonal neighbours.
        particle_dlon += U_dir * unbeach_U * particle.dt
        particle_dlat += V_dir * unbeach_V * particle.dt




def unbeachingBySampling(particle, fieldset, time):
    """Unbeaching particles by sampling the velocity field and moving in the direction of velocities!
    """
    # Measure the velocity field at the final particle location
    # Note, we've already handled the surface and bathymetry boundary conditions, so the particle is within the water column
    new_lon = particle.lon + particle_dlon
    new_lat = particle.lat + particle_dlat
    new_depth = particle.depth + particle_ddepth
    (vel_u, vel_v) = fieldset.UV[time, new_depth, new_lat, new_lon]  # noqa


    # If the future position of the particle's velocity is (nearly) zero, this means we will be on land
    if math.sqrt(vel_u**2 + vel_v**2) < 1e-14: 
        # Case 1: we move onto land -> most likely case because particles should be initialised in water!
        #if math.fabs(particle_dlon) > 1e-14 or math.fabs(particle_dlat) > 1e-14:
        #    # Case 2: We have moved onto land, so rather, don't move, but still allow for the vertical movement.
        #    particle_dlon = 0.
        #    particle_dlat = 0.

        ##Case 2: In the rare case that we start on land -> We haven't moved and our velocity is zero - we need to find which direction to move in.
        #else: # if math.fabs(particle_dlon) == 0. and math.fabs(particle_dlat) == 0.:
        displacement = 1./8. # Degree displacement to sample the velocity field

        unbeach_U = 1. / (1852. * 60. * math.cos(particle.lat * math.pi / 180.)) # Convert 1m/2s to degrees/s at the particle latitude in zonal direction
        unbeach_V = 1. / (1852. * 60.) # Convert 1m/2s to degrees/s in meridional direction

        (U_left, V_left) = fieldset.UV[time, new_depth, new_lat, new_lon - displacement]
        (U_right, V_right) = fieldset.UV[time, new_depth, new_lat, new_lon + displacement]
        (U_up, V_up) = fieldset.UV[time, new_depth, new_lat + displacement, new_lon]
        (U_down, V_down) = fieldset.UV[time, new_depth, new_lat - displacement, new_lon]

        # Find the direction of the highest velocity
        left = math.sqrt(U_left**2 + V_left**2)
        right = math.sqrt(U_right**2 + V_right**2)
        up = math.sqrt(U_up**2 + V_up**2)
        down = math.sqrt(U_down**2 + V_down**2)
        
        max_vel = 0.
        U_dir = 0.
        V_dir = 0.
        
        if left + right + up + down > 1e-14:
            # Annoyingly in C, we can't do python like comparisons, so be basic.
            # Let's just iterate over the four possibilities
            max_vel = left
            U_dir = -1.
            V_dir = 0.
            if max_vel < right:
                max_vel = right
                U_dir = 1.
                V_dir = 0.
            if max_vel < up:
                max_vel = up
                U_dir = 0.
                V_dir = 1.
            if max_vel < down:
                max_vel = down
                U_dir = 0.
                V_dir = -1.

        else: # In the case that all four directions are zero - i.e., we are stuck in a corner
            (U_left_up, V_left_up) = fieldset.UV[time, new_depth, new_lat + displacement, new_lon - displacement]
            (U_left_down, V_left_down) = fieldset.UV[time, new_depth, new_lat - displacement, new_lon - displacement]
            (U_right_up, V_right_up) = fieldset.UV[time, new_depth, new_lat + displacement, new_lon + displacement]
            (U_right_down, V_right_down) = fieldset.UV[time, new_depth, new_lat - displacement, new_lon + displacement]

            left_up = math.sqrt(U_left_up**2 + V_left_up**2)
            left_down = math.sqrt(U_left_down**2 + V_left_down**2)
            right_up = math.sqrt(U_right_up**2 + V_right_up**2)
            right_down = math.sqrt(U_right_down**2 + V_right_down**2)

            max_vel = left_up
            U_dir = -1.
            V_dir = 1.
            if max_vel < left_down:
                max_vel = left_down
                U_dir = 1.
                V_dir = -1.
            if max_vel < right_up:
                max_vel = right_up
                U_dir = 1.
                V_dir = 1.
            if max_vel < right_down:
                max_vel = right_down
                U_dir = 1.
                V_dir = -1.

        # Move the particle in the direction of the highest velocity
        # First checking immediate neighbours, then checking diagonal neighbours.
        particle_dlon += U_dir * unbeach_U * particle.dt
        particle_dlat += V_dir * unbeach_V * particle.dt


def initialise_algae_amount(particle, fieldset, time):
    """ Kernel to initialise the algae amount on a particle.
    This kernel is to be used for 'restarted' simulations. Because particle.time > 0 for restarted simulations,
    we need to assume that if a particle is below the ocean surface and it's algal amount is zero, then we need
    to give it an initial algal amount such that it is neutrally buoyant.
    """
    if particle.algae_amount == 0. and particle.depth > 0.:
        # For the particle to  be initialised as neutrally buoyant, the effective density of the particle must be equal to the seawater density
        initial_algae_amount = 0.5 * particle.plastic_diameter * (particle.seawater_density - particle.plastic_density) / ( fieldset.algae_cell_volume * (fieldset.biofilm_density - particle.seawater_density) )
        particle.algae_amount = max(initial_algae_amount, 0.)  # Ensure non-negative 

def NWS_kernel(particle, fieldset, time):
    """ remove particle if it exits box inside all meshes"""
    # handle each case separately (in case the or confuses the and...)
    if particle.lat > 72.0:
        particle.state = StatusCode.StopExecution
        particle.delete()
    elif particle.lat < 32.0:
        particle.state = StatusCode.StopExecution
        particle.delete()
    elif particle.lon > 14.0:
        particle.state = StatusCode.StopExecution
        particle.delete()
    elif particle.lon < -26.0:
        particle.state = StatusCode.StopExecution
        particle.delete()

def ARC_kernel(particle, fieldset, time):
    """ remove particle if it goes below specific latitude"""
    if particle.lat < 53. and particle.state == StatusCode.Success:
        particle.state = StatusCode.StopExecution

def depth_delete(particle, fieldset, time):
    """ Kernel to delete a particle if it goes below a certain depth. """
    if particle.depth > 71.:# and particle.state == StatusCode.Success:
        particle.state = StatusCode.StopExecution

def ARC_kernel_Copernicus(particle, fieldset, time):
    """ remove particle if it goes below specific latitude"""
    if particle.lat < 50. and particle.state == StatusCode.Success:
        particle.state = StatusCode.StopExecution