Caching also ocean depth, for calculation of z from sigma. Now profil…

…e test is passing.
OpenDrift · Nov 30, 2020 · bb40c78 · bb40c78
1 parent 307894b
commit bb40c78
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 7 deletions.
diff --git a/examples/example_fvcom.py b/examples/example_fvcom.py
@@ -20,8 +20,8 @@
 
 # Seed elements at defined positions, depth and time
 N = 1000
-z = -10*np.random.uniform(0, 1, N)
-o.seed_elements(lon=17.0, lat=70.0, radius=5000, number=N,
+z = -50*np.random.uniform(0, 1, N)
+o.seed_elements(lon=18.0, lat=69.8, radius=5000, number=N,
                 z=z, time=fvcom.start_time)
 
 #%%
@@ -43,3 +43,4 @@
 # Print and plot results
 print(o)
 o.plot(fast=True, buffer = 1.)
+# Two figure windows appear, this is probably a bug independent of FVCOM
diff --git a/opendrift/readers/reader_netCDF_CF_unstructured.py b/opendrift/readers/reader_netCDF_CF_unstructured.py
@@ -75,9 +75,11 @@ class Reader(BaseReader, UnstructuredReader):
 
     dataset = None
 
-    # For in-memory caching of Sigma-coordinates
+    # For in-memory caching of Sigma-coordinates and ocean depth
     siglay = None
     siglev = None
+    ocean_depth_nele = None
+    ocean_depth_node = None
 
     def __init__(self, filename=None, name=None, proj4=None):
         if filename is None:
@@ -343,12 +345,19 @@ def __nearest_node_sigma__(self, var, nodes, z):
             if self.siglay is None:
                 logger.debug('Reading and storing siglay array...')
                 self.siglay = self.dataset['siglay_center'][:]
-            depths = self.siglay[:, nodes]
+            sigmas = self.siglay[:, nodes]
         else:
             if self.siglev is None:
                 logger.debug('Reading and storing siglev array...')
                 self.siglev = self.dataset['siglev_center'][:]
-            depths = self.siglev[:, nodes]
+            sigmas = self.siglev[:, nodes]
+
+        if self.ocean_depth_node is None:
+            logger.debug('Reading and storing ocean depth...')
+            self.ocean_depth_node = self.dataset['h'][:]
+
+        # Calculating depths from sigmas
+        depths = sigmas*self.ocean_depth_node[nodes]
 
         return self._vector_nearest_(depths, z)
 
@@ -362,12 +371,19 @@ def __nearest_face_sigma__(self, var, el, z):
             if self.siglay is None:
                 logger.debug('Reading and storing siglay array...')
                 self.siglay = self.dataset['siglay_center'][:]
-            depths = self.siglay[:, el]
+            sigmas = self.siglay[:, el]
         else:
             if self.siglev is None:
                 logger.debug('Reading and storing siglev array...')
                 self.siglev = self.dataset['siglev_center'][:]
-            depths = self.siglev[:, el]
+            sigmas = self.siglev[:, el]
+
+        if self.ocean_depth_nele is None:
+            logger.debug('Reading and storing ocean depth...')
+            self.ocean_depth_nele = self.dataset['h_center'][:]
+
+        # Calculating depths from sigmas
+        depths = sigmas*self.ocean_depth_nele[el]
 
         return self._vector_nearest_(depths, z)