Index: TOOLS/MOSAIX/CalvingWeights.py
===================================================================
--- TOOLS/MOSAIX/CalvingWeights.py	(revision 6059)
+++ TOOLS/MOSAIX/CalvingWeights.py	(revision 6064)
@@ -122,5 +122,5 @@
     if dst_Name == 'eORCA025.1' : nperio_dst = 6
 print ('nperio_dst: ' + str(nperio_dst) )
-dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' )
+dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T', sval=0 )
 print ('dst_mskutil : ' + str(np.sum(dst_mskutil)))
 
@@ -131,5 +131,5 @@
 # Preparation by closing some straits
 # -----------------------------------
-if dst_Name in [ 'eORCA025', 'eORCA025.1' ] :
+if dst_Name in [ 'eORCA025', 'eORCA025.1', 'eORCA025.4' ] :
     print ('Filling some seas for ' + dst_Name )
     # Set Gibraltar strait to 0 to fill Mediterranean sea
@@ -144,6 +144,6 @@
     dst_mskutil[997:1012,907] = 0
     
-if dst_Name == 'eORCA1.2' :
-    print ('Filling some seas for eORCA1.2')
+if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4' ] :
+    print ('Filling some seas for ' + dst_Name)
     # Set Gibraltar strait to 0 to fill Mediterrean sea
     dst_mskutil[240, 283] = 0
@@ -157,6 +157,6 @@
     dst_mskutil[284,222:225] = 0
     
-if dst_Name == 'ORCA2.3' :
-    print ('Filling some seas for ORCA2.3')
+if dst_Name in [ 'ORCA2.3', 'ORCA2.4' ] :
+    print ('Filling some seas for ' + dst_Name)
     # Set Gibraltar strait to 0 to fill Mediterrean sea
     dst_mskutil[101,139] = 0
@@ -175,5 +175,5 @@
 # Fill closed seas with image processing library
 # ----------------------------------------------
-dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), nperio=nperio_dst, cd_type='T' )
+dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), nperio=nperio_dst, cd_type='T', sval=0)
 
 # Surfaces
Index: TOOLS/MOSAIX/CreateOasisGrids.bash
===================================================================
--- TOOLS/MOSAIX/CreateOasisGrids.bash	(revision 6059)
+++ TOOLS/MOSAIX/CreateOasisGrids.bash	(revision 6064)
@@ -51,4 +51,5 @@
 if [[ $(hostname) = irene* ]]    ; then arch=irene ; center=tgcc ; fi
 if [[ $(hostname) = lsce3005* ]] ; then arch=spip  ; center=spip ; fi
+if [[ $(hostname) = lsce5138* ]] ; then arch=spip  ; center=spip ; fi
 
 PROGRAM=$(basename ${0})
@@ -152,16 +153,16 @@
     ncks -C --history --append --variable area_grid_${OCEGRID} ${OCE}_coordinates_mask.nc areas_${CplModel}.nc
     # Inverts mask values and switch to integer
-    ncks --history -C --variable mask_${OCEGRID} ${OCE}_coordinates_mask.nc mask_${OCEGRID}_tmp.nc
-    ncatted --history \
-	    --attribute coordinates,mask_${OCEGRID},d,,         \
-	    --attribute online_operation,mask_${OCEGRID},d,,    \
-	    --attribute cell_measures,mask_${OCEGRID},d,,       \
+    ncks --history -C --variable maskutil_${OCEGRID} ${OCE}_coordinates_mask.nc mask_${OCEGRID}_tmp.nc
+    ncatted --history \
+	    --attribute coordinates,maskutil_${OCEGRID},d,,         \
+	    --attribute online_operation,maskutil_${OCEGRID},d,,    \
+	    --attribute cell_measures,maskutil_${OCEGRID},d,,       \
     mask_${OCEGRID}_tmp.nc
     
-    ncap2 --history --append --script "mask_${OCEGRID}=int(1-mask_${OCEGRID});"    mask_${OCEGRID}_tmp.nc masks_${CplModel}.nc
+    ncap2 --history --append --script "maskutil_${OCEGRID}=int(1-maskutil_${OCEGRID});"    mask_${OCEGRID}_tmp.nc masks_${CplModel}.nc
     rm mask_${OCEGRID}_tmp.nc
     ncatted --history \
-	    --attribute long_name,mask_${OCEGRID},o,c,"Land-sea mask" \
-	    --attribute units,mask_${OCEGRID},o,c,"Land:1, Ocean:0"    masks_${CplModel}.nc
+	    --attribute long_name,maskutil_${OCEGRID},o,c,"Land-sea mask" \
+	    --attribute units,maskutil_${OCEGRID},o,c,"Land:1, Ocean:0"    masks_${CplModel}.nc
     # Change order of dimensions
     mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc
@@ -327,23 +328,23 @@
 
 
-echo ${Titre}"Generates grid c${atm}, 'o' meaning 'one'"${Norm}
-# same as t${atm} grid, with surfaces set to grid area
-# and mask to 0 (ocean everywhere, to compute integral over the whole grid))
-# ----------------------------------------------------------------------------
-mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc
-ncap2 --history --script "c${atm}_lon=alon ; c${atm}_lat=alat ; bounds_c${atm}_lon=bounds_lon ; bounds_c${atm}_lat=bounds_lat ; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc
-
-mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc
-ncap2 --history --script "c${atm}_aire=aire"   areas_${CplModel}_tmp.nc areas_${CplModel}.nc
-
-mv masks_${CplModel}.nc masks_${CplModel}_tmp.nc
-ncap2 --history --script "c${atm}_mask=int(Oce2AtmMask)*0+0 ; " masks_${CplModel}_tmp.nc masks_${CplModel}.nc
-
-rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc masks_${CplModel}_tmp.nc
+# echo ${Titre}"Generates grid c${atm}, 'o' meaning 'one'"${Norm}
+# # same as t${atm} grid, with surfaces set to grid area
+# # and mask to 0 (ocean everywhere, to compute integral over the whole grid))
+# # ----------------------------------------------------------------------------
+# mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc
+# ncap2 --history --script "c${atm}_lon=alon ; c${atm}_lat=alat ; bounds_c${atm}_lon=bounds_lon ; bounds_c${atm}_lat=bounds_lat ; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc
+
+# mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc
+# ncap2 --history --script "c${atm}_aire=aire"   areas_${CplModel}_tmp.nc areas_${CplModel}.nc
+
+# mv masks_${CplModel}.nc masks_${CplModel}_tmp.nc
+# ncap2 --history --script "c${atm}_mask=int(Oce2AtmMask)*0+0 ; " masks_${CplModel}_tmp.nc masks_${CplModel}.nc
+
+# rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc masks_${CplModel}_tmp.nc
 
 echo ${Titre}"Generates grid o${oce}, 'o' meaning 'one'"${Norm}
 # same as t${oce} grid, with surfaces set to 1
 # and mask to 0 (ocean everywhere, to compute integral over the whole grid))
-# This grid is used when field are quantities instead of fluxes (i.e river flow)
+# This grid is used when field are quantities instead of fluxes
 # --------------------------------------------------------------------------------------------------------
 mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc
@@ -354,26 +355,26 @@
 
 ncks -C --history --overwrite -v maskutil_T ${OCE}_coordinates_mask.nc maskutil_T.nc
-ncap2 --history --append --script "mask_O=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc
+ncap2 --history --append --script "maskutil_O=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc
 
 rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc 
 
 
-echo ${Titre}"Creates OCEAN C grid : redundant points removed to compute proper integrals"${Norm}
-# --------------------------------------------------------------------------------------------------------
-
-mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc
-ncap2 --history --script "nav_lon_grid_C=nav_lon_grid_T; nav_lat_grid_C=nav_lat_grid_T; bounds_lon_grid_C=bounds_lon_grid_T; bounds_lat_grid_C=bounds_o${oce}_lat=bounds_lat_grid_T; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc
-
-mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc
-ncap2 --history --script "area_grid_C=area_grid_T ; "   areas_${CplModel}_tmp.nc areas_${CplModel}.nc
-
-ncap2 --history --append --script "mask_C=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc
-
-rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc maskutil_T.nc
+# echo ${Titre}"Creates OCEAN C grid : redundant points removed to compute proper integrals"${Norm}
+# # --------------------------------------------------------------------------------------------------------
+
+# mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc
+# ncap2 --history --script "nav_lon_grid_C=nav_lon_grid_T; nav_lat_grid_C=nav_lat_grid_T; bounds_lon_grid_C=bounds_lon_grid_T; bounds_lat_grid_C=bounds_o${oce}_lat=bounds_lat_grid_T; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc
+
+# mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc
+# ncap2 --history --script "area_grid_C=area_grid_T ; "   areas_${CplModel}_tmp.nc areas_${CplModel}.nc
+
+# ncap2 --history --append --script "mask_C=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc
+
+# rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc maskutil_T.nc
 
 #
 echo ${Titre}"Final renaming"${Norm}
 # ----------------------------------------------------------------------------
-for OCEGRID in T U V O C
+for OCEGRID in T U V O
 do
     ocegrid=${OCEGRID~} # To lowercase
@@ -382,5 +383,5 @@
     ncrename --history --variable bounds_lon_grid_${OCEGRID},${ocegrid}${oce}.clo  grids_${CplModel}.nc
     ncrename --history --variable bounds_lat_grid_${OCEGRID},${ocegrid}${oce}.cla  grids_${CplModel}.nc
-    ncrename --history --variable mask_${OCEGRID},${ocegrid}${oce}.msk             masks_${CplModel}.nc
+    ncrename --history --variable maskutil_${OCEGRID},${ocegrid}${oce}.msk         masks_${CplModel}.nc
     ncrename --history --variable area_grid_${OCEGRID},${ocegrid}${oce}.srf        areas_${CplModel}.nc
     
@@ -399,5 +400,5 @@
 ncrename --history --variable bounds_lat,t${atm}.cla          grids_${CplModel}.nc
 
-for ATMGRID in O C ; do
+for ATMGRID in O ; do
     atmgrid=${ATMGRID~} # To lowercase
     ncrename --history --variable ${atmgrid}${atm}_lon,${atmgrid}${atm}.lon         grids_${CplModel}.nc
@@ -409,5 +410,5 @@
 done
     
-for ATMGRID in T O C ; do
+for ATMGRID in T O ; do
     atmgrid=${ATMGRID~} # To lowercase
     ncatted --history \
Index: TOOLS/MOSAIX/CreateWeightsMask.bash
===================================================================
--- TOOLS/MOSAIX/CreateWeightsMask.bash	(revision 6059)
+++ TOOLS/MOSAIX/CreateWeightsMask.bash	(revision 6064)
@@ -68,6 +68,6 @@
 echo ${Titre}"Defines model"${Norm}
 # =================================
-CplModel=ORCA2.3xLMD9695
-#CplModel=ORCA2.3xICO30
+#CplModel=ORCA2.3xLMD9695
+CplModel=ORCA2.3xICO30
 #CplModel=ORCA2.3xICO40
 #CplModel=eORCA1.2xLMD144142
@@ -79,5 +79,6 @@
 #Version="v0" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142"
 #Version="v1" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142"
-Version="NoSearchRadius" ; Comment="For testing new routing"
+#Version="NoSearchRadius" ; Comment="For testing new routing"
+Version="v2" ;  Comment="Correction of ORCA masks to have a perfect conservation of run-off"
 
 # If available, get model name from job name
@@ -104,5 +105,5 @@
 
 # More specific cases
-[[ ${CplModel} = eORCA1.2xLMD144142 ]] && { atmCoastWidth=2 ; oceCoastWidth=2 ; searchRadius=550.0 ; Version="v1" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" ; }
+[[ ${CplModel} = eORCA1.2xLMD144142 ]] && { atmCoastWidth=2 ; oceCoastWidth=2 ; searchRadius=550.0 ; }
 
 if [[ ${CplModel} = eORCA1.2xLMD144142 && Version = "NoSearchRadius" ]] ; then
@@ -178,5 +179,5 @@
 ##
 ## ===========================================================================
-Tag="MOSAIX"
+#Tag="MOSAIX"
 SUBMIT_DIR=$(pwd)
 FMT_XIOS=netcdf4
@@ -589,5 +590,4 @@
 	    --attribute NCO,global,o,c,"NCO netCDF Operator ${NCO} http://nco.sourceforge.net" \
 	    --attribute Python,global,o,c,"Python3 version ${PYTHON_VER}"          \
-	    --attribute OS,global,o,c,"$(uname -o)"                               \
 	    --attribute release,global,o,c,"$(uname -r)"                          \
 	    --attribute directory,global,o,c,"$(pwd)"                             \
@@ -800,5 +800,5 @@
 	      --atmQuantity=${runOff_atmQuantity} --oceQuantity=${runOff_oceQuantity} --ocePerio=${OcePerio}
 fi
-
+exit
 ##
 echo ${Titre}"Calving weights"${Norm}
@@ -909,7 +909,5 @@
 NCO             : NCO netCDF Operator ${NCO} http://nco.sourceforge.net
 Python version  : ${PYTHON_VER}
-OS              : $(uname -o)
 release         : $(uname -r)
-hardware        : $(uname -i)
 EOF
 
Index: TOOLS/MOSAIX/RunoffWeights.py
===================================================================
--- TOOLS/MOSAIX/RunoffWeights.py	(revision 6059)
+++ TOOLS/MOSAIX/RunoffWeights.py	(revision 6064)
@@ -149,10 +149,10 @@
     if oce_jpi == 1442 : oce_perio = 6 # ORCA 025
 print ("Oce NPERIO parameter : {:d}".format(oce_perio))
-oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask, (oce_jpj,oce_jpi)), 'T', oce_perio).ravel()
+oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel()
 oce_address = np.arange(oce_jpj*oce_jpi)
 
 print ("Fill closed sea with image processing library")
 oce_grid_imask2D = np.reshape(oce_grid_pmask,(oce_jpj,oce_jpi))
-oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T' )
+oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T', sval=0 )
 oce_grid_imask = oce_grid_imask2D.ravel()
 ## 
@@ -165,5 +165,5 @@
 oceOceanFiltered2D = ndimage.uniform_filter(oceOcean2D.astype(float), size=NNocean)
 oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D
-oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D,(oce_jpj,oce_jpi)), oce_perio, 'T').ravel()
+oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D,(oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T').ravel()
 
 oceOceanFiltered = oceOceanFiltered2D.ravel()
Index: TOOLS/MOSAIX/nemo.py
===================================================================
--- TOOLS/MOSAIX/nemo.py	(revision 6059)
+++ TOOLS/MOSAIX/nemo.py	(revision 6064)
@@ -18,11 +18,11 @@
 '''
 
-import sys
-
-try    : import numpy as np
-except : pass
-
+import sys, numpy as np
 try    : import xarray as xr
 except : pass
+
+rpi = np.pi ; rad = rpi / 180.0
+
+nperio_valid_range = [0,1,4,5,6]
 
 ## SVN information
@@ -33,5 +33,5 @@
 __HeadURL    = "$HeadURL$"
 
-def __guessNperio__ (jpi, nperio=None) :
+def __guessNperio__ (jpj, jpi, nperio=None) :
     '''
     Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details)
@@ -41,16 +41,50 @@
     nperio : periodicity parameter
     '''
-   
     if nperio == None :
         if jpi ==  182 : nperio = 4 #   ORCA2. We choose legacy orca2.
-        if jpi ==  362 : nperio = 6 #   ORCA1.
-        if jpi == 1442 : nperio = 6 # ORCA025.
+        if jpi ==  362 : nperio = 6 #   eORCA1.
+        if jpi == 1442 : nperio = 6 #   ORCA025.
+
+        if jpj ==  149 : nperio = 4 #   ORCA2. We choose legacy orca2.
+        if jpj ==  332 : nperio = 6 #   eORCA1.
         #
         if nperio == None :
-            sys.exit ('in nemo.lbc : nperio not found, and cannot by guessed' )
+            sys.exit ('in nemo module : nperio not found, and cannot by guessed' )
         else :
-            print ('nperio set as {:d} (deduced from jpi={:d})'.format(nperio, jpi))
+            print ('nperio set as {:d} (deduced from jpi={:d})'.format (nperio, jpi))
 
     return nperio
+
+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
+    if isinstance (ptab, xr.core.dataarray.DataArray) :
+        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 :
+            sys.exit ('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
 
 def fixed_lon (lon) :
@@ -60,8 +94,14 @@
     fixed_lon (lon)
     lon : longitudes of the grid. At least 2D.
+
+    Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7
     '''
     fixed_lon = lon.copy ()
     for i, start in enumerate (np.argmax (np.abs (np.diff (lon, axis=-1)) > 180, axis=-1)) :
-        fixed_lon[...,i, start+1:] += 360
+        fixed_lon[..., i, start+1:] += 360
+
+    # Special case for eORCA025
+    if lon.shape[-1] == 1442 : lon [..., -2, :] =  lon [..., -3, :]
+        
     return fixed_lon
 
@@ -82,5 +122,4 @@
     lat (optionnal) : latitudes of the grid
     '''
-    
     if np.max (lat) != None :
         je     = jeq (lat)
@@ -95,16 +134,19 @@
     return lon1D
 
-def latreg (lat) :
-    '''
-    Returns maximum j index wehre gridlines are along latitude in the norethern hemisphere
-    
-    lon : longitudes of the grid
-    lat (optionnal) : latitudes of the grid
-    '''
-    dy     = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2))))
+def latreg (lat, diff=0.1) :
+    '''
+    Returns maximum j index where gridlines are along latitude in the northern hemisphere
+    
+    lat : latitudes of the grid (2D)
+    diff [optional] : tolerance
+    '''
+    if diff == None :
+        dy = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2))))
+        diff = dy/100.
+    
     je     = jeq (lat)
-    jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< dy/100.)[-1][-1] + je
+    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
+    JREG   = jreg
 
     return jreg, latreg
@@ -114,32 +156,28 @@
     Returns 1D latitudes for zonal means and simple plots.
 
-    lat : latitudes of the grid
+    lat : latitudes of the grid (2D)
     '''
     jpj, jpi = lat.shape[-2:]
     try :
-        if type (lat) ==  xr.core.dataarray.DataArray : math = xr
+        if type (lat) == xr.core.dataarray.DataArray : math = xr
         else : math = np
     except : math = np
 
-    # jreg : index of last uniform latitude, north of equator
-
-    dy = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2))))
+    dy     = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2))))
     je     = jeq (lat)
-    lateq =  np.float64 (lat[...,je,:].mean(axis=-1))
+    lat_eq = np.float64 (lat[...,je,:].mean(axis=-1))
       
-    jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< dy/100.)[-1][-1] + je
-    latreg = np.float64 (lat[...,jreg,:].mean(axis=-1))
-    latave = np.mean (lat, axis=-1)
-
-    if ( np.abs (lateq) <  dy/100. ) : # T, U or W grid
-        dys = (90.-latreg)//(jpj-jreg-1)*0.5
-        yrange = 90.-dys-latreg
-    else                          :  # V or F grid
-        yrange = 90.    -latreg
-        
-    lat1D = math.where (latave<latreg, latave, latreg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) )   
-        
-    if math == xr :
-        lat1D.attrs = lat.attrs
+    jreg, lat_reg = latreg (lat)
+    lat_ave = np.mean (lat, axis=-1)
+
+    if ( np.abs (lat_eq) < dy/100. ) : # T, U or W grid
+        dys    = (90.-lat_reg) / (jpj-jreg-1)*0.5
+        yrange = 90.-dys-lat_reg
+    else                             :  # V or F grid
+        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
 
     return lat1D
@@ -148,38 +186,126 @@
     '''
     Returns simple latitude and longitude (1D) for simple plots.
-    '''
-    laty = lat1D (lat)
-    lonx = lon1D (lon, lat)
-
-    return laty, lonx
-        
-def extend (tab, Lon=False, jplus=25) :
-    '''
-    Returns extended field eastward to have better 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) :
+    try :
+        lon = lon.to_masked_array()
+        lat = lat.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)) )
+    try    : tab = xr.where (mask, ptab, np.nan)
+    except : tab = np.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 (?)
 
     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
     
     '''
-    jpi   = tab.shape[-1]
-    JPLUS = jplus
-    
-    try : 
-        if type (tab) ==  xr.core.dataarray.DataArray : math = xr
-        else : math = np
-    except : math = np
-
-    if Lon : xplus = -360.0
-    else   : xplus =    0.0
-
-    if math == xr : 
-        extend = xr.concat      ( (tab[..., 2:jpi] + xplus, tab[..., 2:jplus]), dim=tab.dims[-1] )
-    if math == np :
-        extend = np.concatenate ( (tab[..., 2:jpi] + xplus, tab[..., 2:jplus]), axis=-1 )
+    if tab.shape[-1] == 1 :
+        extend = tab
+
+    else :
+        if jpi == None : jpi = tab.shape[-1]
+    
+        try : 
+            if type (tab) ==  xr.core.dataarray.DataArray : math = xr
+            else : math = np
+        except : math = np
+
+        if Lon : xplus = -360.0
+        else   : xplus =    0.0
+
+        if tab.shape[-1] > jpi :
+            extend = tab
+        else :
+            if nperio == 1 :
+                istart = 0 ; le=jpi+1 ; la=0
+            if nperio > 1 : # OPA case with to halo points for periodicity
+                istart = 1 ; le=jpi-2 ; la=1  # Perfect, except at the pole that should be masked by lbc_plot
+                
+            if math == xr : 
+                extend = xr.concat      ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), dim=tab.dims[-1] )
+            if math == np :
+                extend = np.concatenate ( (tab[..., istart:istart+le+1] + xplus, tab[..., 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 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, cd_type='T') :
+    """
+    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
+    """
+    jpj, jpi = ptab.shape[-2:]
+    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio)
+    
+    if nperio not in nperio_valid_range :
+        print ( 'nperio=', nperio, ' is not in the valid range', nperio_valid_range )
+        sys.exit()
+
+    if cd_type == None and isinstance (ptab, xr.core.dataarray.DataArray) :
+        cd_type = __guessPoint__ (ptab)
+
+    if cd_type not in ['T', 'U', 'V', 'F', 'W', 'F' ]:
+        print ( 'cd_type=', cd_type, ' is not in the list T, U, V, F, W, F' )
+        sys.exit ()
+
+    return jpj, jpi, nperio, cd_type
+        
         
 def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) :
@@ -189,8 +315,4 @@
       rank 2 at least : ptab[...., lat, lon]
     nperio    : Type of periodicity
-      1, 4, 6 : Cyclic on i dimension (generaly longitudes)
-      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
     psgn      : For change of sign for vector components (1 for scalars, -1 for vector components)
@@ -198,9 +320,7 @@
     See NEMO documentation for further details
     """
-
-    jpi    = ptab.shape[-1]
-    nperio = __guessNperio__ ( jpi, nperio)
-    psgn   = ptab.dtype.type(psgn)
-    
+    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type)
+    psgn   = ptab.dtype.type (psgn)
+
     #
     #> East-West boundary conditions
@@ -210,5 +330,4 @@
         ptab [..., :,  0] = ptab [..., :, -2]
         ptab [..., :, -1] = ptab [..., :,  1]
-   
     #
     #> North-South boundary conditions
@@ -216,13 +335,13 @@
     if nperio in [3, 4] :  # North fold T-point pivot
         if cd_type in [ 'T', 'W' ] : # T-, W-point
-            ptab [..., -1, 1:       ] = psgn * ptab [..., -3, -1:0:-1      ]      
+            ptab [..., -1, 1:       ] = psgn * ptab [..., -3, -1:0:-1]
             ptab [..., -1, 0        ] = psgn * ptab [..., -3, 2            ]
             ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2:0:-1  ]
-                      
+                
         if cd_type == 'U' :
             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ]       
             ptab [..., -1,  0       ] = psgn * ptab [..., -3,  1           ]
             ptab [..., -1, -1       ] = psgn * ptab [..., -3, -2           ]
-
+                
             if nemo_4U_bug :
                 ptab [..., -2, jpi//2+1:-1] = psgn * ptab [..., -2, jpi//2-2:0:-1]
@@ -230,10 +349,10 @@
             else :
                 ptab [..., -2, jpi//2-1:-1] = psgn * ptab [..., -2, jpi//2:0:-1]
-            
+                
         if cd_type == 'V' : 
             ptab [..., -2, 1:       ] = psgn * ptab [..., -3, jpi-1:0:-1   ]
             ptab [..., -1, 1:       ] = psgn * ptab [..., -4, -1:0:-1      ]   
             ptab [..., -1, 0        ] = psgn * ptab [..., -4, 2            ]
-            
+                
         if cd_type == 'F' :
             ptab [..., -2, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ]
@@ -241,5 +360,5 @@
             ptab [..., -1,  0       ] = psgn * ptab [..., -4,  1           ]
             ptab [..., -1, -1       ] = psgn * ptab [..., -4, -2           ] 
-      
+        
     if nperio in [5, 6] :            #  North fold F-point pivot  
         if cd_type in ['T', 'W']  :
@@ -248,18 +367,18 @@
         if cd_type == 'U' :
             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -2, -2::-1       ]       
-            ptab [..., -1, -1       ] = psgn * ptab [..., -2, 0            ]
-             
+            ptab [..., -1, -1       ] = psgn * ptab [..., -2, 0            ] # Bug ?
+                
         if cd_type == 'V' :
             ptab [..., -1, 0:       ] = psgn * ptab [..., -3, -1::-1       ]
             ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2-1::-1 ]
-                             
+                
         if cd_type == 'F' :
             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -2::-1       ]
             ptab [..., -1, -1       ] = psgn * ptab [..., -3, 0            ]
             ptab [..., -2, jpi//2:-1] = psgn * ptab [..., -2, jpi//2-2::-1 ]
-
+                
     return ptab
 
-def lbc_mask (ptab, nperio=None, cd_type='T', sval=None) :
+def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) :
     #
     """
@@ -268,8 +387,4 @@
       rank 2 at least : ptab[...., lat, lon]
     nperio    : Type of periodicity
-      1, 4, 6 : Cyclic on i dimension (generaly longitudes)
-      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
     
@@ -277,8 +392,6 @@
     """
 
-    jpi    = ptab.shape[-1]
-    nperio = __guessNperio__ (jpi, nperio)
-    if sval == None : sval = ptab.dtype.type (0)
-
+    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type)
+         
     #
     #> East-West boundary conditions
@@ -303,22 +416,17 @@
             ptab [..., -1, 0        ] = sval 
             ptab [..., -2, jpi//2:  ] = sval
-                      
+                
         if cd_type == 'U' :
-            ptab [..., -1, 0:-1     ] = sval  
-            ptab [..., -1,  0       ] = sval
-            ptab [..., -1, -1       ] = sval 
-            ptab [..., -2, jpi//2  :] = sval
-            
+            ptab [..., -1,  :       ] = sval  
+            ptab [..., -2, jpi//2:  ] = sval
+                
         if cd_type == 'V' :
-            ptab [..., -2, 1:       ] = sval
-            ptab [..., -1, 1:       ] = sval   
-            ptab [..., -1, 0        ] = sval
-            
+            ptab [..., -2, :       ] = sval
+            ptab [..., -1, :       ] = sval   
+                
         if cd_type == 'F' :
-            ptab [..., -2, 0:-1     ] = sval
-            ptab [..., -1, 0:-1     ] = sval
-            ptab [..., -1,  0       ] = sval
-            ptab [..., -1, -1       ] = sval
-      
+            ptab [..., -2, :     ] = sval
+            ptab [..., -1, :     ] = sval
+    
     if nperio in [5, 6] :            #  North fold F-point pivot
         if cd_type in ['T', 'W']  :
@@ -340,4 +448,106 @@
     return ptab
 
+
+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 (works 
+      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, cd_type = lbc_init (ptab, nperio, cd_type)
+    psgn   = ptab.dtype.type (psgn)
+    
+    #
+    #> East-West boundary conditions
+    # ------------------------------
+    if nperio in [1, 4, 6] :
+        # ... cyclic
+        ptab [..., :,  0] = ptab [..., :, -2]
+        ptab [..., :, -1] = ptab [..., :,  1]
+        
+    #
+    #> North-South boundary conditions
+    # --------------------------------
+    if nperio in [3, 4] :  # North fold T-point pivot
+        if cd_type in [ 'T', 'W' ] : # T-, W-point
+            ptab [..., -1, :] = sval
+            ptab [..., -2, :jpi//2] = sval
+            #ptab [..., -2, :] = sval
+                      
+        if cd_type == 'U' :
+            ptab [..., -1, : ] = sval
+
+        if cd_type == 'V' : 
+            ptab [..., -2, : ] = sval
+            ptab [..., -1, : ] = sval
+            
+        if cd_type == 'F' :
+            ptab [..., -2, : ] = sval
+            ptab [..., -1, : ] = sval
+      
+    if nperio in [5, 6] :            #  North fold F-point pivot  
+        if cd_type in ['T', 'W']  :
+            ptab [..., -1, : ] = sval
+                
+        if cd_type == 'U' :
+            ptab [..., -1, : ] = sval      
+             
+        if cd_type == 'V' :
+            ptab [..., -1, : ] = sval
+            ptab [..., -2, jpi//2:  ] = sval
+                             
+        if cd_type == 'F' :
+            ptab [..., -1, : ] = sval
+            ptab [..., -2, jpi//2:-1] = sval
+
+    return ptab
+
+def geo2oce (pxx, pyy, pzz, glam, gphi) : 
+    '''
+    Change vector from geocentric to east/north
+
+    Input
+        pxx, pyy, pzz : components on the geoce,tric system
+        glam, gphi : longitue 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 oce2geo (pte, ptn, glam, gphi) :
+    ##----------------------------------------------------------------------
+    ##                    ***  ROUTINE oce2geo  ***
+    ##      
+    ## ** Purpose :
+    ##
+    ## ** Method  :   Change vector from east/north to geocentric
+    ##
+    ## History :
+    ##        #         (A. Caubel)  oce2geo - Original code
+    ##----------------------------------------------------------------------
+    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
+
 ## ===========================================================================
 ##
