[run-ex] ROMS native reader: Copied newest version of module depth from roppy package, supporting now Vstretching=5. Manually merged in earlier adaptations for zeta made by Kristen Thyng.

Author: knutfrode

history.rst

@@ -204,7 +204,7 @@ History
 * Using now product_id instead of OPeNDAP URL for CMEMS datasets, and using copernicusmarineclient through new reader_copernicusmarine. username/password can be stored in netrc-file with machine name equal to *copernicusmarine* or *nrt.cmems-du.eu*
 * Model property reguired_profiles_z_range is now replaced with config setting drift:profile_depth, and profiles are retrieved from surface to this depth. profiles_depth is now input parameter to get_environment, and not anymore a property of Environment class. prepare_run must now always call prepare_run of parent class, since profile_depth is copied to object in basemodel.prepare_run
 * get_variables_along_trajectory now also takes depth (z) as input parameter
-* updates to wetting/drying in ROMS reader (Kristin Thyng)
+* updates to wetting/drying in ROMS reader (Kristen Thyng)
 * Fill value in output netCDF files is now set to NaN for floats and -999 for integers
 * Moving basereader.prepare() to variables.prepare(), as the former was overriding structured.prepare() due to multiple inheritance, and thus config *drift:max_speed* was not applied if config setting was made after reader was added. Also increasing *drift:max_speed* of OceanDrift from 1 to 2m/s
 * Leeway model now allows capsizing (and un-capsizing for backwards runs), with given probability and reduction of leeway coefficients when wind exceeds given threshold
@@ -217,8 +217,8 @@ History
 
 2024-04-02 / Release v1.11.2
 ----------------------------
-* Proper handling of sea_surface_height implemented by Kristin Thyng. All subclasses of OceanDrift now have `sea_surface_height` (default 0) as new parameter. z=0 is always sea surface (including sea_surface_height), and seafloor is now where z = -(sea_floor_depth + sea_surface_height)
-* Improvements of ROMS reader by Kristin Thyng:
+* Proper handling of sea_surface_height implemented by Kristen Thyng. All subclasses of OceanDrift now have `sea_surface_height` (default 0) as new parameter. z=0 is always sea surface (including sea_surface_height), and seafloor is now where z = -(sea_floor_depth + sea_surface_height)
+* Improvements of ROMS reader by Kristen Thyng:
 
   * Roppy-method `sdepth` (used by ROMS reader) now accounts for `sea_surface_height` (zeta).
   * Improved handling of rotation of vectors.

opendrift/readers/reader_ROMS_native.py

@@ -142,7 +142,7 @@ def drop_non_essential_vars_pop(ds):
             dropvars = [v for v in ds.variables if v not in
                         list(self.ROMS_variable_mapping.keys()) +
                         ['ocean_time', 'time', 'bulk_time', 's_rho',
-                         'Cs_r', 'Cs_rho', 'hc', 'angle', 'Vtransform']
+                         'Cs_r', 'Cs_rho', 'hc', 'angle', 'Vtransform', 'Vstretching']
                         and v[0:3] not in ['lon', 'lat', 'mas']]
             logger.debug('Dropping variables: %s' % dropvars)
             ds = ds.drop_vars(dropvars)
@@ -181,6 +181,13 @@ def drop_non_essential_vars_pop(ds):
             self.Vtransform = 1
         self.Vtransform = np.asarray(self.Vtransform)
 
+        if 'Vstretching' in self.Dataset.variables:
+            self.Vstretching = self.Dataset.variables['Vstretching'].data  # scalar
+        else:
+            logger.warning('Vstretching not found, using 1')
+            self.Vstretching = 1
+        self.Vstretching = np.asarray(self.Vstretching)
+
         if 's_rho' not in self.Dataset.variables:
             dimensions = 2
         else:
@@ -508,7 +515,8 @@ def get_variables(self, requested_variables, time=None,
                 Htot = self.sea_floor_depth_below_sea_level
                 zeta = self.zeta[indxTime]
                 self.z_rho_tot = depth.sdepth(Htot, zeta, self.hc, self.Cs_r,
-                                              Vtransform=self.Vtransform)
+                                              Vtransform=self.Vtransform,
+					      Vstretching=self.Vstretching)
                 # z_rho is positive relative to mean sea level but z is
                 # 0 at the surface.
                 # Transform z_rho to match convention of z.
@@ -517,7 +525,8 @@ def get_variables(self, requested_variables, time=None,
             H = self.sea_floor_depth_below_sea_level[indy, indx]
             zeta = self.zeta[itxy]
             z_rho = depth.sdepth(H, zeta, self.hc, self.Cs_r,
-                                 Vtransform=self.Vtransform)
+                                 Vtransform=self.Vtransform,
+				 Vstretching=self.Vstretching)
 
             # z_rho is positive relative to mean sea level but z is
             # 0 at the surface.
@@ -622,13 +631,11 @@ def get_mask(mask_name, imask, masks_store):
                             logger.debug('Calculating sigma2z-coefficients for whole domain')
                             starttime = datetime.now()
                             dummyvar = np.ones((O, M, N))
-                            dummy, self.s2z_total = depth.multi_zslice(dummyvar, self.z_rho_tot, self.zlevels)
+                            dummy, self.s2z_A, self.s2z_C, self.s2z_I, self.s2z_kmax = depth.multi_zslice(dummyvar, self.z_rho_tot, self.zlevels)
                             # Store arrays/coefficients
-                            self.s2z_A = self.s2z_total[0].reshape(len(self.zlevels), M, N)
-                            self.s2z_C = self.s2z_total[1].reshape(len(self.zlevels), M, N)
-                            #self.s2z_I = self.s2z_total[2].reshape(M, N)
-                            self.s2z_kmax = self.s2z_total[3]
-                            del self.s2z_total  # Free memory
+                            self.s2z_A = self.s2z_A.reshape(len(self.zlevels), M, N)
+                            self.s2z_C = self.s2z_C.reshape(len(self.zlevels), M, N)
+                            #self.s2z_I = self.s2z_I.reshape(M, N)
                             logger.info('Time: ' + str(datetime.now() - starttime))
                         if 'A' not in locals():
                             logger.debug('Re-using sigma2z-coefficients')
@@ -661,9 +668,8 @@ def get_mask(mask_name, imask, masks_store):
                             kmax = len(zle)  # Must be checked. Or number of sigma-layers?
                     if 'A' not in locals():
                         logger.debug('Calculating new sigma2z-coefficients')
-                        variables[par], s2z = depth.multi_zslice(
+                        variables[par], A,C,I,kmax = depth.multi_zslice(
                             variables[par], z_rho, variables['z'])
-                        A,C,I,kmax = s2z
                         # Reshaping to compare with subset of full array
                         #zle = np.arange(zi1, zi2)
                         #A = A.reshape(len(zle), len(indx), len(indy))