Merge pull request #1762 from knutfrode/dev

Fixed three bugs in OpenBerg: decreasing wave force coefficient from …

Author: knutfrode

opendrift/models/openberg.py

@@ -108,7 +108,7 @@ def ocean_force(iceb_vel, water_vel, Avo, Aho, rho_water, water_form_drag_coef,
         water_vel : Ocean current velocity
         rho_water : Water density
         Avo : Vertical area of the iceberg in the ocean (length x draft)
-        Aho : Horizontal area of the iceberg in contact with the ocean (width x draft)
+        Aho : Horizontal area of the iceberg in contact with the ocean (width x length)
         water_form_drag_coef : Co is the ocean form drag coefficient
         water_skin_drag_coef : Cdo is the ocean skin drag coefficient
     """
@@ -129,7 +129,7 @@ def wind_force(iceb_vel, wind_vel, Ava, Aha, wind_form_drag_coef, wind_skin_drag
         iceb_vel : Iceberg's velocity at time t
         wind_vel : Wind velocity
         Ava : Vertical area of the iceberg in the air (length x sail)
-        Aha : Horizontal area of the iceberg in contact with the air (width x sail)
+        Aha : Horizontal area of the iceberg in contact with the air (width x length)
         wind_form_drag_coef : Ca is the air form drag coefficient
         wind_skin_drag_coef : Cda is the air skin drag coefficient
     """
@@ -154,8 +154,8 @@ def wave_radiation_force(rho_water, wave_height, wave_direction, iceb_length):
         wave_direction : Wave direction
         iceb_length    : Iceberg's length
     """
-    F_wave_x = (0.5 * rho_water * wave_drag_coef * g * iceb_length * (wave_height / 2) ** 2 * np.sin(np.deg2rad(wave_direction)))
-    F_wave_y = (0.5 * rho_water * wave_drag_coef * g * iceb_length * (wave_height / 2) ** 2 * np.cos(np.deg2rad(wave_direction)))
+    F_wave_x = (0.25 * rho_water * wave_drag_coef * g * iceb_length * (wave_height / 2) ** 2 * np.sin(np.deg2rad(wave_direction)))
+    F_wave_y = (0.25 * rho_water * wave_drag_coef * g * iceb_length * (wave_height / 2) ** 2 * np.cos(np.deg2rad(wave_direction)))
     return np.array([F_wave_x, F_wave_y])
 
 
@@ -447,10 +447,10 @@ def advect_iceberg(self):
         sea_ice_thickness = self.environment.sea_ice_thickness
         sea_ice_conc = self.environment.sea_ice_area_fraction
         water_depth = self.environment.sea_floor_depth_below_sea_level
-        Avo = length * abs(draft)
-        Aho = width * abs(draft)
+        Avo = length * draft
+        Aho = width * length
         Ava = length * sail
-        Aha = width * sail
+        Aha = width * length
         Ai = sea_ice_thickness * length
 
         mass = width * (Ava + Avo) * rho_iceb * weight_coef

tests/models/test_openberg.py

@@ -36,13 +36,13 @@ def test_openberg_constant_forcing():
     o.set_config('environment:constant:y_sea_water_velocity', 0)
     o.set_config('environment:constant:x_wind', 0)
     o.set_config('environment:constant:y_wind', 0)
-    o.set_config('environment:constant:sea_surface_wave_significant_height', 3)
+    o.set_config('environment:constant:sea_surface_wave_significant_height', 2)
     o.set_config('environment:constant:sea_surface_wave_from_direction', 270)
     o.set_config('drift:horizontal_diffusivity', 0)
     o.set_config('drift:coriolis', False)
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
-    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.127, 3)
+    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.055, 3)
     np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.0, 3)
 
     # No current and wind, waves from east
@@ -57,7 +57,7 @@ def test_openberg_constant_forcing():
     o.set_config('drift:coriolis', False)
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
-    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 3.873, 3)
+    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 3.915, 3)
     np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.0, 3)
 
     # Northwards current, no wind, no Coriolis
@@ -96,7 +96,7 @@ def test_openberg_constant_forcing():
     o.set_config('drift:coriolis', True)
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
-    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.114, 3)
+    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.083, 3)
     np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.058, 3)
 
     # No current, eastwards wind
@@ -109,7 +109,7 @@ def test_openberg_constant_forcing():
     o.set_config('drift:coriolis', False)
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
-    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.107, 3)
+    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.075, 3)
     np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.000, 3)
 
     # No current, weaker eastwards wind
@@ -123,7 +123,7 @@ def test_openberg_constant_forcing():
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
     # TODO: Half wind gives only 25% of eastwards movement. Should be checked.
-    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.028, 3)
+    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4.022, 3)
     np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.000, 3)
 
     # No current, westwards wind
@@ -136,7 +136,7 @@ def test_openberg_constant_forcing():
     o.set_config('drift:coriolis', False)
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
-    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 3.892, 3)
+    np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 3.924, 3)
     np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.000, 3)
 
     # No current, northwards wind
@@ -150,7 +150,7 @@ def test_openberg_constant_forcing():
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
     np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4, 3)
-    np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.054, 3)
+    np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 60.039, 3)
 
     # No current, southwards wind
     o = OpenBerg(loglevel=50)
@@ -163,7 +163,7 @@ def test_openberg_constant_forcing():
     o.seed_elements(4, 60, time=datetime.now())
     o.run(steps=2)
     np.testing.assert_almost_equal(o.result.lon.isel(time=-1), 4, 3)
-    np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 59.946, 3)
+    np.testing.assert_almost_equal(o.result.lat.isel(time=-1), 59.962, 3)
 
 
 def test_openberg_norkyst():