source: TOOLS/WATER_BUDGET/nemo.py @ 6652

# -*- coding: utf-8 -*-
## ===========================================================================
##
##  This software is governed by the CeCILL  license under French law and
##  abiding by the rules of distribution of free software.  You can  use,
##  modify and/ or redistribute the software under the terms of the CeCILL
##  license as circulated by CEA, CNRS and INRIA at the following URL
##  "http://www.cecill.info".
##
##  Warning, to install, configure, run, use any of Olivier Marti's
##  software or to read the associated documentation you'll need at least
##  one (1) brain in a reasonably working order. Lack of this implement
##  will void any warranties (either express or implied).
##  O. Marti assumes no responsability for errors, omissions,
##  data loss, or any other consequences caused directly or indirectly by
##  the usage of his software by incorrectly or partially configured
##  personal.
##
## ===========================================================================
'''
Utilities to plot NEMO ORCA fields
Periodicity and other stuff

- Lots of tests for xarray object
- Not much testerd for numpy objects

olivier.marti@lsce.ipsl.fr

## SVN information
__Author__   = "$Author$"
__Date__     = "$Date$"
__Revision__ = "$Revision$"
__Id__       = "$Id: $"
__HeadURL    = "$HeadURL$"
'''

import numpy as np
try    : import xarray as xr
except ImportError : pass

#try    : import f90nml
#except : pass

#try : from sklearn.impute import SimpleImputer
#except : pass

rpi = np.pi ; rad = np.deg2rad (1.0) ; dar = np.rad2deg (1.0)

nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2]

rday   = 24.*60.*60.     # Day length [s]
rsiyea = 365.25 * rday * 2. * rpi / 6.283076 # Sideral year length [s]
rsiday = rday / (1. + rday / rsiyea)
raamo  =  12.        # Number of months in one year
rjjhh  =  24.        # Number of hours in one day
rhhmm  =  60.        # Number of minutes in one hour
rmmss  =  60.        # Number of seconds in one minute
omega  = 2. * rpi / rsiday # Earth rotation parameter [s-1]
ra     = 6371229.    # Earth radius [m]
grav   = 9.80665     # Gravity [m/s2]
repsi  = np.finfo (1.0).eps

## Default names of dimensions
dim_names = {'x':'xx', 'y':'yy', 'z':'olevel', 't':None}

## All possibles name of dimensions in Nemo files
xName = [ 'x', 'X', 'X1', 'xx', 'XX', 'x_grid_T', 'x_grid_U', 'x_grid_V', 'x_grid_F', 'x_grid_W', 'lon', 'nav_lon', 'longitude', 'X1', 'x_c', 'x_f', ]
yName = [ 'y', 'Y', 'Y1', 'yy', 'YY', 'y_grid_T', 'y_grid_U', 'y_grid_V', 'y_grid_F', 'y_grid_W', 'lat', 'nav_lat', 'latitude' , 'Y1', 'y_c', 'y_f', ]
zName = [ 'z', 'Z', 'Z1', 'zz', 'ZZ', 'depth', 'tdepth', 'udepth', 'vdepth', 'wdepth', 'fdepth', 'deptht', 'depthu', 'depthv', 'depthw', 'depthf', 'olevel', 'z_c', 'z_f', ]
tName = [ 't', 'T', 'tt', 'TT', 'time', 'time_counter', 'time_centered', ]

## All possibles name of units of dimensions in Nemo files
xUnit = [ 'degrees_east', ]
yUnit = [ 'degrees_north', ]
zUnit = [ 'm', 'meter', ]
tUnit = [ 'second', 'minute', 'hour', 'day', 'month', 'year', ]

## All possibles size of dimensions in Orca files
xLength = [ 180, 182, 360, 362 ]
yLength = [ 148, 149, 331, 332 ]
zLength = [31, 75]

## ===========================================================================
def __mmath__ (tab, default=None) :
    '''
    Determines the type of tab : xarray or numpy object ?
    '''
    mmath = default
    try    :
        if type (tab) == xr.core.dataarray.DataArray : mmath = xr
    except : pass

    try    :
        if type (tab) == np.ndarray : mmath = np
    except : pass
            
    return mmath

def __guessNperio__ (jpj, jpi, nperio=None, out='nperio') :
    '''
    Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details)
    
    Inputs
    jpj    : number of latitudes
    jpi    : number of longitudes
    nperio : periodicity parameter
    '''
    if nperio == None :
        nperio = __guessConfig__ (jpj, jpi, nperio=None, out='nperio')
    
    return nperio

def __guessConfig__ (jpj, jpi, nperio=None, config=None, out='nperio') :
    '''
    Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details)

    Inputs
    jpj    : number of latitudes
    jpi    : number of longitudes
    nperio : periodicity parameter
    '''
    print ( jpi, jpj)
    if nperio == None :
        ## Values for NEMO version < 4.2
        if (jpj ==  149 and jpi == 182) or (jpj == None and jpi == 182) or (jpj == 149 or jpi == None) :
            config = 'ORCA2.3'
            nperio = 4   # ORCA2. We choose legacy orca2.
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'T'
        if (jpj == 332 and jpi == 362) or (jpj == None and jpi == 362) or (jpj ==  332 and jpi == None) : # eORCA1.
            config = 'eORCA1.2'
            nperio = 6  
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
        if jpi == 1442 :  # ORCA025.
            config = 'ORCA025'
            nperio = 6 
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
        if jpj ==  294 : # ORCA1
            config = 'ORCA1'
            nperio = 6
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
            
        ## Values for NEMO version >= 4.2. No more halo points
        if (jpj == 148 and jpi == 180) or (jpj == None and jpi == 180) or (jpj == 148 and jpi == None) :
            config = 'ORCA2.4'
            nperio = 4.2 # ORCA2. We choose legacy orca2.
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
        if (jpj == 331 and jpi == 360) or (jpj == None and jpi == 360) or (jpj == 331 and jpi == None) : # eORCA1.
            config = 'eORCA1.4'
            nperio = 6.2
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
        if jpi == 1440 : # ORCA025.
            config = 'ORCA025'
            nperio = 6.2
            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
            
        if nperio == None :
            raise Exception  ('in nemo module : nperio not found, and cannot by guessed')
        else :
            if nperio in nperio_valid_range :
                print ( f'nperio set as {nperio} (deduced from {jpj=} and {jpi=})' )
            else : 
                raise ValueError ( f'nperio set as {nperio} (deduced from {jpi=} and {jpj=}) : nemo.py is not ready for this value' )

    if out == 'nperio' : return nperio
    if out == 'config' : return config
    if out == 'perio'  : return Iperio, Jperio, NFold, NFtype
    if out in ['full', 'all'] : return {'nperio':nperio, 'Iperio':Iperio, 'Jperio':Jperio, 'NFold':NFold, 'NFtype':NFtype}
        
def __guessPoint__ (ptab) :
    '''
    Tries to guess the grid point (periodicity parameter. See NEMO documentation for details)
    
    For array conforments with xgcm requirements

    Inputs
         ptab : xarray array

    Credits : who is the original author ?
    '''
    
    gP = None
    mmath = __mmath__ (ptab)
    if mmath == xr :
        if 'x_c' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'T'
        if 'x_f' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'U'
        if 'x_c' in ptab.dims and 'y_f' in ptab.dims                        : gP = 'V'
        if 'x_f' in ptab.dims and 'y_f' in ptab.dims                        : gP = 'F'
        if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_c' in ptab.dims : gP = 'T'
        if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'W'
        if 'x_f' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'U'
        if 'x_c' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'V'
        if 'x_f' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'F'
             
        if gP == None :
            raise Exception ('in nemo module : cd_type not found, and cannot by guessed')
        else :
            print ( f'Grid set as {gP} deduced from dims {ptab.dims}' )
            return gP
    else :
         raise Exception  ('in nemo module : cd_type not found, input is not an xarray data')

def get_shape ( ptab ) :
    '''
    Get shape of ptab : 
    shape main contain X, Y, Z or T
    Y is missing for a latitudinal slice
    X is missing for on longitudinal slice
    etc ...
    '''
    
    get_shape = ''
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')
    kz, az = __findAxis__ (ptab, 'z')
    lt, at = __findAxis__ (ptab, 't')
    if ax : get_shape = 'X'
    if ay : get_shape = 'Y' + get_shape
    if az : get_shape = 'Z' + get_shape
    if at : get_shape = 'T' + get_shape
    return get_shape
     
def lbc_diag (nperio) :
    lperio = nperio ; aperio = False
    if nperio == 4.2 :
        lperio = 4 ; aperio = True
    if nperio == 6.2 :
        lperio = 6 ; aperio = True
        
    return lperio, aperio

def __findAxis__ (tab, axis='z') :
    '''
    Find order and name of the requested axis
    '''
    mmath = __mmath__ (tab)
    ix = None ; ax = None

    if axis in xName : axName = xName ; unList = xUnit ; Length = xLength
    if axis in yName : axName = yName ; unList = yUnit ; Length = yLength
    if axis in zName : axName = zName ; unList = zUnit ; Length = zLength
    if axis in tName : axName = tName ; unList = tUnit ; Length = None
    
    if mmath == xr :
        for Name in axName :
            try    :
                ix = tab.dims.index (Name)
                ax = Name
            except : pass

        for i, dim in enumerate (tab.dims) :
            if 'units' in tab.coords[dim].attrs.keys() :
                for name in unList :
                    if name in tab.coords[dim].attrs['units'] :
                        ix = i ; ax = dim
    else :
        #if axis in xName : ix=-1
        #if axis in yName :
        #    if len(tab.shape) >= 2 : ix=-2
        #if axis in zName :
        #    if len(tab.shape) >= 3 : ix=-3
        #if axis in tName :
        #    if len(tab.shape) >=3  : ix=-3
        #    if len(tab.shape) >=4  : ix=-4

        l_shape = tab.shape
        for nn in np.arange ( len(l_shape)) :
            if l_shape[nn] in Length : ix = nn
       
    return ix, ax

def findAxis ( tab, axis= 'z' ) :
  ix, xx = __findAxis__ (tab, axis)
  return xx

def fixed_lon (lon, center_lon=0.0) :
    '''
    Returns corrected longitudes for nicer plots

    lon        : longitudes of the grid. At least 2D.
    center_lon : center longitude. Default=0.

    Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7
    '''
    mmath = __mmath__ (lon)
    
    fixed_lon = lon.copy ()
        
    fixed_lon = mmath.where (fixed_lon > center_lon+180., fixed_lon-360.0, fixed_lon)
    fixed_lon = mmath.where (fixed_lon < center_lon-180., fixed_lon+360.0, fixed_lon)
    
    for i, start in enumerate (np.argmax (np.abs (np.diff (fixed_lon, axis=-1)) > 180., axis=-1)) :
        fixed_lon [..., i, start+1:] += 360.

    # Special case for eORCA025
    if fixed_lon.shape [-1] == 1442 : fixed_lon [..., -2, :] = fixed_lon [..., -3, :]
    if fixed_lon.shape [-1] == 1440 : fixed_lon [..., -1, :] = fixed_lon [..., -2, :]

    if fixed_lon.min () > center_lon : fixed_lon += -360.0
    if fixed_lon.max () < center_lon : fixed_lon +=  360.0
        
    if fixed_lon.min () < center_lon-360.0 : fixed_lon +=  360.0
    if fixed_lon.max () > center_lon+360.0 : fixed_lon += -360.0
                
    return fixed_lon

def bounds_clolon ( bounds_lon, lon, rad=False, deg=True) :
    '''Choose closest to lon0 longitude, adding or substacting 360° if needed'''

    if rad : lon_range = 2.0*np.pi
    if deg : lon_range = 360.0
    bounds_clolon = bounds_lon.copy ()

    bounds_clolon = xr.where ( bounds_clolon < lon-lon_range/2., bounds_clolon+lon_range, bounds_clolon )
    bounds_clolon = xr.where ( bounds_clolon > lon+lon_range/2., bounds_clolon-lon_range, bounds_clolon )

    return bounds_clolon

def UnifyDims ( dd, udims=dim_names, verbose=False ) :
    '''
    Rename dimensions to unify them between NEMO versions
    '''
    
    if udims['x'] :
        for xx in xName :
            if xx in dd.dims and xx != udims['x'] :
                if verbose : print ( f"{xx} renamed to {udims['x']}" )
                dd = dd.rename ( {xx:udims['x']})
    if udims['y'] :
        for yy in yName :
            if yy in dd.dims and yy != udims['y']  :
                if verbose : print ( f"{yy} renamed to {udims['y']}" )
                dd = dd.rename ( {yy:udims['y']} )
    if udims['z'] :
        for zz in zName :
            if zz in dd.dims and zz != udims['z'] :
                if verbose : print ( f"{zz} renamed to {udims['z']}" )
                dd = dd.rename ( {zz:udims['z']} )
    if udims['t'] :
        for tt in tName  :
            if tt in dd.dims and tt != udims['t'] :
                if verbose : print ( f"{tt} renamed to {udims['t']}" )
                dd = dd.rename ( {tt:udims['t']} )

    return dd

def fill_empty (ztab, sval=np.nan, transpose=False) :
    '''
    Fill values

    Useful when NEMO has run with no wet points options : 
    some parts of the domain, with no ocean points, have no
    values
    '''
    from sklearn.impute import SimpleImputer
    mmath = __mmath__ (ztab)

    imp = SimpleImputer (missing_values=sval, strategy='mean')
    if transpose :
        imp.fit (ztab.T)
        ptab = imp.transform (ztab.T).T
    else : 
        imp.fit (ztab)
        ptab = imp.transform (ztab)
   
    if mmath == xr :
        ptab = xr.DataArray (ptab, dims=ztab.dims, coords=ztab.coords)
        ptab.attrs = ztab.attrs
        
    return ptab

def fill_lonlat (lon, lat, sval=-1) :
    '''
    Fill longitude/latitude values

    Useful when NEMO has run with no wet points options : 
    some parts of the domain, with no ocean points, have no
    lon/lat values
    '''
    from sklearn.impute import SimpleImputer
    mmath = __mmath__ (lon)

    imp = SimpleImputer (missing_values=sval, strategy='mean')
    imp.fit (lon)
    plon = imp.transform (lon)
    imp.fit (lat.T)
    plat = imp.transform (lat.T).T

    if mmath == xr :
        plon = xr.DataArray (plon, dims=lon.dims, coords=lon.coords)
        plat = xr.DataArray (plat, dims=lat.dims, coords=lat.coords)
        plon.attrs = lon.attrs ; plat.attrs = lat.attrs
        
    plon = fixed_lon (plon)
    
    return plon, plat

def fill_bounds_lonlat (bounds_lon, bounds_lat, sval=-1) :
    '''
    Fill longitude/latitude bounds values

    Useful when NEMO has run with no wet points options : 
    some parts of the domain, with no ocean points, as no
    lon/lat values
    '''
    mmath = __mmath__ (bounds_lon)

    p_bounds_lon = np.empty ( bounds_lon.shape )
    p_bounds_lat = np.empty ( bounds_lat.shape )

    imp = SimpleImputer (missing_values=sval, strategy='mean')
    
    for n in np.arange (4) : 
        imp.fit (bounds_lon[:,:,n])
        p_bounds_lon[:,:,n] = imp.transform (bounds_lon[:,:,n])
        imp.fit (bounds_lat[:,:,n].T)
        p_bounds_lat[:,:,n] = imp.transform (bounds_lat[:,:,n].T).T
        
    if mmath == xr :
        p_bounds_lon = xr.DataArray (bounds_lon, dims=bounds_lon.dims, coords=bounds_lon.coords)
        p_bounds_lat = xr.DataArray (bounds_lat, dims=bounds_lat.dims, coords=bounds_lat.coords)
        p_bounds_lon.attrs = bounds_lat.attrs ; p_bounds_lat.attrs = bounds_lat.attrs
        
    return p_bounds_lon, p_bounds_lat

def jeq (lat) :
    '''
    Returns j index of equator in the grid
    
    lat : latitudes of the grid. At least 2D.
    '''
    mmath = __mmath__ (lat)
    ix, ax = __findAxis__ (lat, 'x')
    jy, ay = __findAxis__ (lat, 'y')

    if mmath == xr :
        jeq = int ( np.mean ( np.argmin (np.abs (np.float64 (lat)), axis=jy) ) )
    else : 
        jeq = np.argmin (np.abs (np.float64 (lat[...,:, 0])))
    return jeq

def lon1D (lon, lat=None) :
    '''
    Returns 1D longitude for simple plots.
    
    lon : longitudes of the grid
    lat (optionnal) : latitudes of the grid
    '''
    mmath = __mmath__ (lon)
    jpj, jpi  = lon.shape [-2:]
    if np.max (lat) :
        je    = jeq (lat)
        #lon1D = lon.copy() [..., je, :]
        lon0 = lon [..., je, 0].copy()
        dlon = lon [..., je, 1].copy() - lon [..., je, 0].copy()
        lon1D = np.linspace ( start=lon0, stop=lon0+360.+2*dlon, num=jpi )
    else :
        lon0 = lon [..., jpj//3, 0].copy()
        dlon = lon [..., jpj//3, 1].copy() - lon [..., jpj//3, 0].copy()
        lon1D = np.linspace ( start=lon0, stop=lon0+360.+2*dlon, num=jpi )

    #start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1)
    #lon1D [..., start+1:] += 360

    if mmath == xr :
        lon1D = xr.DataArray( lon1D, dims=('lon',), coords=(lon1D,))
        lon1D.attrs = lon.attrs
        lon1D.attrs['units']         = 'degrees_east'
        lon1D.attrs['standard_name'] = 'longitude'
        lon1D.attrs['long_name :']   = 'Longitude'
        
    return lon1D

def latreg (lat, diff=0.1) :
    '''
    Returns maximum j index where gridlines are along latitudes in the northern hemisphere
    
    lat : latitudes of the grid (2D)
    diff [optional] : tolerance
    '''
    mmath = __mmath__ (lat)
    if diff == None :
        dy   = np.float64 (np.mean (np.abs (lat - np.roll (lat,shift=1,axis=-2, roll_coords=False))))
        print ( f'{dy=}' )
        diff = dy/100.

    je     = jeq (lat)
    jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< diff)[-1][-1] + je
    latreg = np.float64 (lat[...,jreg,:].mean(axis=-1))
    JREG   = jreg

    return jreg, latreg

def lat1D (lat) :
    '''
    Returns 1D latitudes for zonal means and simple plots.

    lat : latitudes of the grid (2D)
    '''
    mmath = __mmath__ (lat)
    jpj, jpi = lat.shape[-2:]

    dy     = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2))))
    je     = jeq (lat)
    lat_eq = np.float64 (lat[...,je,:].mean(axis=-1))
     
    jreg, lat_reg = latreg (lat)
    lat_ave = np.mean (lat, axis=-1)

    #print ( f'{dy=} {jpj=} {je=} {lat_eq=} {jreg=} ' )
    
    if (np.abs (lat_eq) < dy/100.) : # T, U or W grid
        if jpj-1 > jreg : dys = (90.-lat_reg) / (jpj-jreg-1)*0.5
        else            : dys = (90.-lat_reg) / 2.0
        yrange = (90.-dys-lat_reg)
    else                           :  # V or F grid
        yrange = 90.-lat_reg

    if jpj-1 > jreg :
        lat1D = mmath.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) )
    else :
        lat1D = lat_ave
    lat1D[-1] = 90.0

    if mmath == xr :
        lat1D = xr.DataArray( lat1D.values, dims=('lat',), coords=(lat1D,))
        lat1D.attrs = lat.attrs
        lat1D.attrs ['units']         = 'degrees_north'
        lat1D.attrs ['standard_name'] = 'latitude'
        lat1D.attrs ['long_name :']   = 'Latitude'
        
    return lat1D

def latlon1D (lat, lon) :
    '''
    Returns simple latitude and longitude (1D) for simple plots.

    lat, lon : latitudes and longitudes of the grid (2D)
    '''
    return lat1D (lat),  lon1D (lon, lat)

def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) :
    mmath = __mmath__ (ptab)
    try :
        lon = lon.copy().to_masked_array()
        lat = lat.copy().to_masked_array()
    except : pass
            
    mask = np.logical_and (np.logical_and(lat>y0, lat<y1), 
            np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)),
                                      np.logical_and(lon-360>x0, lon-360<x1)))
    tab = mmath.where (mask, ptab, np.nan)
    
    return tab

def extend (tab, Lon=False, jplus=25, jpi=None, nperio=4) :
    '''
    Returns extended field eastward to have better plots, and box average crossing the boundary
    Works only for xarray and numpy data (?)

    Useful for vertical sections in OCE and ATM.

    tab : field to extend.
    Lon : (optional, default=False) : if True, add 360 in the extended parts of the field
    jpi : normal longitude dimension of the field. exrtend does nothing it the actual
        size of the field != jpi (avoid to extend several times)
    jplus (optional, default=25) : number of points added on the east side of the field
    
    '''
    mmath = __mmath__ (tab)
    
    if tab.shape[-1] == 1 : extend = tab

    else :
        if jpi == None : jpi = tab.shape[-1]

        if Lon : xplus = -360.0
        else   : xplus =    0.0

        if tab.shape[-1] > jpi :
            extend = tab
        else :
            if nperio == 0 or nperio == 4.2 :
                istart = 0 ; le=jpi+1 ; la=0
            if nperio == 1 :
                istart = 0 ; le=jpi+1 ; la=0
            if nperio == 4 or nperio == 6 : # OPA case with two halo points for periodicity
                istart = 1 ; le=jpi-2 ; la=1  # Perfect, except at the pole that should be masked by lbc_plot
           
            if mmath == xr :
                extend = np.concatenate ((tab.values[..., istart   :istart+le+1    ] + xplus,
                                          tab.values[..., istart+la:istart+la+jplus]         ), axis=-1)
                lon    = tab.dims[-1]
                new_coords = []
                for coord in tab.dims :
                    if coord == lon : new_coords.append ( np.arange( extend.shape[-1]))
                    else            : new_coords.append ( tab.coords[coord].values)
                extend = xr.DataArray ( extend, dims=tab.dims, coords=new_coords )
            else : 
                extend = np.concatenate ((tab [..., istart   :istart+le+1    ] + xplus,
                                          tab [..., istart+la:istart+la+jplus]          ), axis=-1)
    return extend

def orca2reg (ff, lat_name='nav_lat', lon_name='nav_lon', y_name='y', x_name='x') :
    '''
    Assign an ORCA dataset on a regular grid.
    For use in the tropical region.
    
    Inputs : 
      ff : xarray dataset
      lat_name, lon_name : name of latitude and longitude 2D field in ff
      y_name, x_name     : namex of dimensions in ff
      
      Returns : xarray dataset with rectangular grid. Incorrect above 20°N
    '''
    # Compute 1D longitude and latitude
    (lat, lon) = latlon1D (ff[lat_name], ff[lon_name])

    # Assign lon and lat as dimensions of the dataset
    if y_name in ff.dims : 
        lat = xr.DataArray (lat, coords=[lat,], dims=['lat',])     
        ff  = ff.rename_dims ({y_name: "lat",}).assign_coords (lat=lat)
    if x_name in ff.dims :
        lon = xr.DataArray (lon, coords=[lon,], dims=['lon',])
        ff  = ff.rename_dims ({x_name: "lon",}).assign_coords (lon=lon)
    # Force dimensions to be in the right order
    coord_order = ['lat', 'lon']
    for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z',
                 'time_counter', 'time', 'tbnds', 
                 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] :
        if dim in ff.dims : coord_order.insert (0, dim)
        
    ff = ff.transpose (*coord_order)
    return ff

def lbc_init (ptab, nperio=None) :
    '''
    Prepare for all lbc calls
    
    Set periodicity on input field
    nperio    : Type of periodicity
      0       : No periodicity
      1, 4, 6 : Cyclic on i dimension (generaly longitudes) with 2 points halo
      2       : Obsolete (was symmetric condition at southern boundary ?)
      3, 4    : North fold T-point pivot (legacy ORCA2)
      5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo)
    cd_type   : Grid specification : T, U, V or F

    See NEMO documentation for further details
    '''
    jpi = None ; jpj = None
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')
    if ax : jpi = ptab.shape[ix]
    if ay : jpj = ptab.shape[jy]
        
    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio)
    
    if nperio not in nperio_valid_range :
        raise Exception ( f'{nperio=} is not in the valid range {nperio_valid_range}' )

    return jpj, jpi, nperio
        
def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) :
    '''
    Set periodicity on input field
    ptab      : Input array (works for rank 2 at least : ptab[...., lat, lon])
    nperio    : Type of periodicity
    cd_type   : Grid specification : T, U, V or F
    psgn      : For change of sign for vector components (1 for scalars, -1 for vector components)
    
    See NEMO documentation for further details
    '''
    jpj, jpi, nperio = lbc_init (ptab, nperio)
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')
    psgn   = ptab.dtype.type (psgn)
    mmath = __mmath__ (ptab)
    
    if mmath == xr : ztab = ptab.values.copy ()
    else           : ztab = ptab.copy ()
        
    if ax : 
        #
        #> East-West boundary conditions
        # ------------------------------
        if nperio in [1, 4, 6] :
        # ... cyclic
            ztab [...,  0] = ztab [..., -2]
            ztab [..., -1] = ztab [...,  1]

        if ay : 
            #
            #> North-South boundary conditions
            # --------------------------------
            if nperio in [3, 4] :  # North fold T-point pivot
                if cd_type in [ 'T', 'W' ] : # T-, W-point
                    ztab [..., -1, 1:       ] = psgn * ztab [..., -3, -1:0:-1      ]
                    ztab [..., -1, 0        ] = psgn * ztab [..., -3, 2            ]
                    ztab [..., -2, jpi//2:  ] = psgn * ztab [..., -2, jpi//2:0:-1  ]
                    
                if cd_type == 'U' :
                    ztab [..., -1, 0:-1     ] = psgn * ztab [..., -3, -1:0:-1      ]       
                    ztab [..., -1,  0       ] = psgn * ztab [..., -3,  1           ]
                    ztab [..., -1, -1       ] = psgn * ztab [..., -3, -2           ]
                    
                    if nemo_4U_bug :
                        ztab [..., -2, jpi//2+1:-1] = psgn * ztab [..., -2, jpi//2-2:0:-1]
                        ztab [..., -2, jpi//2-1   ] = psgn * ztab [..., -2, jpi//2       ]
                    else :
                        ztab [..., -2, jpi//2-1:-1] = psgn * ztab [..., -2, jpi//2:0:-1]
                        
                if cd_type == 'V' : 
                    ztab [..., -2, 1:       ] = psgn * ztab [..., -3, jpi-1:0:-1   ]
                    ztab [..., -1, 1:       ] = psgn * ztab [..., -4, -1:0:-1      ]   
                    ztab [..., -1, 0        ] = psgn * ztab [..., -4, 2            ]
                    
                if cd_type == 'F' :
                    ztab [..., -2, 0:-1     ] = psgn * ztab [..., -3, -1:0:-1      ]
                    ztab [..., -1, 0:-1     ] = psgn * ztab [..., -4, -1:0:-1      ]
                    ztab [..., -1,  0       ] = psgn * ztab [..., -4,  1           ]
                    ztab [..., -1, -1       ] = psgn * ztab [..., -4, -2           ]
                
            if nperio in [4.2] :  # North fold T-point pivot
                if cd_type in [ 'T', 'W' ] : # T-, W-point
                    ztab [..., -1, jpi//2:  ] = psgn * ztab [..., -1, jpi//2:0:-1  ]
                    
                if cd_type == 'U' :
                    ztab [..., -1, jpi//2-1:-1] = psgn * ztab [..., -1, jpi//2:0:-1]
                    
                if cd_type == 'V' : 
                    ztab [..., -1, 1:       ] = psgn * ztab [..., -2, jpi-1:0:-1   ]
                    
                if cd_type == 'F' :
                    ztab [..., -1, 0:-1     ] = psgn * ztab [..., -2, -1:0:-1      ]
                
            if nperio in [5, 6] :            #  North fold F-point pivot  
                if cd_type in ['T', 'W']  :
                    ztab [..., -1, 0:       ] = psgn * ztab [..., -2, -1::-1       ]
                    
                if cd_type == 'U' :
                    ztab [..., -1, 0:-1     ] = psgn * ztab [..., -2, -2::-1       ]       
                    ztab [..., -1, -1       ] = psgn * ztab [..., -2, 0            ] # Bug ?
                    
                if cd_type == 'V' :
                    ztab [..., -1, 0:       ] = psgn * ztab [..., -3, -1::-1       ]
                    ztab [..., -2, jpi//2:  ] = psgn * ztab [..., -2, jpi//2-1::-1 ]
                    
                if cd_type == 'F' :
                    ztab [..., -1, 0:-1     ] = psgn * ztab [..., -3, -2::-1       ]
                    ztab [..., -1, -1       ] = psgn * ztab [..., -3, 0            ]
                    ztab [..., -2, jpi//2:-1] = psgn * ztab [..., -2, jpi//2-2::-1 ]
                    
            #
            #> East-West boundary conditions
            # ------------------------------
            if nperio in [1, 4, 6] :
                # ... cyclic
                ztab [...,  0] = ztab [..., -2]
                ztab [..., -1] = ztab [...,  1]

    if mmath == xr :
        ztab = xr.DataArray ( ztab, dims=ptab.dims, coords=ptab.coords )
        ztab.attrs = ptab.attrs
        
    return ztab

def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) :
    #
    '''
    Mask fields on duplicated points
    ptab      : Input array. Rank 2 at least : ptab [...., lat, lon]
    nperio    : Type of periodicity
    cd_type   : Grid specification : T, U, V or F
    
    See NEMO documentation for further details
    '''
    jpj, jpi, nperio = lbc_init (ptab, nperio)
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')
    ztab = ptab.copy ()

    if ax : 
        #
        #> East-West boundary conditions
        # ------------------------------
        if nperio in [1, 4, 6] :
        # ... cyclic
            ztab [...,  0] = sval
            ztab [..., -1] = sval

        if ay : 
            #
            #> South (in which nperio cases ?)
            # --------------------------------
            if nperio in [1, 3, 4, 5, 6] :
                ztab [..., 0, :] = sval
        
            #
            #> North-South boundary conditions
            # --------------------------------
            if nperio in [3, 4] :  # North fold T-point pivot
                if cd_type in [ 'T', 'W' ] : # T-, W-point
                    ztab [..., -1,  :         ] = sval
                    ztab [..., -2, :jpi//2  ] = sval
                
                if cd_type == 'U' :
                    ztab [..., -1,  :         ] = sval  
                    ztab [..., -2, jpi//2+1:  ] = sval
                
                if cd_type == 'V' :
                    ztab [..., -2, :       ] = sval
                    ztab [..., -1, :       ] = sval   

                if cd_type == 'F' :
                    ztab [..., -2, :       ] = sval
                    ztab [..., -1, :       ] = sval

            if nperio in [4.2] :  # North fold T-point pivot
                if cd_type in [ 'T', 'W' ] : # T-, W-point
                    ztab [..., -1, jpi//2  :  ] = sval

                if cd_type == 'U' :
                    ztab [..., -1, jpi//2-1:-1] = sval

                if cd_type == 'V' : 
                    ztab [..., -1, 1:       ] = sval

                if cd_type == 'F' :
                    ztab [..., -1, 0:-1     ] = sval

            if nperio in [5, 6] :            #  North fold F-point pivot
                if cd_type in ['T', 'W']  :
                    ztab [..., -1, 0:       ] = sval

                if cd_type == 'U' :
                    ztab [..., -1, 0:-1     ] = sval       
                    ztab [..., -1, -1       ] = sval

                if cd_type == 'V' :
                    ztab [..., -1, 0:       ] = sval
                    ztab [..., -2, jpi//2:  ] = sval

                if cd_type == 'F' :
                    ztab [..., -1, 0:-1       ] = sval
                    ztab [..., -1, -1         ] = sval
                    ztab [..., -2, jpi//2+1:-1] = sval

    return ztab

def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) :
    '''
    Set periodicity on input field, adapted for plotting for any cartopy projection
    ptab      : Input array. Rank 2 at least : ptab[...., lat, lon]
    nperio    : Type of periodicity
    cd_type   : Grid specification : T, U, V or F
    psgn      : For change of sign for vector components (1 for scalars, -1 for vector components)
    
    See NEMO documentation for further details
    '''
    jpj, jpi, nperio = lbc_init (ptab, nperio)
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')
    psgn   = ptab.dtype.type (psgn)
    ztab   = ptab.copy ()

    if ax : 
        #
        #> East-West boundary conditions
        # ------------------------------
        if nperio in [1, 4, 6] :
            # ... cyclic
            ztab [..., :,  0] = ztab [..., :, -2]
            ztab [..., :, -1] = ztab [..., :,  1]

        if ay : 
            #> Masks south
            # ------------
            if nperio in [4, 6] : ztab [..., 0, : ] = sval

            #
            #> North-South boundary conditions
            # --------------------------------
            if nperio in [3, 4] :  # North fold T-point pivot
                if cd_type in [ 'T', 'W' ] : # T-, W-point
                    ztab [..., -1,  :      ] = sval
                    #ztab [..., -2, jpi//2: ] = sval
                    ztab [..., -2, :jpi//2 ] = sval # Give better plots than above
                if cd_type == 'U' :
                    ztab [..., -1, : ] = sval

                if cd_type == 'V' : 
                    ztab [..., -2, : ] = sval
                    ztab [..., -1, : ] = sval

                if cd_type == 'F' :
                    ztab [..., -2, : ] = sval
                    ztab [..., -1, : ] = sval

            if nperio in [4.2] :  # North fold T-point pivot
                if cd_type in [ 'T', 'W' ] : # T-, W-point
                    ztab [..., -1, jpi//2:  ] = sval

                if cd_type == 'U' :
                    ztab [..., -1, jpi//2-1:-1] = sval

                if cd_type == 'V' : 
                    ztab [..., -1, 1:       ] = sval

                if cd_type == 'F' :
                    ztab [..., -1, 0:-1     ] = sval

            if nperio in [5, 6] :            #  North fold F-point pivot  
                if cd_type in ['T', 'W']  :
                    ztab [..., -1, : ] = sval

                if cd_type == 'U' :
                    ztab [..., -1, : ] = sval      

                if cd_type == 'V' :
                    ztab [..., -1, :        ] = sval
                    ztab [..., -2, jpi//2:  ] = sval

                if cd_type == 'F' :
                    ztab [..., -1, :          ] = sval
                    ztab [..., -2, jpi//2+1:-1] = sval

    return ztab

def lbc_add (ptab, nperio=None, cd_type=None, psgn=1, sval=None) :
    '''
    Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2
      Peridodicity halo has been removed
    This routine adds the halos if needed

    ptab      : Input array (works 
      rank 2 at least : ptab[...., lat, lon]
    nperio    : Type of periodicity
 
    See NEMO documentation for further details
    '''
    mmath = __mmath__ (ptab) 
    jpj, jpi, nperio = lbc_init (ptab, nperio)
    lshape = get_shape (ptab)
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')

    t_shape = np.array (ptab.shape)

    if nperio == 4.2 or nperio == 6.2 :
      
        ext_shape = t_shape.copy()
        if 'X' in lshape : ext_shape[ix] = ext_shape[ix] + 2
        if 'Y' in lshape : ext_shape[jy] = ext_shape[jy] + 1

        if mmath == xr :
            ptab_ext = xr.DataArray (np.zeros (ext_shape), dims=ptab.dims)
            if 'X' in lshape and 'Y' in lshape :
                ptab_ext.values[..., :-1, 1:-1] = ptab.values.copy ()
            else :
                if 'X' in lshape     : ptab_ext.values[...,      1:-1] = ptab.values.copy ()
                if 'Y' in lshape     : ptab_ext.values[..., :-1      ] = ptab.values.copy ()
        else           :
            ptab_ext =               np.zeros (ext_shape)
            if 'X' in lshape and 'Y' in lshape : ptab_ext       [..., :-1, 1:-1] = ptab.copy ()
            else :
                if 'X' in lshape     : ptab_ext       [...,      1:-1] = ptab.copy ()
                if 'Y' in lshape     : ptab_ext       [..., :-1      ] = ptab.copy ()            

        if nperio == 4.2 : ptab_ext = lbc (ptab_ext, nperio=4, cd_type=cd_type, psgn=psgn)
        if nperio == 6.2 : ptab_ext = lbc (ptab_ext, nperio=6, cd_type=cd_type, psgn=psgn)
        
        if mmath == xr :
            ptab_ext.attrs = ptab.attrs
            kz, az = __findAxis__ (ptab, 'z')
            it, at = __findAxis__ (ptab, 't')
            if az : ptab_ext = ptab_ext.assign_coords ( {az:ptab.coords[az]} )
            if at : ptab_ext = ptab_ext.assign_coords ( {at:ptab.coords[at]} )

    else : ptab_ext = lbc (ptab, nperio=nperio, cd_type=cd_type, psgn=psgn)
        
    return ptab_ext

def lbc_del (ptab, nperio=None, cd_type='T', psgn=1) :
    '''
    Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2
      Periodicity halo has been removed
    This routine removes the halos if needed

    ptab      : Input array (works 
      rank 2 at least : ptab[...., lat, lon]
    nperio    : Type of periodicity
 
    See NEMO documentation for further details
    '''
    jpj, jpi, nperio = lbc_init (ptab, nperio)
    lshape = get_shape (ptab)
    ix, ax = __findAxis__ (ptab, 'x')
    jy, ay = __findAxis__ (ptab, 'y')

    if nperio == 4.2 or nperio == 6.2 :
        if ax or ay : 
            if ax and ay : 
                return lbc (ptab[..., :-1, 1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn)
            else : 
                if ax :
                    return lbc (ptab[...,      1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn)
                if ay :
                    return lbc (ptab[..., -1], nperio=nperio, cd_type=cd_type, psgn=psgn)
        else :
            return ptab
    else :
        return ptab

def lbc_index (jj, ii, jpj, jpi, nperio=None, cd_type='T') :
    '''
    For indexes of a NEMO point, give the corresponding point inside the util domain
    jj, ii    : indexes
    jpi, jpi  : size of domain
    nperio    : type of periodicity
    cd_type   : grid specification : T, U, V or F
    
    See NEMO documentation for further details
    '''

    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio)
    
    ## For the sake of simplicity, switch to the convention of original lbc Fortran routine from NEMO
    ## : starts indexes at 1
    jy = jj + 1 ; ix = ii + 1

    mmath = __mmath__ (jj)
    if mmath == None : mmath=np

    #
    #> East-West boundary conditions
    # ------------------------------
    if nperio in [1, 4, 6] :
        #... cyclic
        ix = mmath.where (ix==jpi, 2   , ix)
        ix = mmath.where (ix== 1 ,jpi-1, ix)

    #
    def modIJ (cond, jy_new, ix_new) :
        jy_r = mmath.where (cond, jy_new, jy)
        ix_r = mmath.where (cond, ix_new, ix)
        return jy_r, ix_r
    #
    #> North-South boundary conditions
    # --------------------------------
    if nperio in [ 3 , 4 ]  :
        if cd_type in  [ 'T' , 'W' ] :
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix>=2       ), jpj-2, jpi-ix+2)
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1       ), jpj-1, 3       )   
            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy   , jpi-ix+2) 

        if cd_type in [ 'U' ] :
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , np.logical_and (ix>=1, ix <= jpi-1)   ), jy   , jpi-ix+1)
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1  )                               , jpj-2, 2       )
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==jpi)                               , jpj-2, jpi-1   )
            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, np.logical_and (ix>=jpi//2, ix<=jpi-1)), jy   , jpi-ix+1)
          
        if cd_type in [ 'V' ] :
            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=2  ), jpj-2, jpi-ix+2)
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix>=2  ), jpj-3, jpi-ix+2)
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1  ), jpj-3,  3      )
            
        if cd_type in [ 'F' ] :
            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix<=jpi-1), jpj-2, jpi-ix+1)
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix<=jpi-1), jpj-3, jpi-ix+1)
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1    ), jpj-3, 2       )
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==jpi  ), jpj-3, jpi-1   )

    if nperio in [ 5 , 6 ] :
        if cd_type in [ 'T' , 'W' ] :                        # T-, W-point
             (jy, ix) = modIJ (jy==jpj, jpj-1, jpi-ix+1)
 
        if cd_type in [ 'U' ] :                              # U-point
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix<=jpi-1   ), jpj-1, jpi-ix  )
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==jpi     ), jpi-1, 1       )
            
        if cd_type in [ 'V' ] :    # V-point
            (jy, ix) = modIJ (jy==jpj                                 , jy   , jpi-ix+1)
            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy   , jpi-ix+1)
            
        if cd_type in [ 'F' ] :                              # F-point
            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix<=jpi-1   ), jpj-2, jpi-ix  )
            (jy, ix) = modIJ (np.logical_and (ix==jpj  , ix==jpi     ), jpj-2, 1       )
            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy   , jpi-ix  )

    ## Restore convention to Python/C : indexes start at 0
    jy += -1 ; ix += -1

    if isinstance (jj, int) : jy = jy.item ()
    if isinstance (ii, int) : ix = ix.item ()

    return jy, ix
    
def findJI (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False, out=None) :
    '''
    Description: seeks J,I indices of the grid point which is the closest of a given point 
    Usage: go FindJI  <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask]
    <longitude fields> <latitude field> are 2D fields on J/I (Y/X) dimensions
    mask : if given, seek only non masked grid points (i.e with mask=1)
    
    Example : findIJ (40, -20, nav_lat, nav_lon, mask=1.0)

    Note : all longitudes and latitudes in degrees
        
    Note : may work with 1D lon/lat (?)
    '''
    # Get grid dimensions
    if len (lon_grid.shape) == 2 : (jpj, jpi) = lon_grid.shape
    else                         : jpj = len(lat_grid) ; jpi=len(lon_grid)

    mmath = __mmath__ (lat_grid)
        
    # Compute distance from the point to all grid points (in radian)
    arg      = np.sin (rad*lat_data) * np.sin (rad*lat_grid) \
             + np.cos (rad*lat_data) * np.cos (rad*lat_grid) * np.cos(rad*(lon_data-lon_grid))
    distance = np.arccos (arg) + 4.0*rpi*(1.0-mask) # Send masked points to 'infinite' 

    # Truncates to alleviate some precision problem with some grids
    prec = int (1E7)
    distance = (distance*prec).astype(int) / prec

    # Compute minimum of distance, and index of minimum
    #
    distance_min = distance.min    ()
    jimin        = int (distance.argmin ())
    
    # Compute 2D indices 
    jmin = jimin // jpi ; imin = jimin - jmin*jpi
    
    # Result
    if verbose :
        # Compute distance achieved
        mindist = distance [jmin, imin]
        
        # Compute azimuth
        dlon = lon_data-lon_grid[jmin,imin]
        arg  = np.sin (rad*dlon) /  (np.cos(rad*lat_data)*np.tan(rad*lat_grid[jmin,imin]) - np.sin(rad*lat_data)*np.cos(rad*dlon))
        azimuth = dar*np.arctan (arg)
        print ( f'I={imin:d} J={jmin:d} - Data:{lat_data:5.1f}°N {lon_data:5.1f}°E - Grid:{lat_grid[jmin,imin]:4.1f}°N '   \
            +   f'{lon_grid[jmin,imin]:4.1f}°E - Dist: {ra*distance[jmin,imin]:6.1f}km {dar*distance[jmin,imin]:5.2f}° ' \
            +   f'- Azimuth: {rad*azimuth:3.2f}rad - {azimuth:5.1f}°' )

    if   out=='dict'                               : return {'x':imin, 'y':jmin}
    elif out=='array' or out=='numpy'  or out=='np': return np.array ( [jmin, imin] )
    elif out=='xarray' or out=='xr'                : return xr.DataArray ( [jmin, imin] )
    elif out=='list'                               : return [jmin, imin]
    elif out=='tuple'                              : return jmin, imin
    else                                           : return jmin, imin

def geo2en (pxx, pyy, pzz, glam, gphi) : 
    '''
    Change vector from geocentric to east/north

    Inputs :
        pxx, pyy, pzz : components on the geocentric system
        glam, gphi : longitude and latitude of the points
    '''

    gsinlon = np.sin (rad * glam)
    gcoslon = np.cos (rad * glam)
    gsinlat = np.sin (rad * gphi)
    gcoslat = np.cos (rad * gphi)
          
    pte = - pxx * gsinlon            + pyy * gcoslon
    ptn = - pxx * gcoslon * gsinlat  - pyy * gsinlon * gsinlat + pzz * gcoslat

    return pte, ptn

def en2geo (pte, ptn, glam, gphi) :
    '''
    Change vector from east/north to geocentric

    Inputs : 
        pte, ptn   : eastward/northward components
        glam, gphi : longitude and latitude of the points
    '''
    
    gsinlon = np.sin (rad * glam)
    gcoslon = np.cos (rad * glam)
    gsinlat = np.sin (rad * gphi)
    gcoslat = np.cos (rad * gphi)

    pxx = - pte * gsinlon - ptn * gcoslon * gsinlat
    pyy =   pte * gcoslon - ptn * gsinlon * gsinlat
    pzz =   ptn * gcoslat
    
    return pxx, pyy, pzz


def clo_lon (lon, lon0=0., rad=False, deg=True) :
    '''Choose closest to lon0 longitude, adding or substacting 360° if needed'''
    mmath = __mmath__ (lon, np)
    if rad : lon_range = 2.*np.pi
    if deg : lon_range = 360.
    clo_lon = lon
    clo_lon = mmath.where (clo_lon > lon0 + lon_range*0.5, clo_lon-lon_range, clo_lon)
    clo_lon = mmath.where (clo_lon < lon0 - lon_range*0.5, clo_lon+lon_range, clo_lon)
    clo_lon = mmath.where (clo_lon > lon0 + lon_range*0.5, clo_lon-lon_range, clo_lon)
    clo_lon = mmath.where (clo_lon < lon0 - lon_range*0.5, clo_lon+lon_range, clo_lon)
    if clo_lon.shape == () : clo_lon = clo_lon.item ()
    if mmath == xr :
        try : 
            for attr in lon.attrs : clo_lon.attrs[attr] = lon.attrs[attr]
        except :
            pass
    return clo_lon


def index2depth (pk, gdept_0) :
    '''
    From index (real, continuous), get depth
    '''
    jpk = gdept_0.shape[0]
    kk = xr.DataArray(pk)
    k  = np.maximum (0, np.minimum (jpk-1, kk    ))
    k0 = np.floor (k).astype (int)
    k1 = np.maximum (0, np.minimum (jpk-1,  k0+1))
    zz = k - k0
    gz = (1.0-zz)*gdept_0[k0]+ zz*gdept_0[k1]
    return gz.values

def depth2index (pz, gdept_0) :
    '''
    From depth, get index (real, continuous)
    '''
    jpk  = gdept_0.shape[0]
    if type (pz) == xr.core.dataarray.DataArray :
        zz   = xr.DataArray (pz.values, dims=('zz',))
    elif type (pz) == np.ndarray :
        zz   = xr.DataArray (pz.ravel(), dims=('zz',))
    else :
        zz   = xr.DataArray (np.array([pz]).ravel(), dims=('zz',))
    zz   = np.minimum (gdept_0[-1], np.maximum (0, zz))
    
    idk1 = np.minimum ( (gdept_0-zz), 0.).argmax (axis=0).astype(int)
    idk1 = np.maximum (0, np.minimum (jpk-1,  idk1  ))
    idk2 = np.maximum (0, np.minimum (jpk-1,  idk1-1))
    
    ff = (zz - gdept_0[idk2])/(gdept_0[idk1]-gdept_0[idk2])
    idk = ff*idk1 + (1.0-ff)*idk2
    idk = xr.where ( np.isnan(idk), idk1, idk)
    return idk.values

def index2depth_panels (pk, gdept_0, depth0, fact) :
    '''
    From  index (real, continuous), get depth, with bottom part compressed
    '''
    jpk = gdept_0.shape[0]
    kk = xr.DataArray (pk)
    k  = np.maximum (0, np.minimum (jpk-1, kk    ))
    k0 = np.floor (k).astype (int)
    k1 = np.maximum (0, np.minimum (jpk-1,  k0+1))
    zz = k - k0
    gz = (1.0-zz)*gdept_0[k0]+ zz*gdept_0[k1]
    gz = xr.where ( gz<depth0, gz, depth0 + (gz-depth0)*fact)
    return gz.values

def depth2index_panels (pz, gdept_0, depth0, fact) :
    '''
    From  index (real, continuous), get depth, with bottom part compressed
    '''
    jpk = gdept_0.shape[0]
    if type (pz) == xr.core.dataarray.DataArray :
        zz   = xr.DataArray (pz.values , dims=('zz',))
    elif type (pz) == np.ndarray :
        zz   = xr.DataArray (pz.ravel(), dims=('zz',))
    else : 
        zz   = xr.DataArray (np.array([pz]).ravel(), dims=('zz',))
    zz         = np.minimum (gdept_0[-1], np.maximum (0, zz))
    gdept_comp = xr.where ( gdept_0>depth0, (gdept_0-depth0)*fact+depth0, gdept_0)
    zz_comp    = xr.where ( zz     >depth0, (zz     -depth0)*fact+depth0, zz     )
    zz_comp    = np.minimum (gdept_comp[-1], np.maximum (0, zz_comp))

    idk1 = np.minimum ( (gdept_0-zz_comp), 0.).argmax (axis=0).astype(int)
    idk1 = np.maximum (0, np.minimum (jpk-1,  idk1  ))
    idk2 = np.maximum (0, np.minimum (jpk-1,  idk1-1))
     
    ff = (zz_comp - gdept_0[idk2])/(gdept_0[idk1]-gdept_0[idk2])
    idk = ff*idk1 + (1.0-ff)*idk2
    idk = xr.where ( np.isnan(idk), idk1, idk)
    return idk.values

def depth2comp (pz, depth0, fact ) :
    '''
    Form depth, get compressed depth, with bottom part compressed
    '''
    #print ('start depth2comp')
    if type (pz) == xr.core.dataarray.DataArray :
        zz   = pz.values
    elif type(pz) == list :
        zz = np.array (pz)
    else : 
        zz   = pz
    gz = np.where ( zz>depth0, (zz-depth0)*fact+depth0, zz)
    #print ( f'depth2comp : {gz=}' )
    if type (pz) in [int, float] : return gz.item()
    else : return gz
    #print ('end depth2comp')

def comp2depth (pz, depth0, fact ) :
    '''
    Form compressed depth, get depth, with bottom part compressed
    '''
    if type (pz) == xr.core.dataarray.DataArray :
        zz   = pz.values
    elif type(pz) == list :
        zz = np.array (pz)
    else : 
        zz   = pz
    gz = np.where ( zz>depth0, (zz-depth0)/fact+depth0, zz)
    if type (pz) in [int, float] : return gz.item()
    else : return gz

def index2lon (pi, lon1D) :
    '''
    From index (real, continuous), get longitude
    '''
    jpi = lon1D.shape[0]
    ii = xr.DataArray (pi)
    i =  np.maximum (0, np.minimum (jpi-1, ii    ))
    i0 = np.floor (i).astype (int)
    i1 = np.maximum (0, np.minimum (jpi-1,  i0+1))
    xx = i - i0
    gx = (1.0-xx)*lon1D[i0]+ xx*lon1D[i1]
    return gx.values

def lon2index (px, lon1D) :
    '''
    From longitude, get index (real, continuous)
    '''
    jpi  = lon1D.shape[0]
    if type (px) == xr.core.dataarray.DataArray :
        xx   = xr.DataArray (px.values , dims=('xx',))
    elif type (px) == np.ndarray :
        xx   = xr.DataArray (px.ravel(), dims=('xx',))
    else : 
        xx   = xr.DataArray (np.array([px]).ravel(), dims=('xx',))
    xx   = xr.where ( xx>lon1D.max(), xx-360.0, xx)
    xx   = xr.where ( xx<lon1D.min(), xx+360.0, xx)
    xx   = np.minimum (lon1D.max(), np.maximum(xx, lon1D.min() ))
    idi1 = np.minimum ( (lon1D-xx), 0.).argmax (axis=0).astype(int)
    idi1 = np.maximum (0, np.minimum (jpi-1,  idi1  ))
    idi2 = np.maximum (0, np.minimum (jpi-1,  idi1-1))
    
    ff = (xx - lon1D[idi2])/(lon1D[idi1]-lon1D[idi2])
    idi = ff*idi1 + (1.0-ff)*idi2
    idi = xr.where ( np.isnan(idi), idi1, idi)
    return idi.values

def index2lat (pj, lat1D) :
    '''
    From index (real, continuous), get latitude
    '''
    jpj = lat1D.shape[0]
    jj  = xr.DataArray (pj)
    j   = np.maximum (0, np.minimum (jpj-1, jj    ))
    j0  = np.floor (j).astype (int)
    j1  = np.maximum (0, np.minimum (jpj-1,  j0+1))
    yy  = j - j0
    gy  = (1.0-yy)*lat1D[j0]+ yy*lat1D[j1]
    return gy.values

def lat2index (py, lat1D) :
    '''
    From latitude, get index (real, continuous)
    '''
    jpj = lat1D.shape[0]
    if type (py) == xr.core.dataarray.DataArray :
        yy   = xr.DataArray (py.values , dims=('yy',))
    elif type (py) == np.ndarray :
        yy   = xr.DataArray (py.ravel(), dims=('yy',))
    else : 
        yy   = xr.DataArray (np.array([py]).ravel(), dims=('yy',))
    yy   = np.minimum (lat1D.max(), np.minimum(yy, lat1D.max() ))
    idj1 = np.minimum ( (lat1D-yy), 0.).argmax (axis=0).astype(int)
    idj1 = np.maximum (0, np.minimum (jpj-1,  idj1  ))
    idj2 = np.maximum (0, np.minimum (jpj-1,  idj1-1))
    
    ff = (yy - lat1D[idj2])/(lat1D[idj1]-lat1D[idj2])
    idj = ff*idj1 + (1.0-ff)*idj2
    idj = xr.where ( np.isnan(idj), idj1, idj)
    return idj.values

def angle_full (glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif, nperio=None) :
    '''Compute sinus and cosinus of model line direction with respect to east'''
    mmath = __mmath__ (glamt)

    zlamt = lbc_add (glamt, nperio, 'T', 1.)
    zphit = lbc_add (gphit, nperio, 'T', 1.)
    zlamu = lbc_add (glamu, nperio, 'U', 1.)
    zphiu = lbc_add (gphiu, nperio, 'U', 1.)
    zlamv = lbc_add (glamv, nperio, 'V', 1.)
    zphiv = lbc_add (gphiv, nperio, 'V', 1.)
    zlamf = lbc_add (glamf, nperio, 'F', 1.)
    zphif = lbc_add (gphif, nperio, 'F', 1.)
    
    # north pole direction & modulous (at T-point)
    zxnpt = 0. - 2.0 * np.cos (rad*zlamt) * np.tan (rpi/4.0 - rad*zphit/2.0)
    zynpt = 0. - 2.0 * np.sin (rad*zlamt) * np.tan (rpi/4.0 - rad*zphit/2.0)
    znnpt = zxnpt*zxnpt + zynpt*zynpt
    
    # north pole direction & modulous (at U-point)
    zxnpu = 0. - 2.0 * np.cos (rad*zlamu) * np.tan (rpi/4.0 - rad*zphiu/2.0)
    zynpu = 0. - 2.0 * np.sin (rad*zlamu) * np.tan (rpi/4.0 - rad*zphiu/2.0)
    znnpu = zxnpu*zxnpu + zynpu*zynpu
    
    # north pole direction & modulous (at V-point)
    zxnpv = 0. - 2.0 * np.cos (rad*zlamv) * np.tan (rpi/4.0 - rad*zphiv/2.0)
    zynpv = 0. - 2.0 * np.sin (rad*zlamv) * np.tan (rpi/4.0 - rad*zphiv/2.0)
    znnpv = zxnpv*zxnpv + zynpv*zynpv

    # north pole direction & modulous (at F-point)
    zxnpf = 0. - 2.0 * np.cos( rad*zlamf ) * np.tan ( rpi/4. - rad*zphif/2. )
    zynpf = 0. - 2.0 * np.sin( rad*zlamf ) * np.tan ( rpi/4. - rad*zphif/2. )
    znnpf = zxnpf*zxnpf + zynpf*zynpf

    # j-direction: v-point segment direction (around T-point)
    zlam = zlamv  
    zphi = zphiv
    zlan = np.roll ( zlamv, axis=-2, shift=1)  # glamv (ji,jj-1)
    zphh = np.roll ( zphiv, axis=-2, shift=1)  # gphiv (ji,jj-1)
    zxvvt =  2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    zyvvt =  2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    znvvt = np.sqrt ( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )

    # j-direction: f-point segment direction (around u-point)
    zlam = zlamf
    zphi = zphif
    zlan = np.roll (zlamf, axis=-2, shift=1) # glamf (ji,jj-1)
    zphh = np.roll (zphif, axis=-2, shift=1) # gphif (ji,jj-1)
    zxffu =  2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    zyffu =  2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    znffu = np.sqrt ( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )

    # i-direction: f-point segment direction (around v-point)
    zlam = zlamf  
    zphi = zphif
    zlan = np.roll (zlamf, axis=-1, shift=1) # glamf (ji-1,jj)
    zphh = np.roll (zphif, axis=-1, shift=1) # gphif (ji-1,jj)
    zxffv =  2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    zyffv =  2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    znffv = np.sqrt ( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )

    # j-direction: u-point segment direction (around f-point)
    zlam = np.roll (zlamu, axis=-2, shift=-1) # glamu (ji,jj+1) 
    zphi = np.roll (zphiu, axis=-2, shift=-1) # gphiu (ji,jj+1)
    zlan = zlamu
    zphh = zphiu
    zxuuf =  2. * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2. * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    zyuuf =  2. * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \
          -  2. * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )
    znuuf = np.sqrt ( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )

    
    # cosinus and sinus using scalar and vectorial products
    gsint = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt
    gcost = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt
    
    gsinu = ( zxnpu*zyffu - zynpu*zxffu ) / znffu
    gcosu = ( zxnpu*zxffu + zynpu*zyffu ) / znffu
    
    gsinf = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf
    gcosf = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf
    
    gsinv = ( zxnpv*zxffv + zynpv*zyffv ) / znffv
    gcosv =-( zxnpv*zyffv - zynpv*zxffv ) / znffv  # (caution, rotation of 90 degres)
    
    #gsint = lbc (gsint, cd_type='T', nperio=nperio, psgn=-1.)
    #gcost = lbc (gcost, cd_type='T', nperio=nperio, psgn=-1.)
    #gsinu = lbc (gsinu, cd_type='U', nperio=nperio, psgn=-1.)
    #gcosu = lbc (gcosu, cd_type='U', nperio=nperio, psgn=-1.)
    #gsinv = lbc (gsinv, cd_type='V', nperio=nperio, psgn=-1.)
    #gcosv = lbc (gcosv, cd_type='V', nperio=nperio, psgn=-1.)
    #gsinf = lbc (gsinf, cd_type='F', nperio=nperio, psgn=-1.)
    #gcosf = lbc (gcosf, cd_type='F', nperio=nperio, psgn=-1.)

    gsint = lbc_del (gsint, cd_type='T', nperio=nperio, psgn=-1.)
    gcost = lbc_del (gcost, cd_type='T', nperio=nperio, psgn=-1.)
    gsinu = lbc_del (gsinu, cd_type='U', nperio=nperio, psgn=-1.)
    gcosu = lbc_del (gcosu, cd_type='U', nperio=nperio, psgn=-1.)
    gsinv = lbc_del (gsinv, cd_type='V', nperio=nperio, psgn=-1.)
    gcosv = lbc_del (gcosv, cd_type='V', nperio=nperio, psgn=-1.)
    gsinf = lbc_del (gsinf, cd_type='F', nperio=nperio, psgn=-1.)
    gcosf = lbc_del (gcosf, cd_type='F', nperio=nperio, psgn=-1.)

    if mmath == xr :
        gsint = gsint.assign_coords ( glamt.coords )
        gcost = gcost.assign_coords ( glamt.coords )
        gsinu = gsinu.assign_coords ( glamu.coords )
        gcosu = gcosu.assign_coords ( glamu.coords )
        gsinv = gsinv.assign_coords ( glamv.coords )
        gcosv = gcosv.assign_coords ( glamv.coords )
        gsinf = gsinf.assign_coords ( glamf.coords )
        gcosf = gcosf.assign_coords ( glamf.coords )

    return gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf

def angle (glam, gphi, nperio, cd_type='T') :
    '''Compute sinus and cosinus of model line direction with respect to east'''
    mmath = __mmath__ (glam)

    zlam = lbc_add (glam, nperio, cd_type, 1.)
    zphi = lbc_add (gphi, nperio, cd_type, 1.)
    
    # north pole direction & modulous
    zxnp = 0. - 2.0 * np.cos (rad*zlam) * np.tan (rpi/4.0 - rad*zphi/2.0)
    zynp = 0. - 2.0 * np.sin (rad*zlam) * np.tan (rpi/4.0 - rad*zphi/2.0)
    znnp = zxnp*zxnp + zynp*zynp

    # j-direction: segment direction (around point)
    zlan_n = np.roll (zlam, axis=-2, shift=-1) # glam [jj+1, ji]
    zphh_n = np.roll (zphi, axis=-2, shift=-1) # gphi [jj+1, ji]
    zlan_s = np.roll (zlam, axis=-2, shift= 1) # glam [jj-1, ji]
    zphh_s = np.roll (zphi, axis=-2, shift= 1) # gphi [jj-1, ji]
    
    zxff = 2.0 * np.cos (rad*zlan_n) * np.tan (rpi/4.0 - rad*zphh_n/2.0) \
        -  2.0 * np.cos (rad*zlan_s) * np.tan (rpi/4.0 - rad*zphh_s/2.0)
    zyff = 2.0 * np.sin (rad*zlan_n) * np.tan (rpi/4.0 - rad*zphh_n/2.0) \
        -  2.0 * np.sin (rad*zlan_s) * np.tan (rpi/4.0 - rad*zphh_s/2.0)
    znff = np.sqrt (znnp * (zxff*zxff + zyff*zyff) )
 
    gsin = (zxnp*zyff - zynp*zxff) / znff
    gcos = (zxnp*zxff + zynp*zyff) / znff

    gsin = lbc_del (gsin, cd_type=cd_type, nperio=nperio, psgn=-1.)
    gcos = lbc_del (gcos, cd_type=cd_type, nperio=nperio, psgn=-1.)

    if mmath == xr :
        gsin = gsin.assign_coords ( glam.coords )
        gcos = gcos.assign_coords ( glam.coords )
        
    return gsin, gcos

def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type ) :
    '''
    ** Purpose :   Rotate the Repere: Change vector componantes between
    geographic grid --> stretched coordinates grid.
    All components are on the same grid (T, U, V or F) 
    '''

    u_i = + u_e * gcos + v_n * gsin
    v_j = - u_e * gsin + v_n * gcos
    
    u_i = lbc (u_i, nperio=nperio, cd_type=cd_type, psgn=-1.0)
    v_j = lbc (v_j, nperio=nperio, cd_type=cd_type, psgn=-1.0)
    
    return u_i, v_j

def rot_ij2en ( u_i, v_j, gsin, gcos, nperio, cd_type='T' ) :
    '''
    ** Purpose :   Rotate the Repere: Change vector componantes from
    stretched coordinates grid --> geographic grid
    All components are on the same grid (T, U, V or F) 
    '''
    u_e = + u_i * gcos - v_j * gsin
    v_n = + u_i * gsin + v_j * gcos
    
    u_e = lbc (u_e, nperio=nperio, cd_type=cd_type, psgn=1.0)
    v_n = lbc (v_n, nperio=nperio, cd_type=cd_type, psgn=1.0)
    
    return u_e, v_n

def rot_uv2en ( uo, vo, gsint, gcost, nperio, zdim=None ) :
    '''
    ** Purpose :   Rotate the Repere: Change vector componantes from
    stretched coordinates grid --> geographic grid
    uo is on the U grid point, vo is on the V grid point
    east-north components on the T grid point   
    '''
    mmath = __mmath__ (uo)
    
    ut = U2T (uo, nperio=nperio, psgn=-1.0, zdim=zdim)
    vt = V2T (vo, nperio=nperio, psgn=-1.0, zdim=zdim)
    
    u_e = + ut * gcost - vt * gsint
    v_n = + ut * gsint + vt * gcost

    u_e = lbc (u_e, nperio=nperio, cd_type='T', psgn=1.0)
    v_n = lbc (v_n, nperio=nperio, cd_type='T', psgn=1.0)
    
    return u_e, v_n

def rot_uv2enF ( uo, vo, gsinf, gcosf, nperio, zdim=None ) :
    '''
    ** Purpose : Rotate the Repere: Change vector componantes from
    stretched coordinates grid --> geographic grid
    uo is on the U grid point, vo is on the V grid point
    east-north components on the T grid point   
    '''
    mmath = __mmath__ (uo)

    uf = U2F (uo, nperio=nperio, psgn=-1.0, zdim=zdim)
    vf = V2F (vo, nperio=nperio, psgn=-1.0, zdim=zdim)
    
    u_e = + uf * gcosf - vf * gsinf
    v_n = + uf * gsinf + vf * gcosf

    u_e = lbc (u_e, nperio=nperio, cd_type='F', psgn= 1.0)
    v_n = lbc (v_n, nperio=nperio, cd_type='F', psgn= 1.0)
    
    return u_e, v_n

def U2T (utab, nperio=None, psgn=-1.0, zdim=None, action='ave') :
    '''Interpolate an array from U grid to T grid i-mean)'''
    mmath = __mmath__ (utab)
    utab_0 = mmath.where ( np.isnan(utab), 0., utab)
    lperio, aperio = lbc_diag (nperio)
    utab_0 = lbc_add (utab_0, nperio=nperio, cd_type='U', psgn=psgn)
    ix, ax = __findAxis__ (utab_0, 'x')
    kz, az = __findAxis__ (utab_0, 'z')

    if ax : 
        if action == 'ave' : ttab = 0.5 *      (utab_0 + np.roll (utab_0, axis=ix, shift=1))
        if action == 'min' : ttab = np.minimum (utab_0 , np.roll (utab_0, axis=ix, shift=1))
        if action == 'max' : ttab = np.maximum (utab_0 , np.roll (utab_0, axis=ix, shift=1))
        if action == 'mult': ttab =             utab_0 * np.roll (utab_0, axis=ix, shift=1)
        ttab = lbc_del (ttab  , nperio=nperio, cd_type='T', psgn=psgn)
    else : 
        ttab = lbc_del (utab_0, nperio=nperio, cd_type='T', psgn=psgn)
        
    if mmath == xr :
        if ax :
            ttab = ttab.assign_coords({ax:np.arange (ttab.shape[ix])+1.})
        if zdim and az :
            if az != zdim : ttab = ttab.rename( {az:zdim}) 
    return ttab

def V2T (vtab, nperio=None, psgn=-1.0, zdim=None, action='ave') :
    '''Interpolate an array from V grid to T grid (j-mean)'''
    mmath = __mmath__ (vtab)
    lperio, aperio = lbc_diag (nperio)
    vtab_0 = mmath.where ( np.isnan(vtab), 0., vtab)
    vtab_0 = lbc_add (vtab_0, nperio=nperio, cd_type='V', psgn=psgn)
    jy, ay = __findAxis__ (vtab_0, 'y')
    kz, az = __findAxis__ (vtab_0, 'z')
    if ay : 
        if action == 'ave'  : ttab = 0.5 *      (vtab_0 + np.roll (vtab_0, axis=jy, shift=1))
        if action == 'min'  : ttab = np.minimum (vtab_0 , np.roll (vtab_0, axis=jy, shift=1))
        if action == 'max'  : ttab = np.maximum (vtab_0 , np.roll (vtab_0, axis=jy, shift=1))
        if action == 'mult' : ttab =             vtab_0 * np.roll (vtab_0, axis=jy, shift=1)
        ttab = lbc_del (ttab  , nperio=nperio, cd_type='T', psgn=psgn)
    else :
        ttab = lbc_del (vtab_0, nperio=nperio, cd_type='T', psgn=psgn)
        
    if mmath == xr :
        if ay :
            ttab = ttab.assign_coords({ay:np.arange(ttab.shape[jy])+1.})
        if zdim and az :
            if az != zdim : ttab = ttab.rename( {az:zdim}) 
    return ttab

def F2T (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') :
    '''Interpolate an array from F grid to T grid (i- and j- means)'''
    mmath = __mmath__ (ftab)
    ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab)
    ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
    ttab = V2T (F2V (ftab_0, nperio=nperio, psgn=psgn, zdim=zdim, action=action), nperio=nperio, psgn=psgn, zdim=zdim, action=action)
    return lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn)

def T2U (ttab, nperio=None, psgn=1.0, zdim=None, action='ave') :
    '''Interpolate an array from T grid to U grid (i-mean)'''
    mmath = __mmath__ (ttab)
    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
    ttab_0 = lbc_add (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn)
    ix, ax = __findAxis__ (ttab_0, 'x')
    kz, az = __findAxis__ (ttab_0, 'z')
    if ix : 
        if action == 'ave'  : utab = 0.5 *      (ttab_0 + np.roll (ttab_0, axis=ix, shift=-1))
        if action == 'min'  : utab = np.minimum (ttab_0 , np.roll (ttab_0, axis=ix, shift=-1))
        if action == 'max'  : utab = np.maximum (ttab_0 , np.roll (ttab_0, axis=ix, shift=-1))
        if action == 'mult' : utab =             ttab_0 * np.roll (ttab_0, axis=ix, shift=-1)
        utab = lbc_del (utab  , nperio=nperio, cd_type='U', psgn=psgn)
    else :
        utab = lbc_del (ttab_0, nperio=nperio, cd_type='U', psgn=psgn)
        
    if mmath == xr :    
        if ax : 
            utab = ttab.assign_coords({ax:np.arange(utab.shape[ix])+1.})
        if zdim and az :
            if az != zdim : utab = utab.rename( {az:zdim}) 
    return utab

def T2V (ttab, nperio=None, psgn=1.0, zdim=None, action='ave') :
    '''Interpolate an array from T grid to V grid (j-mean)'''
    mmath = __mmath__ (ttab)
    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
    ttab_0 = lbc_add (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn)
    jy, ay = __findAxis__ (ttab_0, 'y')
    kz, az = __findAxis__ (ttab_0, 'z')
    if jy : 
        if action == 'ave'  : vtab = 0.5 *      (ttab_0 + np.roll (ttab_0, axis=jy, shift=-1))
        if action == 'min'  : vtab = np.minimum (ttab_0 , np.roll (ttab_0, axis=jy, shift=-1))
        if action == 'max'  : vtab = np.maximum (ttab_0 , np.roll (ttab_0, axis=jy, shift=-1))
        if action == 'mult' : vtab =             ttab_0 * np.roll (ttab_0, axis=jy, shift=-1)
        vtab = lbc_del (vtab  , nperio=nperio, cd_type='V', psgn=psgn)
    else :
        vtab = lbc_del (ttab_0, nperio=nperio, cd_type='V', psgn=psgn)

    if mmath == xr :
        if ay : 
            vtab = vtab.assign_coords({ay:np.arange(vtab.shape[jy])+1.})
        if zdim and az :
            if az != zdim : vtab = vtab.rename( {az:zdim}) 
    return vtab

def V2F (vtab, nperio=None, psgn=-1.0, zdim=None, action='ave') :
    '''Interpolate an array from V grid to F grid (i-mean)'''
    mmath = __mmath__ (vtab)
    vtab_0 = mmath.where ( np.isnan(vtab), 0., vtab)
    vtab_0 = lbc_add (vtab_0 , nperio=nperio, cd_type='V', psgn=psgn)
    ix, ax = __findAxis__ (vtab_0, 'x')
    kz, az = __findAxis__ (vtab_0, 'z')
    if ix : 
        if action == 'ave'  : 0.5 *      (vtab_0 + np.roll (vtab_0, axis=ix, shift=-1))
        if action == 'min'  : np.minimum (vtab_0 , np.roll (vtab_0, axis=ix, shift=-1))
        if action == 'max'  : np.maximum (vtab_0 , np.roll (vtab_0, axis=ix, shift=-1))
        if action == 'mult' :             vtab_0 * np.roll (vtab_0, axis=ix, shift=-1)
        ftab = lbc_del (ftab  , nperio=nperio, cd_type='F', psgn=psgn)
    else :
         ftab = lbc_del (vtab_0, nperio=nperio, cd_type='F', psgn=psgn)
   
    if mmath == xr :
        if ax : 
            ftab = ftab.assign_coords({ax:np.arange(ftab.shape[ix])+1.})
        if zdim and az :
            if az != zdim : ftab = ftab.rename( {az:zdim}) 
    return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)

def U2F (utab, nperio=None, psgn=-1.0, zdim=None, action='ave') :
    '''Interpolate an array from U grid to F grid i-mean)'''
    mmath = __mmath__ (utab)
    utab_0 = mmath.where ( np.isnan(utab), 0., utab)
    utab_0 = lbc_add (utab_0 , nperio=nperio, cd_type='U', psgn=psgn)
    jy, ay = __findAxis__ (utab_0, 'y')
    kz, az = __findAxis__ (utab_0, 'z')
    if jy : 
        if action == 'ave'  :    ftab = 0.5 *      (utab_0 + np.roll (utab_0, axis=jy, shift=-1))
        if action == 'min'  :    ftab = np.minimum (utab_0 , np.roll (utab_0, axis=jy, shift=-1))
        if action == 'max'  :    ftab = np.maximum (utab_0 , np.roll (utab_0, axis=jy, shift=-1))
        if action == 'mult' :    ftab =             utab_0 * np.roll (utab_0, axis=jy, shift=-1)
        ftab = lbc_del (ftab  , nperio=nperio, cd_type='F', psgn=psgn)
    else :
        ftab = lbc_del (utab_0, nperio=nperio, cd_type='F', psgn=psgn)
  
    if mmath == xr :
        if ay : 
            ftab = ftab.assign_coords({'y':np.arange(ftab.shape[jy])+1.})
        if zdim and az :
            if az != zdim : ftab = ftab.rename( {az:zdim}) 
    return ftab

def F2T (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') :
    '''Interpolate an array on F grid to T grid (i- and j- means)'''
    mmath = __mmath__ (ftab)
    ftab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
    ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
    ttab = U2T(F2U(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim, action=action), nperio=nperio, psgn=psgn, zdim=zdim, action=action)
    return lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn)

def T2F (ttab, nperio=None, psgn=1.0, zdim=None, action='mean') :
    '''Interpolate an array on T grid to F grid (i- and j- means)'''
    mmath = __mmath__ (ttab)
    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
    ttab_0 = lbc_add (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn)
    ftab = T2U (U2F (ttab, nperio=nperio, psgn=psgn, zdim=zdim, action=action), nperio=nperio, psgn=psgn, zdim=zdim, action=action)
    
    return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)

def F2U (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') :
    '''Interpolate an array on F grid to FUgrid (i-mean)'''
    mmath = __mmath__ (ftab)
    ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab)
    ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
    jy, ay = __findAxis__ (ftab_0, 'y')
    kz, az = __findAxis__ (ftab_0, 'z')
    if jy : 
        if action == 'ave'  : utab = 0.5 *      (ftab_0 + np.roll (ftab_0, axis=jy, shift=-1))
        if action == 'min'  : utab = np.minimum (ftab_0 , np.roll (ftab_0, axis=jy, shift=-1))
        if action == 'max'  : utab = np.maximum (ftab_0 , np.roll (ftab_0, axis=jy, shift=-1))
        if action == 'mult' : utab =             ftab_0 * np.roll (ftab_0, axis=jy, shift=-1)
        utab = lbc_del (utab  , nperio=nperio, cd_type='U', psgn=psgn)
    else :
        utab = lbc_del (ftab_0, nperio=nperio, cd_type='U', psgn=psgn)

    if mmath == xr :
        utab = utab.assign_coords({ay:np.arange(ftab.shape[jy])+1.})
        if zdim and zz :
            if az != zdim : utab = utab.rename( {az:zdim}) 
    return utab

def F2V (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') :
    '''Interpolate an array from F grid to V grid (i-mean)'''
    mmath = __mmath__ (ftab)
    ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab)
    ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
    ix, ax = __findAxis__ (ftab_0, 'x')
    kz, az = __findAxis__ (ftab_0, 'z')
    if ix : 
        if action == 'ave'  : vtab = 0.5 *      (ftab_0 + np.roll (ftab_0, axis=ix, shift=-1))
        if action == 'min'  : vtab = np.minimum (ftab_0 , np.roll (ftab_0, axis=ix, shift=-1))
        if action == 'max'  : vtab = np.maximum (ftab_0 , np.roll (ftab_0, axis=ix, shift=-1))
        if action == 'mult' : vtab =             ftab_0 * np.roll (ftab_0, axis=ix, shift=-1)
        vtab = lbc_del (vtab  , nperio=nperio, cd_type='V', psgn=psgn)
    else : 
        vtab = lbc_del (ftab_0, nperio=nperio, cd_type='V', psgn=psgn)

    if mmath == xr :
        vtab = vtab.assign_coords({ax:np.arange(ftab.shape[ix])+1.})
        if zdim and az :
            if az != zdim : vtab = vtab.rename( {az:zdim}) 
    return vtab

def W2T (wtab, zcoord=None, zdim=None, sval=np.nan) :
    '''
    Interpolate an array on W grid to T grid (k-mean)
    sval is the bottom value
    '''
    mmath = __mmath__ (wtab)
    wtab_0 = mmath.where ( np.isnan(wtab), 0., wtab)

    kz, az = __findAxis__ (wtab_0, 'z')

    if kz : 
        ttab = 0.5 * ( wtab_0 + np.roll (wtab_0, axis=kz, shift=-1) )
    else :
        ttab = wtab_0

    if mmath == xr :
        ttab[{az:kz}] = sval
        if zdim and az :
            if az != zdim : ttab = ttab.rename ( {az:zdim} )
        if zcoord is not None :
            ttab = ttab.assign_coords ( {zdim:zcoord} )
    else :
        ttab[..., -1, :, :] = sval

    return ttab

def T2W (ttab, zcoord=None, zdim=None, sval=np.nan, extrap_surf=False) :
    '''Interpolate an array from T grid to W grid (k-mean)
    sval is the surface value
    if extrap_surf==True, surface value is taken from 1st level value.
    '''
    mmath = __mmath__ (ttab)
    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
    kz, az = __findAxis__ (ttab_0, 'z')
    wtab = 0.5 * ( ttab_0 + np.roll (ttab_0, axis=kz, shift=1) )

    if mmath == xr :
        if extrap_surf : wtab[{az:0}] = ttabb[{az:0}]
        else           : wtab[{az:0}] = sval
    else : 
        if extrap_surf : wtab[..., 0, :, :] = ttab[..., 0, :, :]
        else           : wtab[..., 0, :, :] = sval

    if mmath == xr :
        if zdim and az :
            if az != zdim : wtab = wtab.rename ( {az:zdim})
        if zcoord is not None :
            wtab = wtab.assign_coords ( {zdim:zcoord})
        else :
            ztab = wtab.assign_coords ( {zdim:np.arange(ttab.shape[kz])+1.} )
    return wtab

def fill (ptab, nperio, cd_type='T', npass=1, sval=0.) :
    '''
    Fill sval values with mean of neighbours
   
    Inputs :
       ptab : input field to fill
       nperio, cd_type : periodicity characteristics
    '''       

    mmath = __mmath__ (ptab)

    DoPerio = False ; lperio = nperio
    if nperio == 4.2 :
        DoPerio = True ; lperio = 4
    if nperio == 6.2 :
        DoPerio = True ; lperio = 6
        
    if DoPerio :
        ztab = lbc_add (ptab, nperio=nperio, sval=sval)
    else : 
        ztab = ptab
        
    if np.isnan (sval) : 
        ztab   = mmath.where (np.isnan(ztab), np.nan, ztab)
    else :
        ztab   = mmath.where (ztab==sval    , np.nan, ztab)
   
    for nn in np.arange (npass) : 
        zmask = mmath.where ( np.isnan(ztab), 0., 1.   )
        ztab0 = mmath.where ( np.isnan(ztab), 0., ztab )
        # Compte du nombre de voisins
        zcount = 1./6. * ( zmask \
          + np.roll(zmask, shift=1, axis=-1) + np.roll(zmask, shift=-1, axis=-1) \
          + np.roll(zmask, shift=1, axis=-2) + np.roll(zmask, shift=-1, axis=-2) \
          + 0.5 * ( \
                + np.roll(np.roll(zmask, shift= 1, axis=-2), shift= 1, axis=-1) \
                + np.roll(np.roll(zmask, shift=-1, axis=-2), shift= 1, axis=-1) \
                + np.roll(np.roll(zmask, shift= 1, axis=-2), shift=-1, axis=-1) \
                + np.roll(np.roll(zmask, shift=-1, axis=-2), shift=-1, axis=-1) ) )

        znew =1./6. * ( ztab0 \
           + np.roll(ztab0, shift=1, axis=-1) + np.roll(ztab0, shift=-1, axis=-1) \
           + np.roll(ztab0, shift=1, axis=-2) + np.roll(ztab0, shift=-1, axis=-2) \
           + 0.5 * ( \
                + np.roll(np.roll(ztab0 , shift= 1, axis=-2), shift= 1, axis=-1) \
                + np.roll(np.roll(ztab0 , shift=-1, axis=-2), shift= 1, axis=-1) \
                + np.roll(np.roll(ztab0 , shift= 1, axis=-2), shift=-1, axis=-1) \
                + np.roll(np.roll(ztab0 , shift=-1, axis=-2), shift=-1, axis=-1) ) )

        zcount = lbc (zcount, nperio=lperio, cd_type=cd_type)
        znew   = lbc (znew  , nperio=lperio, cd_type=cd_type)
        
        ztab = mmath.where (np.logical_and (zmask==0., zcount>0), znew/zcount, ztab)

    ztab = mmath.where (zcount==0, sval, ztab)
    if DoPerio : ztab = lbc_del (ztab, nperio=lperio)

    return ztab

def correct_uv (u, v, lat) :
    '''
    Correct a Cartopy bug in orthographic projection

    See https://github.com/SciTools/cartopy/issues/1179

    The correction is needed with cartopy <= 0.20
    It seems that version 0.21 will correct the bug (https://github.com/SciTools/cartopy/pull/1926)
    Later note : the bug is still present in Cartopy 0.22

    Inputs :
       u, v : eastward/northward components
       lat  : latitude of the point (degrees north)

    Outputs : 
       modified eastward/nothward components to have correct polar projections in cartopy
    '''
    uv = np.sqrt (u*u + v*v)           # Original modulus
    zu = u
    zv = v * np.cos (rad*lat)
    zz = np.sqrt ( zu*zu + zv*zv )     # Corrected modulus
    uc = zu*uv/zz ; vc = zv*uv/zz      # Final corrected values
    return uc, vc

def norm_uv (u, v) :
    '''
    Return norm of a 2 components vector
    '''
    return np.sqrt (u*u + v*v)

def normalize_uv (u, v) :
    '''
    Normalize 2 components vector
    '''
    uv = norm_uv (u, v)
    return u/uv, v/uv

def msf (vv, e1v_e3v, lat1d, depthw) :
    '''
    Computes the meridonal stream function
    First input is meridional_velocity*e1v*e3v
    '''
 
    v_e1v_e3v = vv * e1v_e3v
    v_e1v_e3v.attrs = vv.attrs
    
    ix, ax = __findAxis__ (v_e1v_e3v, 'x')
    kz, az = __findAxis__ (v_e1v_e3v, 'z')
    if az == 'olevel' : new_az = 'olevel'
    else              : new_az = 'depthw'

    zomsf = -v_e1v_e3v.cumsum ( dim=az, keep_attrs=True).sum (dim=ax, keep_attrs=True)*1.E-6
    zomsf = zomsf - zomsf.isel ( { az:-1} )
    
    jy, ay = __findAxis__ (zomsf, 'y' )
    zomsf = zomsf.assign_coords ( {az:depthw.values, ay:lat1d.values})
    
    zomsf = zomsf.rename ( {ay:'lat'})
    if az != new_az : zomsf = zomsf.rename ( {az:new_az} )
    zomsf.attrs['standard_name'] = 'Meridional stream function'
    zomsf.attrs['long_name'] = 'Meridional stream function'
    zomsf.attrs['units'] = 'Sv'
    zomsf[new_az].attrs  = depthw.attrs
    zomsf.lat.attrs=lat1d.attrs
        
    return zomsf

def bsf (uu, e2u_e3u, mask, nperio=None, bsf0=None ) :
    '''
    Computes the barotropic stream function
    First input is zonal_velocity*e2u*e3u
    bsf0 is the point with bsf=0 
    (ex: bsf0={'x':3, 'y':120} for orca2, 
         bsf0={'x':5, 'y':300} for eeORCA1
    '''
    u_e2u_e3u       = uu * e2u_e3u
    u_e2u_e3u.attrs = uu.attrs

    iy, ay = __findAxis__ (u_e2u_e3u, 'y')
    kz, az = __findAxis__ (u_e2u_e3u, 'z')
    
    bsf = -u_e2u_e3u.cumsum ( dim=ay, keep_attrs=True )
    bsf = bsf.sum (dim=az, keep_attrs=True)*1.E-6
        
    if bsf0 :
        bsf = bsf - bsf.isel (bsf0)
       
    bsf = bsf.where (mask !=0, np.nan)
    for attr in uu.attrs :
        bsf.attrs[attr] = uu.attrs[attr]
    bsf.attrs['standard_name'] = 'barotropic_stream_function'
    bsf.attrs['long_name']     = 'Barotropic stream function'
    bsf.attrs['units']         = 'Sv'
    bsf = lbc (bsf, nperio=nperio, cd_type='F')
       
    return bsf

def namelist_read (ref=None, cfg=None, out='dict', flat=False, verbose=False) :
    '''
    Read NEMO namelist(s) and return either a dictionnary or an xarray dataset

    ref : file with reference namelist, or a f90nml.namelist.Namelist object
    cfg : file with config namelist, or a f90nml.namelist.Namelist object
    At least one namelist neaded

    out: 
        'dict' to return a dictonnary
        'xr'   to return an xarray dataset
    flat : only for dict output. Output a flat dictionnary with all values.
    
    '''
    
    import f90nml
    if ref :
        if isinstance (ref, str) : nml_ref = f90nml.read (ref)
        if isinstance (ref, f90nml.namelist.Namelist) : nml_ref = ref
        
    if cfg :
        if isinstance (cfg, str) : nml_cfg = f90nml.read (cfg)
        if isinstance (cfg, f90nml.namelist.Namelist) : nml_cfg = cfg
    
    if out == 'dict' : dict_namelist = {}
    if out == 'xr'   : xr_namelist = xr.Dataset ()

    list_nml = [] ; list_comment = []

    if ref : list_nml.append (nml_ref) ; list_comment.append ('ref')
    if cfg : list_nml.append (nml_cfg) ; list_comment.append ('cfg')

    for nml, comment in zip (list_nml, list_comment) :
        if verbose : print (comment)
        if flat and out =='dict' :
            for nam in nml.keys () :
                if verbose : print (nam)
                for value in nml[nam] :
                     if out == 'dict' : dict_namelist[value] = nml[nam][value]
                     if verbose : print (nam, ':', value, ':', nml[nam][value])
        else :
            for nam in nml.keys () :
                if verbose : print (nam)
                if out == 'dict' :
                    if nam not in dict_namelist.keys () : dict_namelist[nam] = {}
                for value in nml[nam] :
                    if out == 'dict' : dict_namelist[nam][value] = nml[nam][value]
                    if out == 'xr'   : xr_namelist[value] = nml[nam][value]
                    if verbose : print (nam, ':', value, ':', nml[nam][value])

    if out == 'dict' : return dict_namelist
    if out == 'xr'   : return xr_namelist

def fill_closed_seas (imask, nperio=None,  cd_type='T') :
    '''Fill closed seas with image processing library
    imask : mask, 1 on ocean, 0 on land
    '''
    from scipy import ndimage

    imask_filled = ndimage.binary_fill_holes ( lbc (imask, nperio=nperio, cd_type=cd_type))
    imask_filled = lbc ( imask_filled, nperio=nperio, cd_type=cd_type)

    return imask_filled

'''
Sea water state function parameters from NEMO code
'''
rdeltaS = 32. ; r1_S0  = 0.875/35.16504 ; r1_T0  = 1./40. ; r1_Z0  = 1.e-4

EOS000 =  8.0189615746e+02 ; EOS100 =  8.6672408165e+02 ; EOS200 = -1.7864682637e+03 ; EOS300 =  2.0375295546e+03 ; EOS400 = -1.2849161071e+03 ; EOS500 =  4.3227585684e+02 ; EOS600 = -6.0579916612e+01
EOS010 =  2.6010145068e+01 ; EOS110 = -6.5281885265e+01 ; EOS210 =  8.1770425108e+01 ; EOS310 = -5.6888046321e+01 ; EOS410 =  1.7681814114e+01 ; EOS510 = -1.9193502195
EOS020 = -3.7074170417e+01 ; EOS120 =  6.1548258127e+01 ; EOS220 = -6.0362551501e+01 ; EOS320 =  2.9130021253e+01 ; EOS420 = -5.4723692739     ; EOS030 =  2.1661789529e+01 
EOS130 = -3.3449108469e+01 ; EOS230 =  1.9717078466e+01 ; EOS330 = -3.1742946532
EOS040 = -8.3627885467     ; EOS140 =  1.1311538584e+01 ; EOS240 = -5.3563304045
EOS050 =  5.4048723791e-01 ; EOS150 =  4.8169980163e-01
EOS060 = -1.9083568888e-01
EOS001 =  1.9681925209e+01 ; EOS101 = -4.2549998214e+01 ; EOS201 =  5.0774768218e+01 ; EOS301 = -3.0938076334e+01 ; EOS401 =   6.6051753097    ; EOS011 = -1.3336301113e+01
EOS111 = -4.4870114575     ; EOS211 =  5.0042598061     ; EOS311 = -6.5399043664e-01 ; EOS021 =  6.7080479603     ; EOS121 =   3.5063081279
EOS221 = -1.8795372996     ; EOS031 = -2.4649669534     ; EOS131 = -5.5077101279e-01 ; EOS041 =  5.5927935970e-01
EOS002 =  2.0660924175     ; EOS102 = -4.9527603989     ; EOS202 =  2.5019633244     ; EOS012 =  2.0564311499     ; EOS112 = -2.1311365518e-01 ; EOS022 = -1.2419983026
EOS003 = -2.3342758797e-02 ; EOS103 = -1.8507636718e-02 ; EOS013 =  3.7969820455e-01 

def rhop ( ptemp, psal ) :
    '''
    Potential density referenced to surface
    Computation from NEMO code
    '''
    zt  = ptemp * r1_T0                                  # Temperature (°C)
    zs  = np.sqrt ( np.abs( psal + rdeltaS ) * r1_S0 )   # Square root of salinity (PSS)
    #
    prhop = (((((EOS060*zt   \
             + EOS150*zs     + EOS050)*zt   \
             + (EOS240*zs    + EOS140)*zs + EOS040)*zt   \
             + ((EOS330*zs   + EOS230)*zs + EOS130)*zs + EOS030)*zt   \
             + (((EOS420*zs  + EOS320)*zs + EOS220)*zs + EOS120)*zs + EOS020)*zt   \
             + ((((EOS510*zs + EOS410)*zs + EOS310)*zs + EOS210)*zs + EOS110)*zs + EOS010)*zt   \
             + (((((EOS600*zs+ EOS500)*zs + EOS400)*zs + EOS300)*zs + EOS200)*zs + EOS100)*zs + EOS000
    #
    return prhop

def rho ( pdep, ptemp, psal ) :
    '''
    In situ density
    Computation from NEMO code
    '''
    zh  = pdep  * r1_Z0                                  # Depth (m)
    zt  = ptemp * r1_T0                                  # Temperature (°C)
    zs  = np.sqrt ( np.abs( psal + rdeltaS ) * r1_S0 )   # Square root salinity (PSS)
    #
    zn3 = EOS013*zt + EOS103*zs+EOS003
    #
    zn2 = (EOS022*zt + EOS112*zs+EOS012)*zt + (EOS202*zs+EOS102)*zs+EOS002
    #
    zn1 = (((EOS041*zt   \
         + EOS131*zs   + EOS031)*zt   \
         + (EOS221*zs  + EOS121)*zs + EOS021)*zt   \
        + ((EOS311*zs  + EOS211)*zs + EOS111)*zs + EOS011)*zt   \
        + (((EOS401*zs + EOS301)*zs + EOS201)*zs + EOS101)*zs + EOS001
    #
    zn0 = (((((EOS060*zt   \
             + EOS150*zs      + EOS050)*zt   \
             + (EOS240*zs     + EOS140)*zs + EOS040)*zt   \
             + ((EOS330*zs    + EOS230)*zs + EOS130)*zs + EOS030)*zt   \
             + (((EOS420*zs   + EOS320)*zs + EOS220)*zs + EOS120)*zs + EOS020)*zt   \
             + ((((EOS510*zs  + EOS410)*zs + EOS310)*zs + EOS210)*zs + EOS110)*zs + EOS010)*zt   \
             + (((((EOS600*zs + EOS500)*zs + EOS400)*zs + EOS300)*zs + EOS200)*zs + EOS100)*zs + EOS000
    #
    prho  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
    #
    return prho

## ===========================================================================
##
##                               That's all folk's !!!
##
## ===========================================================================

def __is_orca_north_fold__ ( Xtest, cname_long='T' ) :
    '''
    Ported (pirated !!?) from Sosie

    Tell if there is a 2/point band overlaping folding at the north pole typical of the ORCA grid

    0 => not an orca grid (or unknown one)
    4 => North fold T-point pivot (ex: ORCA2)
    6 => North fold F-point pivot (ex: ORCA1)

    We need all this 'cname_long' stuff because with our method, there is a
    confusion between "Grid_U with T-fold" and "Grid_V with F-fold"
    => so knowing the name of the longitude array (as in namelist, and hence as
    in netcdf file) might help taking the righ decision !!! UGLY!!!
    => not implemented yet
    '''
    
    ifld_nord =  0 ; cgrd_type = 'X'
    ny, nx = Xtest.shape[-2:]

    if ny > 3 : # (case if called with a 1D array, ignoring...)
        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. :
          ifld_nord = 4 ; cgrd_type = 'T' # T-pivot, grid_T      

        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2  :-1] ).sum() == 0. :
            if cnlon == 'U' : ifld_nord = 4 ;  cgrd_type = 'U' # T-pivot, grid_T
                ## LOLO: PROBLEM == 6, V !!!

        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. :
            ifld_nord = 4 ; cgrd_type = 'V' # T-pivot, grid_V

        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-2, nx-1-1:nx-nx//2:-1] ).sum() == 0. :
            ifld_nord = 6 ; cgrd_type = 'T'# F-pivot, grid_T

        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-1, nx-1:nx-nx//2-1:-1] ).sum() == 0. :
            ifld_nord = 6 ;  cgrd_type = 'U' # F-pivot, grid_U

        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2  :-1] ).sum() == 0. :
            if cnlon == 'V' : ifld_nord = 6 ; cgrd_type = 'V' # F-pivot, grid_V
                ## LOLO: PROBLEM == 4, U !!!

    return ifld_nord, cgrd_type