Merge pull request #1264 from kthyng/zeta_seafloor_checks

Updates to ROMS zeta, seafloor checks, tests

opendrift/models/basemodel/__init__.py

@@ -695,8 +695,25 @@ def interact_with_seafloor(self):
             return
         if 'sea_floor_depth_below_sea_level' not in self.env.priority_list:
             return
+
+        if not hasattr(self, 'environment') or not hasattr(
+                self.environment, 'sea_surface_height'):
+            logger.warning('Seafloor check not being run because sea_surface_height is missing. '
+                           'This will happen the first time the function is run but if it happens '
+                           'subsequently there is probably a problem.')
+            return
+
+        # the shape of these is different than the original arrays
+        # because it is for active drifters
         sea_floor_depth = self.sea_floor_depth()
-        below = np.where(self.elements.z < -sea_floor_depth)[0]
+        sea_surface_height = self.sea_surface_height()
+
+        # Check if any elements are below sea floor
+        # But remember that the water column is the sea floor depth + sea surface height + parameter Dcrit
+        Dcrit = self.get_config('general:seafloor_action_dcrit')
+        ibelow = self.elements.z < -(sea_floor_depth + sea_surface_height + Dcrit)
+        below = np.where(ibelow)[0]
+
         if len(below) == 0:
             logger.debug('No elements hit seafloor.')
             return
@@ -706,8 +723,7 @@ def interact_with_seafloor(self):
             logger.debug('Lifting %s elements to seafloor.' % len(below))
             self.elements.z[below] = -sea_floor_depth[below]
         elif i == 'deactivate':
-            self.deactivate_elements(self.elements.z < -sea_floor_depth,
-                                     reason='seafloor')
+            self.deactivate_elements(ibelow, reason='seafloor')
             self.elements.z[below] = -sea_floor_depth[below]
         elif i == 'previous':  # Go back to previous position (in water)
             logger.warning('%s elements hit seafloor, '

opendrift/models/chemicaldrift.py

@@ -93,6 +93,7 @@ class ChemicalDrift(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': None},
         'y_sea_water_velocity': {'fallback': None},
+        'sea_surface_height': {'fallback': 0},
         'x_wind': {'fallback': 0},
         'y_wind': {'fallback': 0},
         'land_binary_mask': {'fallback': None},

opendrift/models/larvalfish.py

@@ -68,6 +68,7 @@ class LarvalFish(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         'sea_surface_wave_significant_height': {'fallback': 0},
         'x_wind': {'fallback': 0},
         'y_wind': {'fallback': 0},

opendrift/models/model_template.py

@@ -83,6 +83,7 @@ class ModelTemplate(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         'x_wind': {'fallback': 0, 'important': False},
         'y_wind': {'fallback': 0, 'important': False},
         'ocean_vertical_diffusivity': {'fallback': 0.02, 'important': False, 'profiles': True},

opendrift/models/oceandrift.py

@@ -69,6 +69,7 @@ class OceanDrift(OpenDriftSimulation):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         'x_wind': {'fallback': 0},
         'y_wind': {'fallback': 0},
         'upward_sea_water_velocity': {'fallback': 0},
@@ -158,6 +159,10 @@ def __init__(self, *args, **kwargs):
                 'enum': ['none', 'lift_to_seafloor', 'deactivate', 'previous'],
                 'description': '"deactivate": elements are deactivated; "lift_to_seafloor": elements are lifted to seafloor level; "previous": elements are moved back to previous position; "none"; seafloor is ignored.',
                 'level': CONFIG_LEVEL_ADVANCED},
+            'general:seafloor_action_dcrit': {'type': 'float', 'default': 0.0,
+                'min': 0.0, 'max': 10000, 'units': 'm',
+                'description': 'This parameter represents the amount of water left in a grid cell to keep it wet in a wet/dry simulation for numerical stability. The condition checked for seafloor_action is z < -(sea_floor_depth + sea_surface_height + `general:seafloor_action_dcrit`. It is 0 by default to assume that a wet/dry case is not being run, however, if it is and the correct value is not known 0.1 is a good default to use.',
+                'level': CONFIG_LEVEL_ADVANCED},
             'drift:truncate_ocean_model_below_m': {'type': 'float', 'default': None,
                 'min': 0, 'max': 10000, 'units': 'm',
                 'description': 'Ocean model data are only read down to at most this depth, and extrapolated below. May be specified to read less data to improve performance.',
@@ -411,13 +416,15 @@ def vertical_buoyancy(self):
             self.elements.z[in_ocean] = np.minimum(0,
                 self.elements.z[in_ocean] + self.elements.terminal_velocity[in_ocean] * self.time_step.total_seconds())
 
-        # check for minimum height/maximum depth for each particle
-        Zmin = -1.*self.environment.sea_floor_depth_below_sea_level
+        # check for minimum height/maximum depth for each particle accouting also for
+        # the sea surface height and the critical/min wet cell depth
+        Dcrit = self.get_config('general:seafloor_action_dcrit')
+        Zmin = -1.*(self.environment.sea_floor_depth_below_sea_level + self.environment.sea_surface_height + Dcrit)
 
         # Let particles stick to bottom
         bottom = np.where(self.elements.z < Zmin)
         if len(bottom[0]) > 0:
-            logger.debug('%s elements reached seafloor, set to bottom' % len(bottom[0]))
+            logger.debug('%s elements reached seafloor, interacting with bottom' % len(bottom[0]))
             self.interact_with_seafloor()
             self.bottom_interaction(Zmin)
 
@@ -489,8 +496,10 @@ def vertical_mixing(self, store_depths=False):
 
         dt_mix = self.get_config('vertical_mixing:timestep')
 
-        # minimum height/maximum depth for each particle
-        Zmin = -1.*self.environment.sea_floor_depth_below_sea_level
+        # minimum height/maximum depth for each particle accouting also for
+        # the sea surface height and the critical/min wet cell depth
+        Dcrit = self.get_config('general:seafloor_action_dcrit')
+        Zmin = -1.*(self.environment.sea_floor_depth_below_sea_level + self.environment.sea_surface_height + Dcrit)
 
         # Eventual model specific preparions
         self.prepare_vertical_mixing()
@@ -609,6 +618,7 @@ def vertical_mixing(self, store_depths=False):
                 self.elements.z[reflect] = -self.elements.z[reflect]
 
             # Reflect elements going below seafloor
+            # Zmin accounts for sea_surface_height
             bottom = np.where(np.logical_and(self.elements.z < Zmin, self.elements.moving == 1))
             if len(bottom[0]) > 0:
                 logger.debug('%s elements penetrated seafloor, lifting up' % len(bottom[0]))
@@ -629,7 +639,7 @@ def vertical_mixing(self, store_depths=False):
             # Let particles stick to bottom
             bottom = np.where(self.elements.z < Zmin)
             if len(bottom[0]) > 0:
-                logger.debug('%s elements reached seafloor, set to bottom' % len(bottom[0]))
+                logger.debug('%s elements reached seafloor, interacting with bottom' % len(bottom[0]))
                 self.interact_with_seafloor()
                 self.bottom_interaction(Zmin)
 

opendrift/models/openberg.py

@@ -59,6 +59,7 @@ class OpenBerg(OceanDrift):
     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},

opendrift/models/openhns.py

@@ -109,6 +109,7 @@ class OpenHNS(OceanDrift):
         'y_sea_water_velocity': {
             'fallback': None
         },
+        'sea_surface_height': {'fallback': 0},
         'x_wind': {
             'fallback': None
         },

opendrift/models/openoil/openoil.py

@@ -229,6 +229,7 @@ class OpenOil(OceanDrift):
         'y_wind': {
             'fallback': None
         },
+        'sea_surface_height': {'fallback': 0},
         'upward_sea_water_velocity': {
             'fallback': 0,
             'important': False

opendrift/models/pelagicegg.py

@@ -59,6 +59,7 @@ class PelagicEggDrift(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         'sea_surface_wave_significant_height': {'fallback': 0},
         'sea_ice_area_fraction': {'fallback': 0},
         'x_wind': {'fallback': 0},

opendrift/models/physics_methods.py

@@ -1114,6 +1114,25 @@ def sea_floor_depth(self):
                 env['sea_floor_depth_below_sea_level'].astype('float32')
         return sea_floor_depth
 
+    def sea_surface_height(self):
+        '''Sea surface height (positive/negative) for presently active elements'''
+
+        if hasattr(self, 'environment') and \
+                hasattr(self.environment, 'sea_surface_height'):
+            if len(self.environment.sea_surface_height) == \
+                    self.num_elements_active():
+                sea_surface_height = \
+                    self.environment.sea_surface_height
+        if 'sea_surface_height' not in locals():
+            env, env_profiles, missing = \
+                self.env.get_environment(['sea_surface_height'],
+                                     time=self.time, lon=self.elements.lon,
+                                     lat=self.elements.lat,
+                                     z=0*self.elements.lon, profiles=None)
+            sea_surface_height = \
+                env['sea_surface_height'].astype('float32')
+        return sea_surface_height
+
 
 def wind_drag_coefficient(windspeed):
     '''Large and Pond (1981), J. Phys. Oceanog., 11, 324-336.'''

opendrift/models/plastdrift.py

@@ -44,6 +44,7 @@ class PlastDrift(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0},
         'sea_surface_wave_stokes_drift_y_velocity': {'fallback': 0},
         'sea_surface_wave_significant_height': {'fallback': 0},

opendrift/models/radionuclides.py

@@ -75,6 +75,7 @@ class RadionuclideDrift(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': None},
         'y_sea_water_velocity': {'fallback': None},
+        'sea_surface_height': {'fallback': 0},
         'x_wind': {'fallback': 0},
         'y_wind': {'fallback': 0},
         'land_binary_mask': {'fallback': None},

opendrift/models/sealice.py

@@ -84,6 +84,7 @@ class SeaLice(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         # 'sea_surface_wave_significant_height': {'fallback': 0},
         # 'x_wind': {'fallback': 0},
         # 'y_wind': {'fallback': 0},

opendrift/models/sedimentdrift.py

@@ -45,6 +45,7 @@ class SedimentDrift(OceanDrift):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0},
         'upward_sea_water_velocity': {'fallback': 0},
         'x_wind': {'fallback': 0},
         'y_wind': {'fallback': 0},

opendrift/readers/basereader/consts.py

@@ -17,7 +17,7 @@
                       '(northwards if projection is lonlat/Mercator)'},
     'land_binary_mask': {'valid_min': 0, 'valid_max': 1,
                          'long_name': '1 is land, 0 is sea'},
-    'sea_floor_depth_below_sea_level': {'valid_min': -10, 'valid_max': 12000,
+    'sea_floor_depth_below_sea_level': {'valid_min': -20, 'valid_max': 12000,
                          'long_name': 'Depth of seafloor'},
     'ocean_vertical_diffusivity': {'valid_min': 0, 'valid_max': 1}}