General:

• change name g_default(e.g. to grav). If we are to change this name, we should make sure the modification is also made in gridop and vmodes

In vmodes:

• method for saving and loading modes

@apatlpo , would you mind fixing this (that should be easy for you)?
I think that's a typo. When we specify chunks when calling  postp.Run, the corresponding dataset is chunkmethod... See this line, I guess:

via XIOS

interp2z and interp2z_3d should be updated in gridop.py to be compatible with the fast_interp routine. The latter now includes nearest or linear extrapolation at the top and the bottom, and hence requires 2 more arguments (integers -- 5 args in total). I think the corrected versions are the one in proj_modes.ipynb (in tutorials). Aurélien did the modification of the fast_interp routine, can you confirm and make the changes accordingly?

if v (the DataArrayto interpolate) has no name, this raises the error:

ValueError: cannot reset_coords with drop=False on an unnamed DataArrray

If I remember correctly, the oldest strategy was to load his by default, then switched to grid, and now we try to be compatible. Here are a few questions regarding these few lines:

• first line: shouldn't it be t1/{}grid.nc'.format(prefix)? Otherwise Run.__init__most certainly misses the grid file and switched to old behavior, i.e. load his(and only his) by default. I made the corresponding change in the branch "vmodes".
• Would that be compatible to use _nc_default only if open_nc is empty?

I think l. 321 should read "_f0, _beta = None, None" istead of "f0, beta = None, none"

Puisque f est présent dans le netCDF grid, pourquoi ne pas l'utiliser comme coordonnée ? Ca évite d'avoir à donner f0, beta s'ils ne sont pas dans les fichiers de sortie (ce ne sont que des attributs) et de faire potentiellement des erreurs. Par contre on perdrai f_u, f_v, mais ils se recalculent très facilement avec grid.interp

Si OK je ferai les modifs nécessaires et ferai un pull request

Hi I noticed something that I was not expecting: it seems like interp2z requires both vertical grids to be ordered the same way. It looks like this behaviour is inherited from fast_interp3D.interp routine. Would you mind check this?
https://nbviewer.jupyter.org/github/NoeLahaye/bacasable/blob/master/debug_interpz_order.ipynb
@apatlpo , @slgentil

What is strange is that in my examples, sometimes I get Nan (I would get wrong values if I use boundary extrapolation), while I never have any target value that is outside the range of initial values (for the z grid), so I hope there is not something worse behind this

All the best,

Noé

List of items to address with fsturb configuration:

• check forcing on density works for a purely barotropic forcing: it doesn't. Forcing on velocity was implemented instead

• investigate the forcing amplitude to energy response relationship

• figure out if relaxation needs to be turned back on

• implement simple movie generator and make pretty movies

• implement modal decomposition of the flow and modal energy temporal evolution

• implement atmospheric forcing for NIW generation

• turn internal tide back on and adjust sponge layers

we realized that today with Zoé, this time looking at iwave polarization ... amazing !

Can somebody take a look at this?

@apatlpo, @slgentil
I spotted a few bugs / shortcomings, regarding how xgrid (do not) adapt to changes in datasets, e.g. when selecting portions of a dataset. This is illustrated in this notebook

Projections.

API could be somewhat clarified I believe. When one wants to project onto modes there are several cases:

• 1D vertical modes: all goes well, in other words, the data horizontal grid does not matter

• horizontally varying modes: in that case, data needs to be interpolated horizontally first on the vertical mode grid, otherwise the projection will broadcast horizontal dimensions that differ. If this statement is correct we need to more explicitly state this in the docstring. Right now it is unclear whether this horizontal interpolation may be performed within the call to projector not. The notebook for example falsely lead to believe it is the case, as z_rho (on hgrid rho) is passed to project.

very minor detail:

• vmodes.ipynb / cell 28: set aspect to equal (subplot_kws={'aspect':'equal'})

We need to find a good way to pass options about chunking when opening CROCOrun objects.
Maybe something along the lines of :

r = pp.CROCOrun(root_path+run, prefix='file_', open_nc=['surf'], ...
 chunks={'his': {'time':1}, 'surf': {'time':None}})

any other ideas?

it probably should look like something closer to:

def get_z(run, zeta=None, h=None, vgrid='r', hgrid='r'):
    ''' compute vertical coordinates
        zeta should have the size of the final output
        vertical coordinate should be first
    '''

    ds = run['his']
    N = run.N
    hc = run.params['Hc']

    if zeta is not None:
        _zeta = zeta
    else:
        _zeta = ds.ssh
    if h is not None:
        _h = h
    else:
        _h = ds.h
    #
    if hgrid is 'u':
        _zeta = rho2u(_zeta, ds)
        _h = rho2u(_h, ds)
    elif hgrid is 'v':
        _zeta = rho2v(_zeta, ds)
        _h = rho2v(_h, ds)
    #
    sc=run['his']['sc_'+vgrid]
    cs=run['his']['Cs_'+vgrid]

    #
    z0 = (hc * sc + _h * cs) / (hc + _h)
    z = _zeta + (_zeta + _h) * z0
    #z = np.squeeze(_zeta + (_zeta + _h) * z0)  # this should not be here
    # manually swap dims, could also perform align with T,S
    #if z.ndim ==4:
    #    z = z.transpose(z.dims[0], z.dims[3], z.dims[1], z.dims[2])
    #elif z.ndim == 3:
    #    z = z.transpose(sc.dims[0], _zeta.dims[0], _zeta.dims[1])
    #elif z.ndim == 2:
    #    z = z.transpose(sc.dims[0], _zeta.dims[0])
    return z.rename('z_'+vgrid)

I am not sure transposition is necessary, to discuss

I added a version of get_p in gridop. What is implemented is a "finite volume" implementation, where the density is integrated level-wise (assuming rho is constant in a given sigma rho-layer) from the surface. i.e. it is strictly equivalent to cumsum on w-levels. This is probably not what is actually implemented in croco, but it gives something satisfying.

Please check that it looks fine (and that it works if you want, but I am currently using it so...) and delete the old version get_p_old if it is OK

List of issues that would be worth discussing in a meeting

• how to contribute~: git branch, merge, pull request, using issue, tracking modifications, making suggestions, asking questions (on github)
• strategies for postp and gridop
  ** what to do with grid (xgcm-Grid) ? where to attach it ?