From ae071e4f74b97f39b502423de71bb182e116be2c Mon Sep 17 00:00:00 2001 From: Knut-Frode Dagestad Date: Fri, 12 Apr 2024 14:48:23 +0200 Subject: [PATCH] Levels of diffusivity for vertical mixing can be independent of z from environment_profiles, if analytical method is used --- examples/example_vertical_mixing.py | 1 + opendrift/models/oceandrift.py | 31 ++++++++++++++--------------- tests/models/test_readers.py | 10 ++++++++++ tests/models/test_run.py | 3 ++- 4 files changed, 28 insertions(+), 17 deletions(-) diff --git a/examples/example_vertical_mixing.py b/examples/example_vertical_mixing.py index 5dc773427..0363c98fb 100755 --- a/examples/example_vertical_mixing.py +++ b/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': diff --git a/opendrift/models/oceandrift.py b/opendrift/models/oceandrift.py index 79ffe980e..9850c3cad 100644 --- a/opendrift/models/oceandrift.py +++ b/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): diff --git a/tests/models/test_readers.py b/tests/models/test_readers.py index 913cbca69..d6c42d667 100644 --- a/tests/models/test_readers.py +++ b/tests/models/test_readers.py @@ -408,6 +408,16 @@ def test_vertical_profiles(self): 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') diff --git a/tests/models/test_run.py b/tests/models/test_run.py index 4638ba1f0..fc1ad9e76 100644 --- a/tests/models/test_run.py +++ b/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()