Passing multiple files to FieldSet in the case of vertically adaptive meshes

[tiagoams]

Parcels version

3.0.4

Description

I have managed to run Parcels with NEMO AMM7 output](#1339), using the original time varying sigma levels, when a single netcdf file per variable is passed. But when I pass a list of daily files to FieldSet.from_nemo() I get the error below.
Since passing a list of files works well for FieldSet.from_nemo() for non-time varying depths, and time varying depths works well when single files are used, I wonder if this something that could easily be modified in the code with some guidance. At the moment this is forcing me to concatenate grid_W netcdf files over the time dimension, hurting performance.

Small input files to reproduce the error are available here:
https://we.tl/t-mi0vuJogZe

Code sample

from parcels import FieldSet, ParticleSet, JITParticle, ScipyParticle, AdvectionRK4, AdvectionRK4_3D,  ParcelsRandom

from glob import glob
import numpy as np
from datetime import timedelta as delta
from datetime import datetime
from os import path
import xarray as xr
import math

data_dir = '/data/smallexample/'
datestr = '*1996010*.nc'     # multiple files
# datestr = '*19960101.nc'       # single files


def prep_fieldset(datestr, mesh_mask):
    
    ufiles = sorted(glob(data_dir + '*grid_U_' + datestr))
    vfiles = sorted(glob(data_dir + '*grid_V_' + datestr))
    wfiles = sorted(glob(data_dir + '*grid_W_' + datestr))

    filenames = {'U':        {'lon': mesh_mask, 'lat': mesh_mask, 'depth': ufiles, 'data': ufiles},
                 'V':        {'lon': mesh_mask, 'lat': mesh_mask, 'depth': vfiles, 'data': vfiles},
                 'W':        {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                 'Kh_v3d':   {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                 'dKh_v3d':  {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                 'depthw4d': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                 }

    variables = {'U': 'vozocrtx',
                 'V': 'vomecrty',
                 'W': 'vomecrtz',
                 'Kh_v3d': 'avm',
                 'dKh_v3d': 'avmgrad',
                 'depthw4d': 'gdepw_n',
                }

    dimensions = {'U':        {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'V':        {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'W':        {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'Kh_v3d':   {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'dKh_v3d':  {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'depthw4d': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  }

    fieldset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation=False)
    fieldset.U.set_depth_from_field(fieldset.depthw4d)
    fieldset.V.set_depth_from_field(fieldset.depthw4d)
    fieldset.W.set_depth_from_field(fieldset.depthw4d)
    fieldset.Kh_v3d.set_depth_from_field(fieldset.depthw4d)
    fieldset.dKh_v3d.set_depth_from_field(fieldset.depthw4d)
    fieldset.Kh_v3d.interp_method = 'cgrid_velocity'
    
    return fieldset

fieldset = prep_fieldset(datestr, mesh_path)

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
/tmp/ipykernel_57228/3516857616.py in <module>
     53     return fieldset
     54 
---> 55 fieldset = prep_fieldset(datestr, mesh_path)

/tmp/ipykernel_57228/3516857616.py in prep_fieldset(datestr, mesh_mask)
     41                   }
     42 
---> 43     fieldset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation=False)
     44     fieldset.U.set_depth_from_field(fieldset.depthw4d)
     45     fieldset.V.set_depth_from_field(fieldset.depthw4d)

~/.conda/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldset.py in from_nemo(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, chunksize, **kwargs)
    545         if kwargs.pop('gridindexingtype', 'nemo') != 'nemo':
    546             raise ValueError("gridindexingtype must be 'nemo' in FieldSet.from_nemo(). Use FieldSet.from_c_grid_dataset otherwise")
--> 547         fieldset = cls.from_c_grid_dataset(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
    548                                            allow_time_extrapolation=allow_time_extrapolation, tracer_interp_method=tracer_interp_method,
    549                                            chunksize=chunksize, gridindexingtype='nemo', **kwargs)

~/.conda/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldset.py in from_c_grid_dataset(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, gridindexingtype, chunksize, **kwargs)
    673             kwargs['creation_log'] = 'from_c_grid_dataset'
    674 
--> 675         return cls.from_netcdf(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
    676                                allow_time_extrapolation=allow_time_extrapolation, interp_method=interp_method,
    677                                chunksize=chunksize, gridindexingtype=gridindexingtype, **kwargs)

~/.conda/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldset.py in from_netcdf(cls, filenames, variables, dimensions, indices, fieldtype, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, chunksize, **kwargs)
    454                             dFiles = fields[procvar].dataFiles
    455                             break
--> 456             fields[var] = Field.from_netcdf(paths, (var, name), dims, inds, grid=grid, mesh=mesh, timestamps=timestamps,
    457                                             allow_time_extrapolation=allow_time_extrapolation,
    458                                             time_periodic=time_periodic, deferred_load=deferred_load,

~/.conda/envs/py3_parcels/lib/python3.9/site-packages/parcels/field.py in from_netcdf(cls, filenames, variable, dimensions, indices, grid, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, **kwargs)
    405             depth_filename = cls.get_dim_filenames(filenames, 'depth')
    406             if isinstance(filenames, dict) and len(depth_filename) != 1:
--> 407                 raise NotImplementedError('Vertically adaptive meshes not implemented for from_netcdf()')
    408             depth_filename = depth_filename[0]
    409 

NotImplementedError: Vertically adaptive meshes not implemented for from_netcdf()

[erikvansebille]

Hi @tiagoams, thanks for raising. To be honest, I'm not sure how easy/challenging it would be to implement this. You can give it a try (start with removing that raise NotImplementedError()), but what you will uncover I'm not sure

Note that we are currently working on a completely different approach to Sigma-grids, at #1641. While this is targeted for CROCO velocities, the approach/algorithm might also work for AMM7 output? Keen to hear your thoughts about that!

[tiagoams]

Hi @erikvansebille , many thanks for your help. I removed the check and my simple example with vertical diffusion kernel is working as expected.
I didn't manage to understand the CROCO compliant kernel, so I will wait for more documentation to emerge.