OpenDrift / opendrift

Open source framework for ocean trajectory modelling
https://opendrift.github.io
GNU General Public License v2.0
247 stars 121 forks source link

SCHISM native reader- profile variable #707

Open mymyrere opened 3 years ago

mymyrere commented 3 years ago

Hi, When using OpenOil model and the schism reader I get a crash line 1072 of the basereader because self. env_profiles[var][np.ix_(z_ind, missing_indices)] = \ np.ma.masked_invalid(env_profiles_tmp[var][z_ind,0:len(missing_indices)]).astype('float32')

This is because the shape of env_profiles_tmp[var] is (1,len(particle)) and env_profiles_tmp[z] is the number of node in the schism file.

What about doing an interpolation for fix depth in the schism_reader:

    self.z=[0,-5,10,-15,-20]
    if profiles is not None and varname in profiles:
          profile_data=np.zeros((len(self.z),len(x)))

          for ii,zz in enumerate(self.z):
              dist,i=self.block_KDtree_3d.query((np.vstack((x,y,np.ones(z.shape)*zz)).T).astype(np.float32),nb_closest_nodes)
              dist[dist<DMIN]=DMIN
              fac=(1./dist)
              profile_data[ii,:] = (fac*data.take(i)).sum(-1)/fac.sum(-1)
         profiles_dict[varname] = profile_data # horizontal

I don't really understand how the profile works, so not sure what would happen if all particle are below -20m? Does it need a profile from the minimum to the maximum depth of the particle?

Thank you very much

knutfrode commented 3 years ago

Hi,

I am tagging here @simonweppe , who is the author of the SCHISM reader, and who might be able to help. I am not into all the details of this reader, and how it deals with environment profiles.

Anyway, the point of the environment profiles, is that when calculating vertical turbulent mixing and buoyancy, you need to know not only the environment variables (temp, salt, ...) at the actual position/depth of the particle, but also the vertical density gradient and diffusivity, which is calculated. Thus, in addition to returning a 1D array of environment variables (env) at the positions and depths of the particles, readers should also return a 2D array (env_profiles) with dimension (particles, depths).

I am not sure how much effort is required to make SCHISM reader provide profiles (perhaps only a quick fix?), but in the meantime, it could perhaps work to disable the vertical motion of oil by o.disable_vertical_motion()

mymyrere commented 3 years ago

Thank you, I will work with Simon on this

simonweppe commented 1 month ago

@knutfrode - coming back to this...late !

Is there an example somewhere where vertical diffusion values are interpolated to vertical profiles (e.g. from 3D diffusion fields) to then be used in the vertical mixing i.e. as per : o.set_config('vertical_mixing:diffusivitymodel','environment')

Will be useful to understand the process and adapt to unstructured reader block

Thanks

simonweppe commented 1 month ago

@calquigs just adding you here so we can keep track

simonweppe commented 1 month ago

@knutfrode sorry to bug you with this. Is there some an example , or sample file, that could be used to test the 3D interpolation of diffusion field on regular grid , so I can them mirror that in the SCHISM reader. Thanks !

knutfrode commented 1 month ago

Hi, I apologize for the very late reply. I wanted to give a proper answer, and needed first to dig into this, which I did not work with for a long time. So far I did not find time for a proper investigation, and am this week on a workshop.

But some quick comments in the meantime: There are no test-datasets within OpenDrift repository which contain 3D-values of vertical diffusivity, but here is one example on thredds: https://thredds.met.no/thredds/dodsC/fou-hi/norkystv3_his_files/2024/09/25/norkyst800_his_zdepth_20240925T00Z_m00_AN.nc

Below is a small script performing vertical mixing with some particles based on the vertical diffusivities from this file. NOTE: CF standard_name is here ocean_vertical_salt_diffusivity and I just now added a mapping from this to ocean_vertical_diffusivity. Thus you must update your OpenDrift to run the example below.

import numpy as np
from datetime import datetime, timedelta
from opendrift.models.oceandrift import OceanDrift
from opendrift.readers.reader_netCDF_CF_generic import Reader

o = OceanDrift(loglevel=0)  # loglevel 0 to see all details about what is hppening
o.set_config('drift:vertical_mixing', True)

r = Reader('https://thredds.met.no/thredds/dodsC/fou-hi/norkystv3_his_files/2024/09/25/norkyst800_his_zdepth_20240925T00Z_m00_AN.nc')
o.add_reader(r)

# Seed one element each m below surface down to 30m
z = np.arange(-1, -30, -1)
o.seed_elements(lon=4, lat=60, time=r.start_time, z=z, number=len(z))

o.run(time_step=timedelta(minutes=10), duration=timedelta(hours=3))
o.plot_property('z', filename='vertical_diffusivity_depth.png')
o.plot(linecolor='ocean_vertical_diffusivity', buffer=0, filename='vertical_diffusivity_horizontal.png')

The output figures are attached.

I am not sure if this example is useful, but you can have a look and provide some feedback - and I can add more information later. vertical_diffusivity_horizontal vertical_diffusivity_depth

simonweppe commented 1 month ago

Thank you very much @knutfrode - I think that will be really useful. I'll keep you posted for a potential update of the schism reader.

knutfrode commented 1 month ago

Good. I can add one thing: Presently OpenDrift internally always uses z as the vertical coordinate, and this e.g. the ROMS native reader must regrid the blocks of data from (x,y,sigma) to (x,y,z). This can be quite slow for very large datasets/files, and also involves hardcoding some z-levels which might not always be appropriate. Thus I am considering having a generic internal vertical coordinate in OpenDrift so that this regridding is not anymore necessary. But this is not a trivial change, so will require some careful planning first. E.g. the vertical mixing algorithm presently calculates dK/dz, but it must then also be possible to use dK/dsigma

I am not sure if unstructured files (e.g. Schism) always have sigma as vertical coordinate, or if they can also have z as vertical coordinate?

simonweppe commented 1 month ago

Yes good point. For SCHISM, there can be z levels, or sigma, or a mix of both...and they typically change over time following sea surface elevation (might be the same for ROMS). There is a coordinate zcor that stores the vertical levels at each step (that my reader requires for any 3D interp) so information is there and could be used to work out the dz or dsigma . I'll see how I can handle that when I properly dive into this.

knutfrode commented 1 month ago

Ok, so this means that z (zcor) is a variable (i.e. 4D variable) and not a coordinate (1D variable)?

That can then probably be used, but presently env_profiles requires all profiles to have the same z-levels - thus a regridding is (unfortunately) necessary. But could be an idea to make z a variable also internally in OpenDrift with regards to env_profiles, instead of a coordinate as presently.

Ideally one would like to have a generic joint reader for all unstructured models, but I guess there are still a bit of differences in the actual output files.

simonweppe commented 1 month ago

so this means that z (zcor) is a variable (i.e. 4D variable) >> Yes that's right.

Ok noted about the regridding, I agree ideally we'd pull out the full profile with no re-interpolate if it is already available in the source files. I'll see if I can come up with something that could be more generic.