Merge pull request #1289 from kthyng/fix_z_rho_d_crit

fix z_rho and remove dcrit

opendrift/models/basemodel/__init__.py

@@ -701,22 +701,24 @@ def interact_with_seafloor(self):
         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. '
+        if not hasattr(self, 'environment'):
+            logger.warning('Seafloor check not being run because environment is missing. '
                            'This will happen the first time the function is run but if it happens '
                            'subsequently there is probably a problem.')
             return
 
+        if not hasattr(self.environment, 'sea_surface_height'):
+            logger.warning('Seafloor check not being run because sea_surface_height is missing. ')
+            return
+
         # the shape of these is different than the original arrays
         # because it is for active drifters
         sea_floor_depth = self.sea_floor_depth()
         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)
+        # But remember that the water column is the sea floor depth + sea surface height
+        ibelow = self.elements.z < -(sea_floor_depth + sea_surface_height)
         below = np.where(ibelow)[0]
 
         if len(below) == 0:

opendrift/models/oceandrift.py

@@ -160,10 +160,6 @@ 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.',
@@ -421,9 +417,8 @@ def vertical_buoyancy(self):
                 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 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)
+        # the sea surface height
+        Zmin = -1.*(self.environment.sea_floor_depth_below_sea_level + self.environment.sea_surface_height)
 
         # Let particles stick to bottom
         bottom = np.where(self.elements.z < Zmin)
@@ -500,9 +495,8 @@ def vertical_mixing(self, store_depths=False):
         dt_mix = self.get_config('vertical_mixing:timestep')
 
         # 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)
+        # the sea surface height
+        Zmin = -1.*(self.environment.sea_floor_depth_below_sea_level + self.environment.sea_surface_height)
 
         # Eventual model specific preparions
         self.prepare_vertical_mixing()

opendrift/readers/reader_ROMS_native.py

@@ -122,7 +122,8 @@ def __init__(self, filename=None, name=None, gridfile=None, standard_name_mappin
             'uwind': 'x_wind',
             'vwind': 'y_wind',
             'Uwind': 'x_wind',
-            'Vwind': 'y_wind'}
+            'Vwind': 'y_wind',
+            }
 
         # Add user provided variable mappings
         self.ROMS_variable_mapping.update(standard_name_mapping)
@@ -483,8 +484,13 @@ def get_variables(self, requested_variables, time=None,
             # Transform z_rho to match convention of z.
             z_rho -= np.asarray(zeta)[np.newaxis]
 
+            # Check for positive values in z_rho
+            # if there are any, nan them out since they are above
+            # the surface and therefore the cell is dry
+            # this should only come up for wet/dry simulations
             if (np.nanmax(z_rho) > 0).any():
-                logger.warning('z_rho is positive, but should be negative or 0.')
+                logger.info(f'z_rho had positive values that are now nans.')
+                z_rho[z_rho>0] = np.nan
 
             # Element indices must be relative to extracted subset
             indx_el = np.clip(indx_el - indx.min(), 0, z_rho.shape[2]-1)
@@ -543,7 +549,7 @@ def get_mask(mask_name, imask, masks_store):
                 variables[par] = var[itzxy]
             else:
                 raise Exception('Wrong dimension of variable: ' +
-                                self.variable_mapping[par])
+                                self.ROMS_variable_mapping[par])
 
             # If Xarray, load in so can be used in future loop iterations too
             variables[par] = np.asarray(variables[par])