Merge pull request #1273 from knutfrode/dev

Dev

examples/example_long_cmems.py

@@ -11,9 +11,9 @@
 
     machine copernicusmarine login <your username> password <your password>
 
-This file must be stored in your home folder or in the main OpenDrift folder
+This file must be stored in your home folder (and unreadable by others) or in the main OpenDrift folder
 
-Alternatively, an Xarray dataset can be created explicitly with the copoernicusmarine client, and provided to reader_netCDF_CF_generic:
+Alternatively, an Xarray dataset can be created explicitly with the copernicusmarine client, and provided to reader_netCDF_CF_generic:
 https://opendrift.github.io/gallery/example_long_cmems_new.html
 """
 

examples/example_thredds_resources.py

@@ -16,7 +16,7 @@
 #%%
 # Open each thredds dataset to check contents and spatial coverage
 for t in thredds_resources:
-    if t.startswith('http') and 'nrt.cmems' not in t:
+    if t.startswith('http') and not t.startswith('cmems'):
         start = datetime.now()
         print('\n#%%\n%s\n' % t)
         r = Reader(t)

examples/example_vertical_mixing.py

@@ -41,6 +41,7 @@
 o.environment = np.array(list(zip(np.ones(N)*sea_floor_depth, np.zeros(N))),
                 dtype=[('sea_floor_depth_below_sea_level', np.float32),
                        ('sea_surface_height', np.float32)]).view(np.recarray)
+o.environment.ocean_mixed_layer_thickness = np.ones(N)*50
 o.environment_profiles = {
         'z': z,
         'ocean_vertical_diffusivity':

opendrift/models/oceandrift.py

@@ -138,6 +138,9 @@ def __init__(self, *args, **kwargs):
             'vertical_mixing:TSprofiles': {'type': 'bool', 'default': False, 'level':
                 CONFIG_LEVEL_ADVANCED,
                 'description': 'Update T and S profiles within inner loop of vertical mixing. This takes more time, but may be slightly more accurate.'},
+            'drift:profiles_depth': {'type': 'float', 'default': 50, 'min': 0, 'max': None,
+                'level': CONFIG_LEVEL_ADVANCED, 'units': 'meters', 'description':
+                'Environment profiles will be retrieved from surface and down to this depth'},
             'drift:wind_drift_depth': {'type': 'float', 'default': 0.1,
                 'min': 0, 'max': 10, 'units': 'meters',
                 'description': 'The direct wind drift (windage) is linearly decreasing from the surface value (wind_drift_factor) until 0 at this depth.',
@@ -445,8 +448,7 @@ def surface_wave_mixing(self, time_step_seconds):
 
         pass
 
-    def get_diffusivity_profile(self, model):
-        depths = np.abs(self.environment_profiles['z'])
+    def get_diffusivity_profile(self, model, depths):
         wind, depth = np.meshgrid(self.wind_speed(), depths)
         MLD = self.environment.ocean_mixed_layer_thickness
         background_diffusivity = self.get_config('vertical_mixing:background_diffusivity')
@@ -508,6 +510,8 @@ def vertical_mixing(self, store_depths=False):
         # get vertical eddy diffusivity from environment or specific model
         diffusivity_model = self.get_config('vertical_mixing:diffusivitymodel')
         diffusivity_fallback = self.get_config('environment:fallback:ocean_vertical_diffusivity')
+        # Mixing levels to be used if analytical expression:
+        mixing_z_analytical = -np.arange(0, self.environment.ocean_mixed_layer_thickness.max() + 2)
         if diffusivity_model == 'environment':
             if 'ocean_vertical_diffusivity' in self.environment_profiles and not (
                     self.environment_profiles['ocean_vertical_diffusivity'].min() ==
@@ -516,26 +520,22 @@ def vertical_mixing(self, store_depths=False):
                     diffusivity_fallback):
                 Kprofiles = self.environment_profiles[
                     'ocean_vertical_diffusivity']
+                mixing_z = self.environment_profiles['z']
                 logger.debug('Using diffusivity from ocean model')
             else:
                 logger.debug('Using diffusivity from Large1994 since model diffusivities not available')
-                # Using higher vertical resolution when analytical
-                self.environment_profiles['z'] = \
-                        -np.arange(0, self.environment.ocean_mixed_layer_thickness.max() + 2)
-                Kprofiles = self.get_diffusivity_profile('windspeed_Large1994')
+                mixing_z = mixing_z_analytical
+                Kprofiles = self.get_diffusivity_profile('windspeed_Large1994', np.abs(mixing_z))
         elif diffusivity_model == 'constant':
             logger.debug('Using constant diffusivity specified by fallback_values[''ocean_vertical_diffusivity''] = %s m2.s-1' % (diffusivity_fallback))
+            mixing_z = self.environment_profiles['z']
             Kprofiles = diffusivity_fallback*np.ones(
                     self.environment_profiles['ocean_vertical_diffusivity'].shape) # keep constant value for ocean_vertical_diffusivity
         else:
             logger.debug('Using functional expression for diffusivity')
-            # Using higher vertical resolution when analytical
-            if self.environment_profiles is None:
-                self.environment_profiles = {}
-            self.environment_profiles['z'] = \
-                    -np.arange(0, self.environment.ocean_mixed_layer_thickness.max() + 2)
             # Note: although analytical functions, z is discretised
-            Kprofiles = self.get_diffusivity_profile(diffusivity_model)
+            mixing_z = mixing_z_analytical
+            Kprofiles = self.get_diffusivity_profile(diffusivity_model, np.abs(mixing_z))
 
         logger.debug('Diffusivities are in range %s to %s' %
                       (Kprofiles.min(), Kprofiles.max()))
@@ -561,12 +561,11 @@ def vertical_mixing(self, store_depths=False):
             Tprofiles = None
 
         # prepare vertical interpolation coordinates
-        z_i = range(self.environment_profiles['z'].shape[0])
+        z_i = range(mixing_z.shape[0])
         if len(z_i) == 1:
             z_index = 0
         else:
-            z_index = interp1d(-self.environment_profiles['z'],
-                               z_i, bounds_error=False,
+            z_index = interp1d(-mixing_z, z_i, bounds_error=False,
                                fill_value=(0,len(z_i)-1))  # Extrapolation
 
         # Internal loop for fast time step of vertical mixing model.
@@ -583,7 +582,7 @@ def vertical_mixing(self, store_depths=False):
             depths[0, :] = self.elements.z
 
         # Calculating dK/dz for all profiles before the loop
-        gradK = -np.gradient(Kprofiles, self.environment_profiles['z'], axis=0)
+        gradK = -np.gradient(Kprofiles, mixing_z, axis=0)
         gradK[np.abs(gradK)<1e-10] = 0
 
         for i in range(0, ntimes_mix):

tests/models/test_readers.py

@@ -381,10 +381,43 @@ def test_vertical_profiles(self):
                                        time=norkyst3d.start_time,
                                        x=x, y=y, z=[0, -100])
         self.assertEqual(data['z'][4], -25)
-        self.assertEqual(data['z'][4], -25)
         self.assertAlmostEqual(data['sea_water_temperature'][:,0,0][7],
                          9.220000267028809)
 
+        # test also profiles depth
+        data, profiles = norkyst3d.get_variables_interpolated(
+                variables, profiles=['sea_water_temperature'],
+                profiles_depth = [-100, 0],
+                time = norkyst3d.start_time,
+                lon=lon, lat=lat, z=0)
+        assert profiles['z'].min() == -150
+        # Request profiles with less depth
+        # NB TODO: must use later time, otherwise cache (with deeper profile) is reused
+        data, profiles = norkyst3d.get_variables_interpolated(
+                variables, profiles=['sea_water_temperature'],
+                profiles_depth = [-30, 0],
+                time = norkyst3d.start_time + timedelta(hours=1),
+                lon=lon, lat=lat, z=0)
+        assert profiles['z'].min() == -75
+        # Setting vertical buffer to 0, and checking that now less extra layers are read in the vertical
+        norkyst3d.verticalbuffer=0
+        data, profiles = norkyst3d.get_variables_interpolated(
+                variables, profiles=['sea_water_temperature'],
+                profiles_depth = [-30, 0],
+                time = norkyst3d.start_time + timedelta(hours=2),
+                lon=lon, lat=lat, z=0)
+        assert profiles['z'].min() == -50
+
+    def test_vertical_profile_from_simulation(self):
+        norkyst3d = reader_netCDF_CF_generic.Reader(o.test_data_folder() +
+            '14Jan2016_NorKyst_z_3d/NorKyst-800m_ZDEPTHS_his_00_3Dsubset.nc')
+        lon = np.array([4.73])
+        lat = np.array([62.35])
+        variables = ['x_sea_water_velocity', 'x_sea_water_velocity',
+                     'sea_water_temperature']
+        # TODO TO BE COMPLETED
+
+
     def test_vertical_interpolation(self):
         norkyst3d = reader_netCDF_CF_generic.Reader(o.test_data_folder() +
             '14Jan2016_NorKyst_z_3d/NorKyst-800m_ZDEPTHS_his_00_3Dsubset.nc')
@@ -406,10 +439,6 @@ def test_vertical_interpolation(self):
         # Check interpolated temperature at 33 m depth
         self.assertAlmostEqual(data['sea_water_temperature'][1],
                                8.32, 2)
-        #import matplotlib.pyplot as plt
-        #plt.plot(profiles['sea_water_temperature'][:,0])
-        #plt.plot(profiles['sea_water_temperature'][:,1], 'r')
-        #plt.show()
 
     def test_vertical_interpolation_sigma(self):
         nordic3d = reader_ROMS_native.Reader(o.test_data_folder() +

tests/models/test_run.py

@@ -393,7 +393,8 @@ def test_vertical_mixing_profiles(self):
             o.environment = np.array([(100, 0) for _ in range(N)],
                          dtype=[('sea_floor_depth_below_sea_level', np.float32),
                                 ('sea_surface_height', np.float32)]).view(np.recarray)
-            o.environment_profiles = { 'z': z, 'ocean_vertical_diffusivity':
+            o.environment.ocean_mixed_layer_thickness = np.ones(N)*50
+            o.environment_profiles = {'z': z, 'ocean_vertical_diffusivity':
                                       np.tile(diffusivity, (N, 1)).T}
             o.env.finalize()
             o.vertical_mixing()