From f472b3f73f5a1a8a36e05b57b33b79bfcd4586ae Mon Sep 17 00:00:00 2001 From: Achref Othmani Date: Fri, 28 Jun 2024 09:53:24 +0200 Subject: [PATCH] #1341 updated icebergs module --- examples/example_openberg_new.py | 18 +- opendrift/models/openberg.py | 551 +++++++++++++++++++++++++------ 2 files changed, 456 insertions(+), 113 deletions(-) diff --git a/examples/example_openberg_new.py b/examples/example_openberg_new.py index f810e8909..ed7882bb8 100755 --- a/examples/example_openberg_new.py +++ b/examples/example_openberg_new.py @@ -1,20 +1,24 @@ #!/usr/bin/env python """ -Ice bergs (openberg) +Icebergs (openberg) ==================== """ -from datetime import datetime, timedelta from opendrift.models.openberg import OpenBerg - -o = OpenBerg() +from datetime import datetime,timedelta + +o = OpenBerg(add_stokes_drift = True, + wave_rad = True, + grounding = False, + vertical_profile = False, + melting = False, + choose_melting = {"wave": True, "lateral": True, "basal": True}) o.add_readers_from_list([ 'https://thredds.met.no/thredds/dodsC/cmems/topaz6/dataset-topaz6-arc-15min-3km-be.ncml', 'https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/ncep_global/NCEP_Global_Atmospheric_Model_best.ncd']) -o.seed_elements(lon=28, lat=77, time=datetime.now(), number=1000, radius=1000) - -o.run(duration=timedelta(days=2)) +o.seed_elements(lon= -57,lat= 69,time=datetime.now(), number=1000, radius=1000) +o.run(duration=timedelta(days=2)) o.plot(fast=True) diff --git a/opendrift/models/openberg.py b/opendrift/models/openberg.py index 9de178abf..2a959f4d9 100644 --- a/opendrift/models/openberg.py +++ b/opendrift/models/openberg.py @@ -1,106 +1,445 @@ -# This file is part of OpenDrift. -# -# OpenDrift is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, version 2 -# -# OpenDrift is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with OpenDrift. If not, see . -# -# Copyright 2015, 2023, Knut-Frode Dagestad, MET Norway -# Copyright 2023, Achref Othmani, NERSC, Norway - - -import numpy as np -import logging; logger = logging.getLogger(__name__) -from opendrift.models.oceandrift import OceanDrift, Lagrangian3DArray - -class IcebergObj(Lagrangian3DArray): - """Extending LagrangianArray with relevant properties for an Iceberg""" - - variables = Lagrangian3DArray.add_variables([ - ('sail', {'dtype': np.float32, # Sail of Iceberg (part above waterline ) - 'units': 'm', - 'default': 10}), - ('draft', {'dtype': np.float32, # Draft of Iceberg (part below waterline) - 'units': 'm', - 'default': 90}), - ('length', {'dtype': np.float32, # length of Iceberg - 'units': 'm', - 'default': 100}), - ('width', {'dtype': np.float32, # width of Iceberg - 'units': 'm', - 'default': 30}), - ('weight_coeff', {'dtype': np.float32, # Proportion of the mass of rectangular Iceberg, weight_coeff = 1.0 if 100%, 0.5 if 50% and 0.35 if 35% - 'units': '1', - 'default': 1}), - ('water_drag_coeff', {'dtype': np.float32, # cdo - 'units': '1', - 'default': 0.25}), - ('wind_drag_coeff', {'dtype': np.float32, # cda - 'units': '1', - 'default': 0.7}) - ]) - - -class OpenBerg(OceanDrift): - - ElementType = IcebergObj - - # Specify which environment variables (e.g. wind, waves, currents...) - # are needed/required by the present model, to be used for updating - # the element properties (including propagation). - - required_variables = { - 'x_sea_water_velocity': {'fallback': None}, - 'y_sea_water_velocity': {'fallback': None}, - 'sea_surface_height': {'fallback': 0}, - 'sea_floor_depth_below_sea_level': {'fallback': 10000}, - 'x_wind': {'fallback': None}, - 'y_wind': {'fallback': None}, - #'sea_surface_wave_significant_height': {'fallback': None}, # Needed for melting - 'land_binary_mask': {'fallback': None}, - } - - # Configuration - def __init__(self, *args, **kwargs): - - # The constructor of parent class must always be called - # to perform some necessary common initialisation tasks: - super(OpenBerg, self).__init__(*args, **kwargs) - - - def thermodynamics(self): - pass # not yet implemented - - def update(self): - """Update positions and properties of icebergs""" - - self.thermodynamics() - - # Constants - rho_water = 1027 - rho_air = 1.293 - rho_ice = 917 - rho_iceb = 900 - - - # Areas exposed - Ao = abs(self.elements.draft) * self.elements.length # Area_wet - Aa = self.elements.sail * self.elements.length # Area_dry - - # See ACCIBERG presentation - # https://docs.google.com/presentation/d/1O5C2v7PA3PW8a93IAGU-aS6BSKt3s-Fw/edit#slide=id.p1 - k = rho_air*self.elements.wind_drag_coeff*Aa / (rho_water*self.elements.water_drag_coeff*Ao) - f = np.sqrt(k)/(1+np.sqrt(k)) - vx = (1-f)*self.environment.x_sea_water_velocity + f*self.environment.x_wind - vy = (1-f)*self.environment.y_sea_water_velocity + f*self.environment.y_wind - self.update_positions(vx, vy) - - # Grounding - self.deactivate_elements(self.elements.draft > self.environment.sea_floor_depth_below_sea_level, reason='Grounded iceberg') +# This file is part of OpenDrift. +# +# OpenDrift is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, version 2 +# +# OpenDrift is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with OpenDrift. If not, see . +# +# Copyright 2015, 2023, Knut-Frode Dagestad, MET Norway +# Copyright 2023, Achref Othmani & 2024, Lenny Hucher, NERSC, Norway + +import numpy as np +import logging; logger = logging.getLogger(__name__) +from opendrift.models.oceandrift import OceanDrift, Lagrangian3DArray +from opendrift.models.physics_methods import PhysicsMethods +from scipy.integrate import solve_ivp + + + +# Constants +rho_water = 1027 +rho_air = 1.293 +rho_ice = 917 +rho_iceb = 900 +g = 9.81 +csi = 1 # Sea ice coefficient of resistance +wave_drag_coef = 0.3 + + +class IcebergObj(Lagrangian3DArray): + """ Extending Lagrangian3DArray with relevant properties for an Iceberg """ + + variables = Lagrangian3DArray.add_variables([ + ('sail', {'dtype': np.float32, # Sail of iceberg (part above waterline ) + 'units': 'm', + 'default': 10}), + ('draft', {'dtype': np.float32, # Draft of iceberg (part below waterline) + 'units': 'm', + 'default': 90}), + ('length', {'dtype': np.float32, # length of iceberg + 'units': 'm', + 'default': 100}), + ('width', {'dtype': np.float32, # width of iceberg + 'units': 'm', + 'default': 30}), + ('weight_coeff', {'dtype': np.float32, # Relative to the shape of iceberg (e.g. 1 for tabular; 0.3 for pinnacle) + 'units': '1', + 'default': 1}), + ('water_drag_coeff', {'dtype': np.float32, # Ocean drag coeff. + 'units': '1', + 'default': 0.25}), + ('wind_drag_coeff', {'dtype': np.float32, # Wind drag coeff. + 'units': '1', + 'default': 0.7}), + ("iceb_x_velocity", {"dtype": np.float32, # Iceberg velocity in the x-direction + "units": "m/s", + "default": 0.0}), + ("iceb_y_velocity", {"dtype": np.float32, # Iceberg velocity in the y-direction + "units": "m/s", + "default": 0.0}), + ]) + + + +# Define the functions needed +def water_drag(iceb_vel, water_vel, Ao, rho_water, water_drag_coef): + """ Ocean force + Args: + iceb_vel : Iceberg velocity at time t + water_vel : Ocean current velocity + Ao : Iceberg's area exposed to the ocean current(length x draft) + rho_water : Water density + water_drag_coef : Co is the drag coefficient applied on the iceberg's draft + Returns: + Ocean force + """ + vxo, vyo = water_vel[0], water_vel[1] # water velocity + x_vel, y_vel = iceb_vel[0], iceb_vel[1] # iceberg velocity + rel_water_x_vel = vxo - x_vel + rel_water_y_vel = vyo - y_vel + rel_water_norm = np.sqrt(rel_water_x_vel**2 + rel_water_y_vel**2) + F_ocean_x = (0.5 * rho_water * water_drag_coef * Ao * rel_water_norm * rel_water_x_vel) + F_ocean_y = (0.5 * rho_water * water_drag_coef * Ao * rel_water_norm * rel_water_y_vel) + return np.array([F_ocean_x, F_ocean_y]) + + +def wind_drag(iceb_vel, wind_vel, Aa, wind_drag_coef): + """ Wind force + Args: + iceb_vel : Iceberg velocity at time t + wind_vel : Wind velocity + Aa : Iceberg's area exposed to the wind (length x sail) + rho_air : Air density + wind_drag_coef : Ca is the drag coefficient applied on the iceberg's sail + Returns: + Wind force + """ + vxa, vya = wind_vel[0], wind_vel[1] # wind velocity + x_vel, y_vel = iceb_vel[0], iceb_vel[1] # iceberg velocity + rel_wind_x_vel = vxa - x_vel + rel_wind_y_vel = vya - y_vel + rel_wind_norm = np.sqrt(rel_wind_x_vel**2 + rel_wind_y_vel**2) + F_wind_x = 0.5 * rho_air * wind_drag_coef * Aa * rel_wind_norm * rel_wind_x_vel + F_wind_y = 0.5 * rho_air * wind_drag_coef * Aa * rel_wind_norm * rel_wind_y_vel + return np.array([F_wind_x, F_wind_y]) + + +def wave_radiation_force(rho_water, wave_height, wave_direction, iceb_length): + """ Wave radiation force + Args: + iceb_vel : Iceberg velocity at time t + rho_water : Water density + wave_height : Wave significant height + wave_direction : Wave direction + iceb_length : Iceberg's length + Returns: + Wave radiation force + """ + F_wave_x = (0.5 * rho_water * wave_drag_coef * g * iceb_length * (wave_height / 2) ** 2 * np.sin(np.deg2rad((wave_direction + 180) % 360))) + F_wave_y = (0.5 * rho_water * wave_drag_coef * g * iceb_length * (wave_height / 2) ** 2 * np.cos(np.deg2rad((wave_direction + 180) % 360))) + return np.array([F_wave_x, F_wave_y]) + + +def advect_iceberg_no_acc(f, water_vel, wind_vel): + """ Advect iceberg without acceleration + Args: + f : Wind drift factor + water_vel : Ocean current velocity + wind_vel : Wind velocity + + Returns: + Iceberg velocity without acceleration + """ + vxo, vyo = water_vel[0], water_vel[1] + vxa, vya = wind_vel[0], wind_vel[1] + no_acc_model_x = (1 - f) * vxo + f * vxa + no_acc_model_y = (1 - f) * vyo + f * vya + V = np.array([no_acc_model_x, no_acc_model_y]) + if not np.isfinite(V).all(): + logger.error("Infinite value in iceberg's velocity without acceleration: Please check the wind drift factor f ") + return V + + +def sea_ice_force(iceb_vel, sea_ice_conc, sea_ice_thickness, sea_ice_vel, iceb_width, sum_force): + """ Sea ice force + Args: + iceb_vel : Iceberg velocity at time t + sea_ice_conc : Sea ice concentration [%] + sea_ice_thickness : Sea ice thickness [m] + sea_ice_vel : Sea ice velocity [m/s] + rho_ice : Sea ice density + iceb_width : Iceberg width + sum_force : Effect of all other forces exerted on the iceberg (apart from the sea ice force) + Returns: + Sea ice force + """ + ice_x, ice_y = sea_ice_vel + x_vel, y_vel = iceb_vel[0], iceb_vel[1] + diff_vel = np.sqrt((ice_x - x_vel) ** 2 + (ice_y - y_vel) ** 2) + force_x, force_y = sum_force + F_ice_x = np.zeros_like(x_vel) + F_ice_y = np.zeros_like(y_vel) + F_ice_x = (0.5 * (rho_ice * csi * sea_ice_thickness * iceb_width) * diff_vel * (ice_x - x_vel)) + F_ice_y = (0.5 * (rho_ice * csi * sea_ice_thickness * iceb_width) * diff_vel * (ice_y - y_vel)) + F_ice_x[sea_ice_conc <= 0.15] = 0 + F_ice_y[sea_ice_conc <= 0.15] = 0 + F_ice_x[sea_ice_conc >= 0.9] = -force_x[sea_ice_conc >= 0.9] + F_ice_y[sea_ice_conc >= 0.9] = -force_y[sea_ice_conc >= 0.9] + return np.array([F_ice_x, F_ice_y]) + + +def melwav(iceb_length, iceb_width, x_wind, y_wind, sst, sea_ice_conc, dt): + """ Update the iceberg's dimensions (length and width) due to wave erosion + Parameterization from Keghouche et al. 2009 + Args: + iceb_length : Iceberg's length + iceb_width : Iceberg's width + x_wind : Wind speed in the x-direction + y_wind : Wind speed in the y-direction + sst : Sea surface temperature + sea_ice_conc : Sea ice concentration + dt : Timestep of the simulation + """ + Ss = -5 + np.sqrt(32 + 2 * np.sqrt(x_wind**2 + y_wind**2)) + Vsst = (1 / 6.0) * (sst + 2) * Ss + Vwe = Vsst * 0.5 * (1 + np.cos(np.pi * sea_ice_conc**3)) / 86400 + new_iceb_length = np.zeros_like(iceb_length) + new_iceb_width = np.zeros_like(iceb_width) + new_iceb_length[iceb_length != 0] = (iceb_length[iceb_length != 0] - Vwe[iceb_length != 0] * dt) + new_iceb_width[iceb_length != 0] = (iceb_width[iceb_length != 0] / iceb_length[iceb_length != 0] * new_iceb_length[iceb_length != 0]) + new_iceb_length[new_iceb_length < 0] = 0 + new_iceb_width[new_iceb_width < 0] = 0 + return new_iceb_length, new_iceb_width + +def mellat(iceb_length, iceb_width, tempib, salnib, dt): + """ Update the iceberg's dimensions (length and width) due to lateral melting + Parameterization from Keghouche et al. 2009 + Args: + iceb_length : Iceberg's length + iceb_width : Iceberg's width + tempib : Water temperature + salnib : Water salinity + dt : Timestep of the simulation + """ + TfS = -0.036 - 0.0499 * salnib - 0.000112 * salnib**2 + Tfp = TfS * np.exp(-0.19 * (tempib - TfS)) + deltaT = tempib - Tfp + deltaT = np.concatenate([2.78 * deltaT, 0.47 * deltaT**2], axis=0) + sumVb = np.nansum(deltaT, axis=0) + dx = sumVb / 365 / 86400 * dt # Convert sumVb from (meter/year) to (meter per second) + + new_iceb_length = np.zeros_like(iceb_length) + new_iceb_width = np.zeros_like(iceb_width) + new_iceb_length[iceb_length != 0] = (iceb_length[iceb_length != 0] - 2 * dx[iceb_length != 0]) + new_iceb_width[iceb_length != 0] = (iceb_width[iceb_length != 0] / iceb_length[iceb_length != 0] * new_iceb_length[iceb_length != 0]) + new_iceb_length[new_iceb_length < 0] = 0 + new_iceb_width[new_iceb_width < 0] = 0 + return new_iceb_length, new_iceb_width + +def melbas(iceb_draft, iceb_sail, iceb_length, salnib, tempib, x_water_vel, y_water_vel, x_iceb_vel, y_iceb_vel, dt): + """ Update the iceberg's dimensions (draft and sail) due to forced convection + Parameterization from Keghouche et al. 2009 + """ + # Temperature at the base layer of the iceberg + absv = np.sqrt(((x_water_vel - x_iceb_vel) ** 2 + (y_water_vel - y_iceb_vel) ** 2)) + TfS = -0.036 - 0.0499 * salnib - 0.000112 * salnib**2 + Tfp = TfS * 2.71828 ** (-0.19 * (tempib - TfS)) + deltat = tempib - Tfp + Vf = 0.58 * absv**0.8 * deltat / (iceb_length**0.2) + Vf = Vf / 86400 # converted to m/s + # Update the draft + new_iceb_draft = np.zeros_like(iceb_draft) + new_iceb_draft[iceb_draft != 0] = (abs(iceb_draft[iceb_draft != 0]) - Vf[iceb_draft != 0] * dt) + # Melt at base of the iceberg + new_iceb_draft[iceb_draft < 0] = 0 + return new_iceb_draft, iceb_sail + + +class OpenBerg(OceanDrift): + + ElementType = IcebergObj + + required_variables = { + "x_sea_water_velocity": {"fallback": None, "profiles": True}, + "y_sea_water_velocity": {"fallback": None, "profiles": True}, + "sea_floor_depth_below_sea_level": {"fallback": 10000}, + "x_wind": {"fallback": 0, "important": False}, + "y_wind": {"fallback": 0, "important": False}, + "sea_surface_wave_significant_height": {"fallback": 0}, + "sea_surface_wave_from_direction": {"fallback": 0}, + "sea_surface_wave_stokes_drift_x_velocity": {"fallback": 0}, + "sea_surface_wave_stokes_drift_y_velocity": {"fallback": 0}, + "sea_water_temperature": {"fallback": 2, "profiles": True}, + "sea_water_salinity": {"fallback": 35, "profiles": True}, + "sea_ice_area_fraction": {"fallback": 0}, + "sea_ice_thickness": {"fallback": 0}, + "sea_ice_x_velocity": {"fallback": 0, "important": False}, + "sea_ice_y_velocity": {"fallback": 0, "important": False}, + "land_binary_mask": {"fallback": None}, + } + + def get_profile_masked(self, variable): + """ Apply a mask to extract data from the surface down to the iceberg's draft """ + draft = self.elements.draft + profile = self.environment_profiles[variable] + z = self.environment_profiles["z"] + if len(z) == 1: + if z == np.array([None]): # Convert None values to zeros for surface current + z = np.zeros_like(z) + mask = draft[:, np.newaxis] < -z + mask = mask.T + mask[np.argmax(mask, axis=0), np.arange(mask.shape[1])] = False + return np.ma.masked_array(profile, mask, fill_value=np.nan) + + def get_basal_env(self, variable): + """ Get the basal layer of the variable for the icebergs """ + profile = self.get_profile_masked(variable) + last = np.argmin(np.logical_not(profile.mask), axis=0) - 1 + return profile[last, np.arange(profile.shape[1])] + + + # Configuration + def __init__( + self, + add_stokes_drift: bool = True, + wave_rad: bool = True, + grounding: bool = False, + vertical_profile: bool = False, + melting: bool = False, + choose_melting: dict[bool] = {"wave": True, "lateral": True, "basal": True}, + *args, + **kwargs, + ): + + # The constructor of parent class must always be called + # to perform some necessary common initialisation tasks: + super(OpenBerg, self).__init__(*args, **kwargs) + self.wave_rad = wave_rad + self.add_stokes_drift = add_stokes_drift + self.grounding = grounding + self.vertical_profile = vertical_profile + self.melting = (melting) #Activate the melting + self.choose_melting = (choose_melting) #Choose how it melts + + def advect_iceberg(self, rho_water, stokes_drift=True, wave_rad=True, grounding=False, vertical_profile=False): + draft = self.elements.draft + length = self.elements.length + Ao = abs(draft) * length ### Area_wet + Aa = self.elements.sail * length ### Area_dry + mass = self.elements.width * (Aa + Ao) * rho_iceb + k = (rho_air * self.elements.wind_drag_coeff * Aa / (rho_water * self.elements.water_drag_coeff * Ao)) + f = np.sqrt(k) / (1 + np.sqrt(k)) + + if not vertical_profile: + water_vel = np.array([self.environment.x_sea_water_velocity + int(stokes_drift) * self.environment.sea_surface_wave_stokes_drift_x_velocity, + self.environment.y_sea_water_velocity + int(stokes_drift) * self.environment.sea_surface_wave_stokes_drift_y_velocity]) + else: + logger.info("Calculating the Depth Integrated Currents ...") + uprof = self.get_profile_masked("x_sea_water_velocity") + vprof = self.get_profile_masked("y_sea_water_velocity") + + z = self.environment_profiles["z"] + thickness = -(z[1:] - z[:-1]).reshape((-1, 1)) + mask = uprof.mask + uprof_mean_inter = (uprof[1:] + uprof[:-1]) / 2 + vprof_mean_inter = (vprof[1:] + vprof[:-1]) / 2 + mask = mask[:-1] + thickness_reshaped = np.tile(thickness, (1, mask.shape[1])) + thickness_reshaped[mask] = np.nan + umean = np.nansum(thickness_reshaped * uprof_mean_inter, axis=0) / np.nansum(thickness_reshaped, axis=0) + vmean = np.nansum(thickness_reshaped * vprof_mean_inter, axis=0) / np.nansum(thickness_reshaped, axis=0) + water_vel = np.array([umean, vmean]) + + water_depth = self.environment.sea_floor_depth_below_sea_level + wind_vel = np.array([self.environment.x_wind, self.environment.y_wind]) + wave_height = self.environment.sea_surface_wave_significant_height + wave_direction = self.environment.sea_surface_wave_from_direction + sea_ice_conc = self.environment.sea_ice_area_fraction + sea_ice_thickness = self.environment.sea_ice_thickness + sea_ice_vel = np.array([self.environment.sea_ice_x_velocity, self.environment.sea_ice_y_velocity]) + + def dynamic(t,iceb_vel, water_vel, wind_vel, wave_height, wave_direction, Ao, + Aa, rho_water, water_drag_coef, wind_drag_coef, iceb_length, mass): + """ Function required by solve_ivp. The t and iceb_vel parameters are required by solve_ivp, shouldn't be deleted """ + sum_force = (water_drag(iceb_vel, water_vel, Ao, rho_water, water_drag_coef) + + wind_drag(iceb_vel,wind_vel, Aa, wind_drag_coef) + + int(wave_rad) * wave_radiation_force(rho_water, wave_height, wave_direction, iceb_length)) + + sum_force = sum_force + sea_ice_force(iceb_vel, sea_ice_conc, sea_ice_thickness, sea_ice_vel, self.elements.width, sum_force) + return 1 / mass * sum_force + + V0 = advect_iceberg_no_acc(f, water_vel, wind_vel) # Approximation of the solution of the dynamic equation for the iceberg velocity + V0[:, sea_ice_conc >= 0.9] = sea_ice_vel[:, sea_ice_conc >= 0.9] # With this criterium, the iceberg moves with the sea ice + V0 = V0.flatten() + hwall = draft - water_depth + grounded = np.logical_and(hwall >= 0, grounding) + if any(grounded) and grounding: + logger.info(f"Grounding condition : Icebergs grounded = {len(hwall[hwall>0])}, hwall={np.round(hwall[hwall>0],3)} meters") + + sol = solve_ivp(dynamic, [0, self.time_step.total_seconds()], V0, + args=(water_vel, wind_vel, wave_height, wave_direction, + Ao, Aa, rho_water, + self.elements.water_drag_coeff, + self.elements.wind_drag_coeff, + self.elements.length, + mass), + vectorized=True, + t_eval=np.array([self.time_step.total_seconds()]), + ) + + V = sol.y.reshape((2, -1)) + Vx, Vy = V[0], V[1] + Vx[grounded] = 0 + Vy[grounded] = 0 + self.update_positions(Vx, Vy) + self.elements.iceb_x_velocity, self.elements.iceb_y_velocity = Vx, Vy + + def melt(self): + """ Activate melting """ + x_wind = self.environment.x_wind + y_wind = self.environment.y_wind + uoib = self.get_basal_env("x_sea_water_velocity") + voib = self.get_basal_env("y_sea_water_velocity") + T_profile = self.environment_profiles["sea_water_temperature"] + S_profile = self.environment_profiles["sea_water_salinity"] + Tn = self.get_basal_env("sea_water_temperature") + Sn = self.get_basal_env("sea_water_salinity") + sea_ice_conc = self.environment.sea_ice_area_fraction + + # Wave melting + if self.choose_melting["wave"]: + self.elements.length, self.elements.width = melwav(self.elements.length, self.elements.width, x_wind, y_wind, T_profile[0], sea_ice_conc, self.time_step.total_seconds()) + # Lateral melting + if self.choose_melting["lateral"]: + self.elements.length, self.elements.width = mellat(self.elements.length, self.elements.width, T_profile, S_profile, self.time_step.total_seconds()) + # Basal melting + if self.choose_melting["basal"]: + self.elements.draft, self.elements.sail = melbas(self.elements.draft, self.elements.sail, self.elements.length, Sn, Tn, uoib, voib, self.elements.iceb_x_velocity, self.elements.iceb_y_velocity, self.time_step.total_seconds()) + + # Deactivate elements less than 1 meter + self.deactivate_elements(self.elements.draft < 1, "Iceberg melted") + self.deactivate_elements(self.elements.length < 1, "Iceberg melted") + self.deactivate_elements(self.elements.width < 1, "Iceberg melted") + self.deactivate_elements(self.elements.sail < 1, "Iceberg melted") + + + def roll_over(self, rho_water): + """ Iceberg's stability criterium: Parameterization from Keghouche et al. 2009 with a correction from Wagner et al. 2017 """ + L = self.elements.length + W = self.elements.width + H = self.elements.draft + self.elements.sail + alpha = rho_iceb / rho_water + crit = np.sqrt(6 * alpha * (1 - alpha)) + W, L = np.min([L, W], axis=0), np.max([L, W], axis=0) + mask = (W / H) < crit + if any(mask): + logger.info(f"Rolling over : {np.sum(mask)} icebergs ...") + nL, nW, nH = (np.max([L[mask], H[mask]], axis=0), np.min([L[mask], H[mask]], axis=0), W[mask]) + L[mask], W[mask], H[mask] = nL, nW, nH + depthib = H * alpha + sailib = H - depthib + self.elements.length = L + self.elements.width = W + self.elements.sail = sailib + self.elements.draft = depthib + + def prepare_run(self): + self.profiles_depth = self.elements_scheduled.draft.max() + logger.info(f"Icebergs max draft is : {self.profiles_depth} meters") + + def update(self): + """ Update positions and properties of particles """ + T = self.environment.sea_water_temperature + S = self.environment.sea_water_salinity + rho_water = PhysicsMethods.sea_water_density(T, S) + self.roll_over(rho_water) + if self.melting: + self.melt() + self.advect_iceberg(rho_water, self.add_stokes_drift, self.wave_rad, self.grounding, self.vertical_profile) \ No newline at end of file