Index: TOOLS/MOSAIX/nemo.py
===================================================================
--- TOOLS/MOSAIX/nemo.py	(revision 6244)
+++ TOOLS/MOSAIX/nemo.py	(revision 6245)
@@ -2,23 +2,18 @@
 ## ===========================================================================
 ##
-##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt"
-##  file for an english version of the licence and
-##  "Licence_CeCILL_V2-fr.txt" for a french version.
+##  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".
 ##
-##  Permission is hereby granted, free of charge, to any person or
-##  organization obtaining a copy of the software and accompanying
-##  documentation covered by this license (the "Software") to use,
-##  reproduce, display, distribute, execute, and transmit the
-##  Software, and to prepare derivative works of the Software, and to
-##  permit third-parties to whom the Software is furnished to do so,
-##  all subject to the following:
-##
-##  Warning, to install, configure, run, use any of MOSAIX 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).  Authors 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
+##  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.
 ##
 ## ===========================================================================
@@ -37,6 +32,15 @@
 __HeadURL    = "$HeadURL$"
 
-import sys, numpy as np
+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
+
+try    : import numba
 except : pass
 
@@ -62,26 +66,21 @@
 tList = [ 't', 'T', 'time'  , ]
 
-## SVN information
-__Author__   = "$Author$"
-__Date__     = "$Date$"
-__Revision__ = "$Revision$"
-__Id__       = "$Id$"
-__HeadURL    = "$HeadURL$"
-
-def __math__ (tab, default=None) :
-    math = default
+## ===========================================================================
+def __mmath__ (tab, default=None) :
+    mmath = default
     try    :
-        if type (tab) == xr.core.dataarray.DataArray : math = xr
+        if type (tab) == xr.core.dataarray.DataArray : mmath = xr
     except :
         pass
 
     try    :
-        if type (tab) == np.ndarray : math = np
+        if type (tab) == np.ndarray : mmath = np
     except :
         pass
             
-    return math
-
-def __guessNperio__ (jpj, jpi, nperio=None) :
+    return mmath
+
+
+def __guessNperio__ (jpj, jpi, nperio=None, out='nperio') :
     '''
     Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details)
@@ -93,48 +92,63 @@
     '''
     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
+    '''
+    if nperio == None :
         ## Values for NEMO version < 4.2
-        if jpi ==  182 :
+        if jpj ==  149 and jpi ==  182 :
+            config = 'ORCA2.3'
             nperio = 4   # ORCA2. We choose legacy orca2.
             Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'T'
-        if jpi ==  362 : # eORCA1.
+        if jpj ==  332 and jpi ==  362 : # 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 ==  149 : # ORCA2. We choose legacy orca2.
-            nperio = 4
-            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
         if jpj ==  294 : # ORCA1
+            config = 'ORCA1'
             nperio = 6
-            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
-        if jpj ==  332 : # eORCA1.
-            nperio = 6  
             Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
             
         ## Values for NEMO version >= 4.2. No more halo points
-        if jpi ==  180 :
+        if jpj == 148 and jpi ==  180 :
+            config = 'ORCA2.4'
             nperio = 4.2 # ORCA2. We choose legacy orca2.
             Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F'
-        if jpi ==  360 : # eORCA1.
+        if jpj == 331  and jpi ==  360 : # 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 jpj ==  148 : nperio = 4.2 # ORCA2. We choose legacy orca2.
-        if jpj ==  330 : nperio = 6.2 # eORCA1.
-
+            
         if nperio == None :
-            sys.exit ('in nemo module : nperio not found, and cannot by guessed')
+            raise Exception  ('in nemo module : nperio not found, and cannot by guessed')
         else :
             if nperio in nperio_valid_range :
-                print ('nperio set as {:} (deduced from jpi={:d})'.format (nperio, jpi))
+                print ('nperio set as {:} (deduced from jpj={:d} jpi={:d})'.format (nperio, jpj, jpi))
             else : 
-                sys.exit  ('nperio set as {:} (deduced from jpi={:d}) : nemo.py is not ready for this value'.format (nperio, jpi))
-            
-    return nperio
-
+                raise ValueError ('nperio set as {:} (deduced from jpi={:d}) : nemo.py is not ready for this value'.format (nperio, jpi))
+
+    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) :
     '''
@@ -149,6 +163,6 @@
     '''
     gP = None
-    math = __math__ (ptab)
-    if math == xr :
+    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'
@@ -162,11 +176,19 @@
              
         if gP == None :
-            sys.exit ('in nemo module : cd_type not found, and cannot by guessed')
+            raise Exception ('in nemo module : cd_type not found, and cannot by guessed')
         else :
             print ('Grid set as', gP, 'deduced from dims ', ptab.dims)
             return gP
     else :
-         sys.exit ('in nemo module : cd_type not found, input is not an xarray data')
-         #return gP
+         raise Exception  ('in nemo module : cd_type not found, input is not an xarray data')
+
+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') :
@@ -174,5 +196,5 @@
     Find number and name of the requested axis
     '''
-    math = __math__ (tab)
+    mmath = __mmath__ (tab)
     ix = None ; ax = None
 
@@ -198,5 +220,5 @@
         unList = [ 'second', 'minute', 'hour', 'day', 'month' ]
     
-    if math == xr :
+    if mmath == xr :
         for Name in axList :
             try    :
@@ -222,5 +244,6 @@
        
     return ix, ax
-                
+
+#@numba.jit(forceobj=True)
 def fixed_lon (lon, center_lon=0.0) :
     '''
@@ -232,13 +255,13 @@
     Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7
     '''
-    math = __math__ (lon)
+    mmath = __mmath__ (lon)
     
     fixed_lon = lon.copy ()
         
-    fixed_lon = math.where (fixed_lon > center_lon+180., fixed_lon-360.0, fixed_lon)
-    fixed_lon = math.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)
+    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
+        fixed_lon [..., i, start+1:] += 360.
 
     # Special case for eORCA025
@@ -246,8 +269,64 @@
     if fixed_lon.shape [-1] == 1440 : fixed_lon [..., -1, :] = fixed_lon [..., -2, :]
 
-    if fixed_lon.min () > center_lon : fixed_lon = fixed_lon-360.0
+    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
 
+#@numba.jit(forceobj=True)
+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, has no
+    lon/lat values
+    '''
+    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
+
+#@numba.jit(forceobj=True)
+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, as no
+    lon/lat values
+    '''
+    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
+
+#@numba.jit(forceobj=True)
 def jeq (lat) :
     '''
@@ -256,9 +335,9 @@
     lat : latitudes of the grid. At least 2D.
     '''
-    math = __math__ (lat)
+    mmath = __mmath__ (lat)
     ix, ax = __findAxis__ (lat, 'x')
     iy, ay = __findAxis__ (lat, 'y')
 
-    if math == xr :
+    if mmath == xr :
         jeq = int ( np.mean ( np.argmin (np.abs (np.float64 (lat)), axis=iy) ) )
     else : 
@@ -266,4 +345,5 @@
     return jeq
 
+#@numba.jit(forceobj=True)
 def lon1D (lon, lat=None) :
     '''
@@ -273,24 +353,30 @@
     lat (optionnal) : latitudes of the grid
     '''
+    mmath = __mmath__ (lon)
     if np.max (lat) != None :
         je    = jeq (lat)
-        lon1D = lon [..., je, :]
+        lon1D = lon.copy() [..., je, :]
     else :
         jpj, jpi = lon.shape [-2:]
-        lon1D    = lon [..., jpj//3, :]
+        lon1D    = lon.copy() [..., jpj//3, :]
 
     start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1)
     lon1D [..., start+1:] += 360
+
+    if mmath == xr :
+        lon1D.attrs = lon.attrs
+        lon1D = lon1D.assign_coords ( {'x':lon1D} )
         
     return lon1D
 
+#@numba.jit(forceobj=True)
 def latreg (lat, diff=0.1) :
     '''
-    Returns maximum j index where gridlines are along latitude in the northern hemisphere
+    Returns maximum j index where gridlines are along latitudes in the northern hemisphere
     
     lat : latitudes of the grid (2D)
     diff [optional] : tolerance
     '''
-    math = __math__ (lat)
+    mmath = __mmath__ (lat)
     if diff == None :
         dy   = np.float64 (np.mean (np.abs (lat - np.roll (lat,shift=1,axis=-2, roll_coords=False))))
@@ -304,4 +390,5 @@
     return jreg, latreg
 
+#@numba.jit(forceobj=True)
 def lat1D (lat) :
     '''
@@ -310,5 +397,5 @@
     lat : latitudes of the grid (2D)
     '''
-    math = __math__ (lat)
+    mmath = __mmath__ (lat)
     jpj, jpi = lat.shape[-2:]
 
@@ -326,10 +413,13 @@
         yrange = 90.    -lat_reg
         
-    lat1D = math.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1))   
-        
-    if math == xr : lat1D.attrs = lat.attrs
+    lat1D = mmath.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1))   
+        
+    if mmath == xr :
+        lat1D.attrs = lat.attrs
+        lat1D = lat1D.assign_coords ( {'y':lat1D} )
 
     return lat1D
 
+#@numba.jit(forceobj=True)
 def latlon1D (lat, lon) :
     '''
@@ -340,6 +430,7 @@
     return lat1D (lat),  lon1D (lon, lat)
 
+##@numba.jit(forceobj=True)
 def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) :
-    math = __math__ (ptab)
+    mmath = __mmath__ (ptab)
     try :
         lon = lon.copy().to_masked_array()
@@ -350,8 +441,9 @@
             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 = math.where (mask, ptab, np.nan)
+    tab = mmath.where (mask, ptab, np.nan)
     
     return tab
-        
+
+#@numba.jit(forceobj=True)      
 def extend (tab, Lon=False, jplus=25, jpi=None, nperio=4) :
     '''
@@ -366,5 +458,5 @@
     
     '''
-    math = __math__ (tab)
+    mmath = __mmath__ (tab)
     
     if tab.shape[-1] == 1 : extend = tab
@@ -386,5 +478,5 @@
                 istart = 1 ; le=jpi-2 ; la=1  # Perfect, except at the pole that should be masked by lbc_plot
            
-            if math == xr :
+            if mmath == xr :
                 extend = np.concatenate ((tab.values[..., istart   :istart+le+1    ] + xplus,
                                           tab.values[..., istart+la:istart+la+jplus]         ), axis=-1)
@@ -402,16 +494,17 @@
 def orca2reg (ff, lat_name='nav_lat', lon_name='nav_lon', y_name='y', x_name='x') :
     '''
-    Assign ORCA dataset on a regular grid.
+    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
+      
+      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 : 
@@ -450,10 +543,9 @@
     
     if nperio not in nperio_valid_range :
-        print ('nperio=', nperio, ' is not in the valid range', nperio_valid_range)
-        sys.exit ()
+        raise Exception ('nperio=', nperio, ' is not in the valid range', nperio_valid_range)
 
     return jpj, jpi, nperio
         
-        
+#@numba.jit(forceobj=True)
 def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) :
     '''
@@ -468,8 +560,9 @@
     jpj, jpi, nperio = lbc_init (ptab, nperio)
     psgn   = ptab.dtype.type (psgn)
-    math = __math__ (ptab)
-    
-    if math == xr : ztab = ptab.values.copy ()
-    else          : ztab = ptab.copy ()
+    mmath = __mmath__ (ptab)
+    
+    if mmath == xr : ztab = ptab.values.copy ()
+    else           : ztab = ptab.copy ()
+        
     #
     #> East-West boundary conditions
@@ -548,9 +641,11 @@
         ztab [..., :, -1] = ztab [..., :,  1]
 
-    if math == xr :
+    if mmath == xr :
         ztab = xr.DataArray ( ztab, dims=ptab.dims, coords=ptab.coords)
+        ztab.attrs = ptab.attrs
         
     return ztab
 
+#@numba.jit(forceobj=True)
 def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) :
     #
@@ -632,4 +727,5 @@
     return ztab
 
+#@numba.jit(forceobj=True)
 def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) :
     '''
@@ -707,8 +803,9 @@
     return ztab
 
-def lbc_add (ptab, nperio=None, cd_type =None) :
+#@numba.jit(forceobj=True)
+def lbc_add (ptab, nperio=None, cd_type=None, psgn=1, sval=None) :
     '''
     Handle NEMO domain changes between NEMO 4.0 to NEMO 4.2
-      Peridodicity halo has been removes
+      Peridodicity halo has been removed
     This routine adds the halos if needed
 
@@ -719,5 +816,5 @@
     See NEMO documentation for further details
     '''
-
+    mmath = __mmath__ (ptab) 
     jpj, jpi, nperio = lbc_init (ptab, nperio)
 
@@ -725,5 +822,4 @@
 
     if nperio == 4.2 or nperio == 6.2 :
-        math = __math__ (ptab) 
       
         ext_shape = t_shape
@@ -731,24 +827,27 @@
         ext_shape[-2] = ext_shape[-2] + 1
 
-        if math == xr : ptab_ext = xr.DataArray (np.empty (ext_shape), dims=ptab.dims) 
-        else          : ptab_ext =               np.empty (ext_shape)
+        if mmath == xr :
+            ptab_ext = xr.DataArray (np.zeros (ext_shape), dims=ptab.dims) 
+            ptab_ext.values[..., :-1, 1:-1] = ptab.values.copy ()
+        else           :
+            ptab_ext =               np.zeros (ext_shape)
+            ptab_ext[..., :-1, 1:-1] = ptab.copy ()
             
-        ptab_ext[..., :-1, 1:-1] = ptab
-        
-        if nperio == 4.2 : ptab_ext = lbc (ptab_ext, 4, cd_type)
-        if nperio == 6.2 : ptab_ext = lbc (ptab_ext, 6, cd_type)
-
-        if math == xr :
-            for attr in ptab.attrs :
-                ptab_ext.attrs [attr] = ptab.attrs [attr]
-
-    else : ptab_ext = ptab
-        
-    return ptab
-
-def lbc_del (ptab, nperio=None) :
+        #if sval != None :  ptab_ext[..., 0, :] = sval
+        
+        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
+
+    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) :
     '''
     Handle NEMO domain changes between NEMO 4.0 to NEMO 4.2
-      Peridodicity halo has been removes
+      Periodicity halo has been removed
     This routine removes the halos if needed
 
@@ -762,9 +861,10 @@
     jpj, jpi, nperio = lbc_init (ptab, nperio)
 
-    if nperio == 4 or nperio == 6 :
-        return ptab[..., :-1, 1:-1]
+    if nperio == 4.2 or nperio == 6.2 :
+        return lbc (ptab[..., :-1, 1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn)
     else :
         return ptab
 
+#@numba.jit(forceobj=True)
 def lbc_index (jj, ii, jpj, jpi, nperio=None, cd_type='T') :
     '''
@@ -784,6 +884,6 @@
     jy = jj + 1 ; ix = ii + 1
 
-    math = __math__ (jj)
-    if math == None : math=np
+    mmath = __mmath__ (jj)
+    if mmath == None : mmath=np
 
     #
@@ -792,11 +892,11 @@
     if nperio in [1, 4, 6] :
         #... cyclic
-        ix = math.where (ix==jpi, 2   , ix)
-        ix = math.where (ix== 1 ,jpi-1, ix)
+        ix = mmath.where (ix==jpi, 2   , ix)
+        ix = mmath.where (ix== 1 ,jpi-1, ix)
 
     #
     def modIJ (cond, jy_new, ix_new) :
-        jy_r = math.where (cond, jy_new, jy)
-        ix_r = math.where (cond, ix_new, ix)
+        jy_r = mmath.where (cond, jy_new, jy)
+        ix_r = mmath.where (cond, ix_new, ix)
         return jy_r, ix_r
     #
@@ -907,5 +1007,5 @@
     else                         : jpj = len(lat_grid) ; jpi=len(lon_grid)
 
-    math = __math__ (lat_grid)
+    mmath = __mmath__ (lat_grid)
         
     # Compute distance from the point to all grid points (in radian)
@@ -943,11 +1043,11 @@
 def clo_lon (lon, lon0) :
     '''Choose closest to lon0 longitude, adding or substacting 360Â° if needed'''
-    math = __math__ (lon, np)
+    mmath = __mmath__ (lon, np)
         
     clo_lon = lon
-    clo_lon = math.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon)
-    clo_lon = math.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon)
-    clo_lon = math.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon)
-    clo_lon = math.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon)
+    clo_lon = mmath.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon)
+    clo_lon = mmath.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon)
+    clo_lon = mmath.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon)
+    clo_lon = mmath.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon)
     if clo_lon.shape == () : clo_lon = clo_lon.item ()
     return clo_lon
@@ -955,31 +1055,40 @@
 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'''
-    math = __math__ (glamt)
+    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*glamt) * np.tan (rpi/4.0 - rad*gphit/2.0)
-    zynpt = 0. - 2.0 * np.sin (rad*glamt) * np.tan (rpi/4.0 - rad*gphit/2.0)
+    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*glamu) * np.tan (rpi/4.0 - rad*gphiu/2.0)
-    zynpu = 0. - 2.0 * np.sin (rad*glamu) * np.tan (rpi/4.0 - rad*gphiu/2.0)
+    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*glamv) * np.tan (rpi/4.0 - rad*gphiv/2.0)
-    zynpv = 0. - 2.0 * np.sin (rad*glamv) * np.tan (rpi/4.0 - rad*gphiv/2.0)
+    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*glamf ) * np.tan ( rpi/4. - rad*gphif/2. )
-    zynpf = 0. - 2.0 * np.sin( rad*glamf ) * np.tan ( rpi/4. - rad*gphif/2. )
+    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 = glamv  
-    zphi = gphiv
-    zlan = np.roll ( glamv, axis=-2, shift=1)  # glamv (ji,jj-1)
-    zphh = np.roll ( gphiv, axis=-2, shift=1)  # gphiv (ji,jj-1)
+    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. )
@@ -989,8 +1098,8 @@
 
     # j-direction: f-point segment direction (around u-point)
-    zlam = glamf
-    zphi = gphif
-    zlan = np.roll (glamf, axis=-2, shift=1) # glamf (ji,jj-1)
-    zphh = np.roll (gphif, axis=-2, shift=1) # gphif (ji,jj-1)
+    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. )
@@ -1000,8 +1109,8 @@
 
     # i-direction: f-point segment direction (around v-point)
-    zlam = glamf  
-    zphi = gphif
-    zlan = np.roll (glamf, axis=-1, shift=1) # glamf (ji-1,jj)
-    zphh = np.roll (gphif, axis=-1, shift=1) # gphif (ji-1,jj)
+    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. )
@@ -1011,8 +1120,8 @@
 
     # j-direction: u-point segment direction (around f-point)
-    zlam = np.roll (glamu, axis=-2, shift=-1) # glamu (ji,jj+1) 
-    zphi = np.roll (gphiu, axis=-2, shift=-1) # gphiu (ji,jj+1)
-    zlan = glamu
-    zphh = gphiu
+    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. )
@@ -1035,14 +1144,23 @@
     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.)
-
-    if math == xr :
+    #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 )
@@ -1058,15 +1176,19 @@
 def angle (glam, gphi, nperio, cd_type='T') :
     '''Compute sinus and cosinus of model line direction with respect to east'''
-    math = __math__ (glam) 
+    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*glam) * np.tan (rpi/4.0 - rad*gphi/2.0)
-    zynp = 0. - 2.0 * np.sin (rad*glam) * np.tan (rpi/4.0 - rad*gphi/2.0)
+    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 (glam, axis=-2, shift=-1) # glam [jj+1, ji]
-    zphh_n = np.roll (gphi, axis=-2, shift=-1) # gphi [jj+1, ji]
-    zlan_s = np.roll (glam, axis=-2, shift= 1) # glam [jj-1, ji]
-    zphh_s = np.roll (gphi, axis=-2, shift= 1) # gphi [jj-1, ji]
+    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) \
@@ -1079,8 +1201,8 @@
     gcos = (zxnp*zxff + zynp*zyff) / znff
 
-    gsin = lbc (gsin, cd_type=cd_type, nperio=nperio, psgn=-1.)
-    gcos = lbc (gcos, cd_type=cd_type, nperio=nperio, psgn=-1.)
-
-    if math == xr :
+    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 )
@@ -1124,5 +1246,5 @@
     east-north components on the T grid point   
     '''
-    math = __math__ (uo)
+    mmath = __mmath__ (uo)
 
     ut = U2T (uo, nperio=nperio, psgn=-1.0, zdim=zdim)
@@ -1139,10 +1261,10 @@
 def rot_uv2enF ( uo, vo, gsinf, gcosf, nperio, zdim='deptht' ) :
     '''
-    ** Purpose :   Rotate the Repere: Change vector componantes from
+    ** 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   
     '''
-    math = __math__ (uo)
+    mmath = __mmath__ (uo)
 
     uf = U2F (uo, nperio=nperio, psgn=-1.0, zdim=zdim)
@@ -1156,15 +1278,21 @@
     
     return u_e, v_n
-             
-def U2T (utab, nperio=None, psgn=-1.0, zdim='deptht') :
+
+#@numba.jit(forceobj=True)
+def U2T (utab, nperio=None, psgn=-1.0, zdim='deptht', action='ave') :
     '''Interpolate an array from U grid to T grid i-mean)'''
-    math = __math__ (utab)
-    utab_0 = math.where ( np.isnan(utab), 0., utab)
-    utab_0 = lbc (utab_0 , nperio=nperio, cd_type='U', psgn=psgn)
+    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')
     iz, az = __findAxis__ (utab_0, 'z')
-    ttab =  lbc ( 0.5 * (utab_0 + np.roll (utab_0, axis=ix, shift=1)), cd_type='T', nperio=nperio, psgn=psgn)
-    
-    if math == xr :
+    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)
+    
+    if mmath == xr :
         if ax != None :
             ttab = ttab.assign_coords({ax:np.arange (ttab.shape[ix])+1.})
@@ -1173,14 +1301,19 @@
     return ttab
 
-def V2T (vtab, nperio=None, psgn=-1.0, zdim='deptht') :
+#@numba.jit(forceobj=True)
+def V2T (vtab, nperio=None, psgn=-1.0, zdim='deptht', action='ave') :
     '''Interpolate an array from V grid to T grid (j-mean)'''
-    math = __math__ (vtab)
-    vtab_0 = math.where ( np.isnan(vtab), 0., vtab)
-    vtab_0 = lbc (vtab_0 , nperio=nperio, cd_type='V', psgn=psgn)
+    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)
     iy, ay = __findAxis__ (vtab_0, 'y')
     iz, az = __findAxis__ (vtab_0, 'z')
-    ttab = lbc ( 0.5 * (vtab_0 + np.roll (vtab_0, axis=iy, shift=1)), cd_type='T', nperio=nperio, psgn=psgn)
-                     
-    if math == xr :
+    if action == 'ave'  : ttab = 0.5 *      (vtab_0 + np.roll (vtab_0, axis=iy, shift=1))
+    if action == 'min'  : ttab = np.minimum (vtab_0 , np.roll (vtab_0, axis=iy, shift=1))
+    if action == 'max'  : ttab = np.maximum (vtab_0 , np.roll (vtab_0, axis=iy, shift=1))
+    if action == 'mult' : ttab =             vtab_0 * np.roll (vtab_0, axis=iy, shift=1)
+    ttab = lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn)
+    if mmath == xr :
         if ay !=None : 
             ttab = ttab.assign_coords({ay:np.arange(ttab.shape[iy])+1.})
@@ -1189,22 +1322,28 @@
     return ttab
 
-def F2T (ftab, nperio=None, psgn=1.0, zdim='depthf') :
+#@numba.jit(forceobj=True)
+def F2T (ftab, nperio=None, psgn=1.0, zdim='depthf', action='ave') :
     '''Interpolate an array from F grid to T grid (i- and j- means)'''
-    math = __math__ (ftab)
-    ftab_0 = math.where ( np.isnan(ftab), 0., ftab)
-    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
-    ttab = V2T(F2V(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim), nperio=nperio, psgn=psgn, zdim=zdim)
-    return ttab
-
-def T2U (ttab, nperio=None, psgn=1.0, zdim='depthu') :
+    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)
+
+#@numba.jit(forceobj=True)
+def T2U (ttab, nperio=None, psgn=1.0, zdim='depthu', action='ave') :
     '''Interpolate an array from T grid to U grid (i-mean)'''
-    math = __math__ (ttab)
-    ttab_0 = math.where ( np.isnan(ttab), 0., ttab)
-    ttab_0 = lbc (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn)
+    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')
     iz, az = __findAxis__ (ttab_0, 'z')
-    utab = lbc ( 0.5 * (ttab_0 + np.roll (ttab_0, axis=ix, shift=-1)), cd_type='U', nperio=nperio, psgn=psgn)
-
-    if math == xr :    
+    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)
+
+    if mmath == xr :    
         if ax != None : 
             utab = ttab.assign_coords({ax:np.arange(utab.shape[ix])+1.})
@@ -1213,14 +1352,19 @@
     return utab
 
-def T2V (ttab, nperio=None, psgn=1.0, zdim='depthv') :
+#@numba.jit(forceobj=True)
+def T2V (ttab, nperio=None, psgn=1.0, zdim='depthv', action='ave') :
     '''Interpolate an array from T grid to V grid (j-mean)'''
-    math = __math__ (ttab)
-    ttab_0 = math.where ( np.isnan(ttab), 0., ttab)
-    ttab_0 = lbc (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn)
+    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)
     iy, ay = __findAxis__ (ttab_0, 'y')
     iz, az = __findAxis__ (ttab_0, 'z')
-    vtab = lbc ( 0.5 * (ttab_0 + np.roll (ttab_0, axis=iy, shift=-1)), cd_type='V', nperio=nperio, psgn=psgn)
-
-    if math == xr :
+    if action == 'ave'  : vtab = 0.5 *      (ttab_0 + np.roll (ttab_0, axis=iy, shift=-1))
+    if action == 'min'  : vtab = np.minimum (ttab_0 , np.roll (ttab_0, axis=iy, shift=-1))
+    if action == 'max'  : vtab = np.maximum (ttab_0 , np.roll (ttab_0, axis=iy, shift=-1))
+    if action == 'mult' : vtab =             ttab_0 * np.roll (ttab_0, axis=iy, shift=-1)
+
+    vtab = lbc_del (vtab, nperio=nperio, cd_type='V', psgn=psgn)
+    if mmath == xr :
         if ay != None : 
             vtab = vtab.assign_coords({ay:np.arange(vtab.shape[iy])+1.})
@@ -1229,30 +1373,40 @@
     return vtab
 
-def V2F (vtab, nperio=None, psgn=-1.0, zdim='depthf') :
+#@numba.jit(forceobj=True)
+def V2F (vtab, nperio=None, psgn=-1.0, zdim='depthf', action='ave') :
     '''Interpolate an array from V grid to F grid (i-mean)'''
-    math = __math__ (vtab)
-    vtab_0 = math.where ( np.isnan(vtab), 0., vtab)
-    vtab_0 = lbc (vtab_0 , nperio=nperio, cd_type='V', psgn=psgn)
+    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')
     iz, az = __findAxis__ (vtab_0, 'z')
-    ftab = lbc ( 0.5 * (vtab_0 + np.roll (vtab_0, axis=ix, shift=-1)), cd_type='F', nperio=nperio, psgn=psgn)
-
-    if math == xr :
+    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)
+    
+    if mmath == xr :
         if ax != None : 
             ftab = ftab.assign_coords({ax:np.arange(ftab.shape[ix])+1.})
         if zdim != None and iz != None and az != 'olevel' :
             ftab = ftab.rename( {az:zdim}) 
-    return ftab
-
-def U2F (utab, nperio=None, psgn=-1.0, zdim='depthf') :
+    return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)
+
+#@numba.jit(forceobj=True)
+def U2F (utab, nperio=None, psgn=-1.0, zdim='depthf', action='ave') :
     '''Interpolate an array from U grid to F grid i-mean)'''
-    math = __math__ (utab)
-    utab_0 = math.where ( np.isnan(utab), 0., utab)
-    utab_0 = lbc (utab_0 , nperio=nperio, cd_type='U', psgn=psgn)
+    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)
     iy, ay = __findAxis__ (utab_0, 'y')
     iz, az = __findAxis__ (utab_0, 'z')
-    ftab = lbc ( 0.5 * (utab_0 + np.roll (utab_0, axis=iy, shift=-1)), cd_type='F', nperio=nperio, psgn=psgn)
-
-    if math == xr :
+    if action == 'ave'  :    ftab = 0.5 *      (utab_0 + np.roll (utab_0, axis=iy, shift=-1))
+    if action == 'min'  :    ftab = np.minimum (utab_0 , np.roll (utab_0, axis=iy, shift=-1))
+    if action == 'max'  :    ftab = np.maximum (utab_0 , np.roll (utab_0, axis=iy, shift=-1))
+    if action == 'mult' :    ftab =             utab_0 * np.roll (utab_0, axis=iy, shift=-1)
+    ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)
+
+    if mmath == xr :
         if ay != None : 
             ftab = ftab.assign_coords({'y':np.arange(ftab.shape[iy])+1.})
@@ -1261,30 +1415,39 @@
     return ftab
 
-def F2T (ftab, nperio=None, psgn=1.0, zdim='deptht') :
+#@numba.jit(forceobj=True)
+def F2T (ftab, nperio=None, psgn=1.0, zdim='deptht', action='ave') :
     '''Interpolate an array on F grid to T grid (i- and j- means)'''
-    math = __math__ (ftab)
-    ftab_0 = math.where ( np.isnan(ttab), 0., ttab)
-    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
-    ttab = U2T(F2U(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim), nperio=nperio, psgn=psgn, zdim=zdim)
-    return ttab
-
-def T2F (ttab, nperio=None, psgn=1.0, zdim='deptht') :
+    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)
+
+#@numba.jit(forceobj=True)
+def T2F (ttab, nperio=None, psgn=1.0, zdim='deptht', action='mean') :
     '''Interpolate an array on T grid to F grid (i- and j- means)'''
-    math = __math__ (ftab)
-    ttab_0 = math.where ( np.isnan(ttab), 0., ttab)
-    ttab_0 = lbc (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn)
-    ftab = T2U(U2F(ttab, nperio=nperio, psgn=psgn, zdim=zdim), nperio=nperio, psgn=psgn, zdim=zdim)
-    return ftab
-
-def F2U (ftab, nperio=None, psgn=1.0, zdim='depthu') :
+    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)
+
+#@numba.jit(forceobj=True)
+def F2U (ftab, nperio=None, psgn=1.0, zdim='depthu', action='ave') :
     '''Interpolate an array on F grid to FUgrid (i-mean)'''
-    math = __math__ (ftab)
-    ftab_0 = math.where ( np.isnan(ftab), 0., ftab)
-    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
+    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)
     iy, ay = __findAxis__ (ftab_0, 'y')
     iz, az = __findAxis__ (ftab_0, 'z')
-    utab = lbc ( 0.5 * (ftab_0 + np.roll (ftab_0, axis=iy, shift=-1)), cd_type='U', nperio=nperio, psgn=psgn)
-    
-    if math == xr :
+    if action == 'ave'  : utab = 0.5 *      (ftab_0 + np.roll (ftab_0, axis=iy, shift=-1))
+    if action == 'min'  : utab = np.minimum (ftab_0 , np.roll (ftab_0, axis=iy, shift=-1))
+    if action == 'max'  : utab = np.maximum (ftab_0 , np.roll (ftab_0, axis=iy, shift=-1))
+    if action == 'mult' : utab =             ftab_0 * np.roll (ftab_0, axis=iy, shift=-1)
+
+    utab = lbc_del (utab, nperio=nperio, cd_type='U', psgn=psgn)
+    
+    if mmath == xr :
         utab = utab.assign_coords({ay:np.arange(ftab.shape[iy])+1.})
         if zdim != None and iz != None and az != 'olevel' :
@@ -1292,14 +1455,19 @@
     return utab
 
-def F2V (ftab, nperio=None, psgn=1.0, zdim='depthv') :
+#@numba.jit(forceobj=True)
+def F2V (ftab, nperio=None, psgn=1.0, zdim='depthv', action='ave') :
     '''Interpolate an array from F grid to V grid (i-mean)'''
-    math = __math__ (ftab)
-    ftab_0 = math.where ( np.isnan(ftab), 0., ftab)
-    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn)
+    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')
     iz, az = __findAxis__ (ftab_0, 'z')
-    vtab = lbc ( 0.5 * (ftab_0 + np.roll (ftab_0, axis=ix, shift=-1)), cd_type='V', nperio=nperio)
-
-    if math == xr :
+    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)
+    if mmath == xr :
         vtab = vtab.assign_coords({ax:np.arange(ftab.shape[ix])+1.})
         if zdim != None and iz != None and az != 'olevel' :
@@ -1307,4 +1475,5 @@
     return vtab
 
+#@numba.jit(forceobj=True)
 def W2T (wtab, zcoord=None, zdim='deptht', sval=np.nan) :
     '''
@@ -1312,6 +1481,6 @@
     sval is the bottom value
     '''
-    math = __math__ (wtab)
-    wtab_0 = math.where ( np.isnan(wtab), 0., wtab)
+    mmath = __mmath__ (wtab)
+    wtab_0 = mmath.where ( np.isnan(wtab), 0., wtab)
 
     iz, az = __findAxis__ (wtab_0, 'z')
@@ -1319,5 +1488,5 @@
     ttab = 0.5 * ( wtab_0 + np.roll (wtab_0, axis=iz, shift=-1) )
     
-    if math == xr :
+    if mmath == xr :
         ttab[{az:iz}] = sval
         if zdim != None and iz != None and az != 'olevel' :
@@ -1330,4 +1499,5 @@
     return ttab
 
+#@numba.jit(forceobj=True)
 def T2W (ttab, zcoord=None, zdim='depthw', sval=np.nan, extrap_surf=False) :
     '''Interpolate an array from T grid to W grid (k-mean)
@@ -1335,10 +1505,10 @@
     if extrap_surf==True, surface value is taken from 1st level value.
     '''
-    math = __math__ (ttab)
-    ttab_0 = math.where ( np.isnan(ttab), 0., ttab)
+    mmath = __mmath__ (ttab)
+    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
     iz, az = __findAxis__ (ttab_0, 'z')
     wtab = 0.5 * ( ttab_0 + np.roll (ttab_0, axis=iz, shift=1) )
 
-    if math == xr :
+    if mmath == xr :
         if extrap_surf : wtab[{az:0}] = ttabb[{az:0}]
         else           : wtab[{az:0}] = sval
@@ -1347,5 +1517,5 @@
         else           : wtab[..., 0, :, :] = sval
 
-    if math == xr :
+    if mmath == xr :
         if zdim != None and iz != None and az != 'olevel' :
                 wtab = wtab.rename ( {az:zdim})
@@ -1354,7 +1524,8 @@
     return wtab
 
-def fill (ptab, nperio, cd_type='T', npass=1) :
-    '''
-    Fill 0. or np.nan values with mean of neighbours
+#@numba.jit(forceobj=True)
+def fill (ptab, nperio, cd_type='T', npass=1, sval=0.) :
+    '''
+    Fill sval values with mean of neighbours
    
     Inputs :
@@ -1363,16 +1534,27 @@
     '''       
 
-    math = __math__ (ptab)
-    ztab   = math.where (np.isnan(ptab), 0., ptab)
-
-    DoPerio = False
-    if nperio == 4.2 or nperio == 6.2 : DoPerio = True
-        
-    if DoPerio :  ztab = lbc_add (ztab)
+    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  = math.where (ztab==0., 0., 1.  )
+        zmask = mmath.where ( np.isnan(ztab), 0., 1.   )
+        ztab0 = mmath.where ( np.isnan(ztab), 0., ztab )
         # Compte du nombre de voisins
-        zcount = zmask \
+        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) \
@@ -1381,26 +1563,26 @@
                 + 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 = zmask \
-          + np.roll(ztab, shift=1, axis=-1) + np.roll(ztab, shift=-1, axis=-1) \
-          + np.roll(ztab, shift=1, axis=-2) + np.roll(ztab, shift=-1, axis=-2) \
-          + 0.5 * ( \
-            + np.roll(np.roll(ztab, shift= 1, axis=-2), shift= 1, axis=-1) \
-            + np.roll(np.roll(ztab, shift=-1, axis=-2), shift= 1, axis=-1) \
-            + np.roll(np.roll(ztab, shift= 1, axis=-2), shift=-1, axis=-1) \
-            + np.roll(np.roll(ztab, shift=-1, axis=-2), shift=-1, axis=-1) )
-
-        zcount = lbc (zcount, nperio=nperio, cd_type=cd_type)
-        znew   = lbc (znew  , nperio=nperio, cd_type=cd_type)
-    
-        ztab = math.where (np.logical_and (zmask==0., zcount>0), znew/zcount, ztab)
-        
-    ztab = math.where (ztab==0., np.nan, ztab)
-
-    if DoPerio : ztab = lbc_del (ztab)
+                + 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
 
+#@numba.jit(forceobj=True)
 def correct_uv (u, v, lat) :
     '''
@@ -1417,7 +1599,6 @@
 
     Outputs : 
-       modified eastward/nothward components have correct polar projection in cartopy
-    '''
-    math = __math__ (u)
+       modified eastward/nothward components to have correct polar projections in cartopy
+    '''
     uv = np.sqrt (u*u + v*v)           # Original modulus
     zu = u
@@ -1427,4 +1608,127 @@
     return uc, vc
 
+def msf (v_e1v_e3v, lat1d, depthw) :
+    '''
+    Computes the meridonal stream function
+    First input is meridional_velocity*e1v*e3v
+    '''
+    @numba.jit(forceobj=True)
+    def iin (tab, dim) :
+        '''
+        Integrate from the bottom
+        '''
+        result = tab * 0.0
+        nlen = len(tab.coords[dim])
+        for jn in np.arange (nlen-2, 0, -1) :
+            result [{dim:jn}] = result [{dim:jn+1}] - tab [{dim:jn}]
+        result = result.where (result !=0, np.nan)
+        return result
+    
+    zomsf = iin ((v_e1v_e3v).sum (dim='x', keep_attrs=True)*1E-6, dim='depthv')
+    zomsf = zomsf.assign_coords ( {'depthv':depthw.values, 'y':lat1d})
+    zomsf = zomsf.rename ( {'depthv':'depthw', 'y':'lat'})
+    zomsf.attrs['long_name'] = 'Meridional stream function'
+
+    zomsf.attrs['units'] = 'Sv'
+    zomsf.depthw.attrs=depthw.attrs
+    zomsf.lat.attrs=lat1d.attrs
+        
+    return zomsf
+
+def bsf (u_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':5, 'y':120} )
+    '''
+    @numba.jit(forceobj=True)
+    def iin (tab, dim) :
+        '''
+        Integrate from the south
+        '''
+        result = tab * 0.0
+        nlen = len(tab.coords[dim])
+        for jn in np.arange (3, nlen) :
+            result [{dim:jn}] = result [{dim:jn-1}] + tab [{dim:jn}]
+        return result
+    
+    bsf = iin ((u_e2u_e3u).sum(dim='depthu', keep_attrs=True)*1E-6, dim='y')
+    bsf.attrs = u_e2u_e3u.attrs
+    if bsf0 != None :
+        bsf = bsf - bsf.isel (bsf0)
+       
+    bsf = bsf.where (mask !=0, np.nan)
+    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.
+    
+    '''
+
+    if ref != None :
+        if isinstance (ref, str) : nml_ref = f90nml.read (ref)
+        if isinstance (ref, f90nml.namelist.Namelist) : nml_ref = ref
+        
+    if cfg != None :
+        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 != None :
+        list_nml.append (nml_ref) ; list_comment.append ('ref')
+    if cfg != None :
+        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
+
 ## ===========================================================================
 ##
@@ -1433,5 +1737,5 @@
 ## ===========================================================================
 
-def is_orca_north_fold ( Xtest, cname_long='T' ) :
+def __is_orca_north_fold__ ( Xtest, cname_long='T' ) :
     '''
     Ported (pirated !!?) from Sosie
