Changeset 6899 for TOOLS

Timestamp:

07/15/24 11:11:01 (5 months ago)

Author:

omamce

Message:

MOSIAX, by OM : add ORCA2.4.2 case

Location:

TOOLS/MOSAIX

Files:

• Build_coordinates_mask.py (modified) (2 diffs)
• CalvingWeights.py (modified) (3 diffs)
• CreateWeightsMask.bash (modified) (3 diffs)
• README.md (modified) (5 diffs)
• RunoffWeights.py (modified) (2 diffs)
• nemo.py (modified) (132 diffs)

TOOLS/MOSAIX/Build_coordinates_mask.py

@@ -88 +88 @@
     if nperio != 0 :
         if nperio != 6 :
-            print(f'Warning ! model = {model} and ocePerio = {nperio} instead of 6 !')
+            print (f'Warning ! model = {model} and ocePerio = {nperio} instead of 6 !')
     else :
         nperio = 6
@@ -95 +95 @@
     if nperio != 0 :
         if nperio != 4 :
-            print(f'Attention ! model = {model} and ocePerio = {nperio} instead of 4 !')
+            print (f'Attention ! model = {model} and ocePerio = {nperio} instead of 4 !')
     else :
         nperio = 4
 
-print(f' model = {model}\n ocePerio = {nperio}\n bathy = {n_bathy}\n coord = {n_coord}\n fout  = {n_out}')
-print(f' Switchs : nemo4Ubug = {nemoUbug}, straits = {straits}, coordPerio = {coordperio}, maskbathy = {maskbathynemo}')
+print (f' model = {model}\n ocePerio = {nperio}\n bathy = {n_bathy}\n coord = {n_coord}\n fout  = {n_out}')
+print (f' Switchs : nemo4Ubug = {nemoUbug}, straits = {straits}, coordPerio = {coordperio}, maskbathy = {maskbathynemo}')
 
 ##

TOOLS/MOSAIX/CalvingWeights.py

@@ -78 +78 @@
 # Adding arguments
 parser.add_argument ('--oce'     , help='oce model name', type=str, default='eORCA1.2',
-                      choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2',
+                      choices=['ORCA2.3', 'ORCA2.4', 'ORCA2.4.2', 
+                               'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2',
                                'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] )
 parser.add_argument ('--atm'     , type=str, default='ICO40',
@@ -178 +179 @@
 # Periodicity masking for NEMO
 if nperio_dst == 0 :
+    if dst_Name in ['ORCA2.4.2',]                    : nperio_dst = 4.2
     if dst_Name in ['ORCA2.3', 'ORCA2.4']            : nperio_dst = 4
     if dst_Name in ['eORCA1.2', 'eORCA1.3']          : nperio_dst = 6
@@ -356 +358 @@
         #
         #print ( 'Sum of remap_matrix : {:12.3e}'.format(np.sum(matrix_local)) )
-        print ( 'Point in atmospheric grid : {index_src:4d} -- num_links: {num_links:6d}' )
+        print ( f'Point in atmospheric grid : {index_src:4d} -- num_links: {num_links:6d}' )
         print ('  ')
 

TOOLS/MOSAIX/CreateWeightsMask.bash

@@ -75 +75 @@
 echo ${Titre}"Defines model"${Norm}
 # =================================
-CplModel=ORCA2.3xLMD9695
+#CplModel=ORCA2.3xLMD9695
+CplModel=ORCA2.4.2xLMD9695
 #CplModel=ORCA2.3xICO30
 #CplModel=ORCA2.3xICO40
@@ -88 +89 @@
 #CplModel=eORCA025.1xLMD256256
 
-#Version="v0" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142"
+Version="v0" ; Comment="Test ORCA 2.4.2"
 #Version="v1" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142"
 #Version="NoSearchRadius" ; Comment="For testing new routing"
-Version="v3" ; Comment="Correction of ORCA masks to have a perfect conservation of run-off"
+#Version="v3" ; Comment="Correction of ORCA masks to have a perfect conservation of run-off"
 
 # If available, get model name from job name
@@ -326 +327 @@
 
 case ${OCE} in # Periodicity type of ORCA grid
+    ( ORCA2.4.2            ) OcePerio=6.2 ;; # 4 ORCA 6 PALEORCA
     ( ORCA2.4              ) OcePerio=4.2 ;; # 4 ORCA 6 PALEORCA
     ( ORCA2*               ) OcePerio=4   ;; # 4 ORCA 6 PALEORCA

TOOLS/MOSAIX/README.md

@@ -1 +1 @@
 # MOSAIX
+
+## MOSAIX
+
+MOSAIX is a small set of fortran routines, bash scripts, python
+scripts and xml files allowing to generate coupling weights for the
+IPSL Earth System Model. It replaces the old MOSAIC
+software. It uses XIOS as a library to compute weights. In the present
+state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids.
+
+MOSAIX can generate second order conservative interpolation
+weights.  They are used in place of bilinear interpolation by
+[OASIS-MCT](https://portal.enes.org/oasis)
+
+As MOSAIX uses XIOS, it works in parallel, using MPI distributed memory.
+
 
 ## SVN information
@@ -37 +52 @@
 to the following warning :
 
-### Warning
+## Warning
 
 Warning, to install, configure, run, use any of MOSAIX software or to
@@ -47 +62 @@
 software by incorrectly or partially configured
 
-## MOSAIX
-
-MOSAIX is a small set of fortran routines, bash scripts, python
-scripts and xml files allowing to generate coupling weights for the
-IPSL Earth System Model. It is under development, and presently not
-used by the production runs. It aims to replace the old MOSAIC
-software. It uses XIOS as a library to compute weights. In the present
-state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids.
-
-MOSAIX can generate second order conservative interpolation
-weights.  They are used in place of bilinear interpolation by
-[OASIS-MCT](https://portal.enes.org/oasis)
-
-As MOSAIX uses XIOS, it works in parallel, using MPI distributed memory.
 
 ## Compatibility with MOSAIC
@@ -83 +84 @@
 
 MOSAIX been developped and extensively used on Irene Skylake computer,
-using the Intel compilation tools. MOSAIX is supposed to work on any
-computer where XIOS works.
+using the Intel compilation tools. Irene AMD computer is also commonly used. MOSAIX is supposed to work on any computer where XIOS works.
 
 It has been tested on Mac OS X, using `gcc` as a compiler. It is
-fragile, as regularly XIOS updates break the compatibility.
+fragile, as regularly XIOS or Mac OS updates break the compatibility.
 
 ## TODO
@@ -231 +231 @@
   -  _oorc_. ORCA grid at _T_ point. Masked on land, duplicated (from
    periodicity) point masked, with area equal 1.
+   
+## Understanding ORCA periodicity parameter in MOSAIX (`OcePerio`)
+
+ORCA grids are global, with an east-west periodicity, and a special
+fold at the north pole. See https://www.nemo-ocean.eu for mode details.
+
+Depending on the NEMO code version, and ot the ORCA
+configuration, ones may have different situations :
+
+### North pole folding
+- Old versions of ORCA have a north pole folding around a T point.
+- Newer versions of ORCA have a north pole folding around a F point.
+
+### Halos
+- NEMO version < 4.2 have grid file that contains the halo points at
+  the east-west periodidicity and at the upper line.
+- Since version >= 4.2, grid files does not include halo points.
+
+### Values MOSAIX periodicity parameter `OcePerio`
+
+
+                  OcePerio         | NEMO<4.2     | NEMO>=4.2      |
+|----------------------------------|--------------|----------------|
+| Old ORCA grids (T-point folding) | `OcePerio=4` | `OcePerio=4.2` |
+| New ORCA grids (F-point folding) | `OcePerio=6` | `OcePerio=6.2` |
+
 
 ## Generating input grid files

TOOLS/MOSAIX/RunoffWeights.py

@@ -101 +101 @@
 # Adding arguments
 parser.add_argument ('--oce'          , help='oce model name', type=str, default='eORCA1.2',
-                     choices=['ORCA2.3',  'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2',
-                              'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] )
+                     choices=['ORCA2.3', 'ORCA2.4', 'ORCA2.4.2',
+                              'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2',
+                              'eORCA1.4', 'eORCA1.4.2',
+                              'eORCA025', 'eORCA025.1', 'eORCA025.4'] )
 parser.add_argument ('--atm'          , help='atm model name', type=str, default='LMD9695'    )
 parser.add_argument ('--atmCoastWidth', type=int, default=1,
@@ -231 +233 @@
 #    if oce_jpi ==  362 : oce_perio = 6 # ORCA 1
 #    if oce_jpi == 1442 : oce_perio = 6 # ORCA 025
-oce_preio = nemo.__guess_nperio__ (oce_jpj, oce_jpi, nperio=oce_perio)
+oce_perio = nemo.__guess_nperio__ (oce_jpj, oce_jpi, nperio=oce_perio)
 
 print ( f'Oce NPERIO parameter : {oce_perio}' )

TOOLS/MOSAIX/nemo.py

@@ -1 @@
-# -*- coding: utf-8 -*-
+7# -*- coding: utf-8 -*-
 ## ===========================================================================
 ##
@@ -18 +18 @@
 ##
 ## ===========================================================================
-'''Utilities to plot NEMO ORCA fields,
+'''
+Utilities to plot NEMO ORCA fields,
 
 Handles periodicity and other stuff
@@ -38 +39 @@
 import xarray as xr
 
-# Tries to import some moldules that are available in all Python installations
-#  and used in only a few specific functions
 try :
     from sklearn.impute import SimpleImputer
@@ -68 +67 @@
 REPSI = np.finfo (1.0).eps
 
+RAAMO  = xr.DataArray (12)        ; RAAMO .name='RAAMO'  ; RAAMO.attrs.update  ({'units':"mth"    , 'long_name':"Number of months in one year" })
+RJJHH  = xr.DataArray (24)        ; RJJHH .name='RJJHH'  ; RJJHH.attrs.update  ({'units':"h"      , 'long_name':"Number of hours in one day"} )
+RHHMM  = xr.DataArray (60)        ; RHHMM .name='RHHMM'  ; RHHMM.attrs.update  ({'units':"min"    , 'long_name':"Number of minutes in one hour"} )
+RMMSS  = xr.DataArray (60)        ; RMMSS .name='RMMSS'  ; RMMSS.attrs.update  ({'units':"s"      , 'long_name':"Number of seconds in one minute"} )
+RA     = xr.DataArray (6371229.0) ; RA    .name='RA'     ; RA.attrs.update     ({'units':"m"      , 'long_name':"Earth radius"} )
+GRAV   = xr.DataArray (9.80665)   ; GRAV  .name='GRAV'   ; GRAV.attrs.update   ({'units':"m/s2"   , 'long_name':"Gravity"} )
+RT0    = xr.DataArray (273.15)    ; RT0   .name='RT0'    ; RT0.attrs.update    ({'units':"K"      , 'long_name':"Freezing point of fresh water"} )
+RAU0   = xr.DataArray (1026.0)    ; RAU0  .name='RAU0'   ; RAU0.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of sea water"} )
+SICE   = xr.DataArray (6.0)       ; SICE  .name='SICE'   ; SICE.attrs.update   ({'units':"psu"    , 'long_name':"Salinity of ice (for pisces)"} )
+SOCE   = xr.DataArray (34.7)      ; SOCE  .name='SOCE'   ; SOCE.attrs.update   ({'units':"psu"    , 'long_name':"Salinity of sea (for pisces and isf)"} )
+RLEVAP = xr.DataArray (2.5e+6)    ; RLEVAP.name='RLEVAP' ; RLEVAP.attrs.update ({'units':"J/K"    , 'long_name':"Latent heat of evaporation (water)"} )
+VKARMN = xr.DataArray (0.4)       ; VKARMN.name='VKARMN' ; VKARMN.attrs.update ({'units':"-"      , 'long_name':"Von Karman constant"} )
+STEFAN = xr.DataArray (5.67e-8)   ; STEFAN.name='STEFAN' ; STEFAN.attrs.update ({'units':"W/m2/K4", 'long_name':"Stefan-Boltzmann constant"} )
+RHOS   = xr.DataArray (330.)      ; RHOS  .name='RHOS'   ; RHOS.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of snow"} )
+RHOI   = xr.DataArray (917.)      ; RHOI  .name='RHOI'   ; RHOI.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of sea ice"} )
+RHOW   = xr.DataArray (1000.)     ; RHOW  .name='RHOW'   ; RHOW.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of freshwater in melt ponds"} )
+RCND_I = xr.DataArray (2.034396)  ; RCND_I.name='RCND_I' ; RCND_I.attrs.update ({'units':"W/m/J"  , 'long_name':"Thermal conductivity of fresh ice"} )
+RCPI   = xr.DataArray (2067.0)    ; RCPI  .name='RCPI'   ; RCPI.attrs.update   ({'units':"J/kg/K" , 'long_name':"Specific heat of fresh ice"} )
+RLSUB  = xr.DataArray (2.834e+6)  ; RLSUB .name='RLSUB'  ; RLSUB.attrs.update  ({'units':"J/kg"   , 'long_name':"Pure ice latent heat of sublimation"} )
+RLFUS  = xr.DataArray (0.334e+6)  ; RLFUS .name='RLFUS'  ; RLFUS.attrs.update  ({'units':"J/kg"   , 'long_name':"Latent heat of fusion of fresh ice"} )
+RTMLT  = xr.DataArray (0.054)     ; RTMLT .name='RTMLT'  ; RTMLT.attrs.update  ({'units':"C/PSU"  , 'long_name':"Decrease of seawater meltpoint with salinity"} )
+
+RDAY     = RJJHH * RHHMM * RMMSS                ; RDAY  .name='RDAY'   ; RDAY.attrs.update   ({'units':"s"      , 'long_name':"Day length"} )
+RSIYEA   = 365.25 * RDAY * 2. * RPI / 6.283076  ; RSIYEA.name='RSIYEA' ; RSIYEA.attrs.update ({'units':"s"      , 'long_name':"Sideral year length"} )
+RSIDAY   = RDAY / (1. + RDAY / RSIYEA)          ; RSIDAY.name='RSIDAY' ; RSIDAY.attrs.update ({'units':"s"      , 'long_name':"Sideral day length"} )
+ROMEGA   = 2. * RPI / RSIDAY                    ; ROMEGA.name='ROMEGA' ; ROMEGA.attrs.update ({'units':"s-1"    , 'long_name':"Earth rotation parameter"} )
+
 NPERIO_VALID_RANGE = [0, 1, 4, 4.2, 5, 6, 6.2]
 
-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
-RA       = 6371229.0        # Earth radius                                  [m]
-GRAV     =       9.80665    # Gravity                                       [m/s2]
-RT0      =     273.15       # Freezing point of fresh water                 [Kelvin]
-RAU0     =    1026.0        # Volumic mass of sea water                     [kg/m3]
-SICE     =       6.0        # Salinity of ice (for pisces)                  [psu]
-SOCE     =      34.7        # Salinity of sea (for pisces and isf)          [psu]
-RLEVAP   =       2.5e+6     # Latent heat of evaporation (water)            [J/K]
-VKARMN   =       0.4        # Von Karman constant
-STEFAN   =       5.67e-8    # Stefan-Boltzmann constant                     [W/m2/K4]
-RHOS     =     330.         # Volumic mass of snow                          [kg/m3]
-RHOI     =     917.         # Volumic mass of sea ice                       [kg/m3]
-RHOW     =    1000.         # Volumic mass of freshwater in melt ponds      [kg/m3]
-RCND_I   =       2.034396   # Thermal conductivity of fresh ice             [W/m/K]
-RCPI     =    2067.0        # Specific heat of fresh ice                    [J/kg/K]
-RLSUB    =       2.834e+6   # Pure ice latent heat of sublimation           [J/kg]
-RLFUS    =       0.334e+6   # Latent heat of fusion of fresh ice            [J/kg]
-RTMLT    =       0.054      # Decrease of seawater meltpoint with salinity
-
-RDAY     = RJJHH * RHHMM * RMMSS                # Day length               [s]
-RSIYEA   = 365.25 * RDAY * 2. * RPI / 6.283076  # Sideral year length      [s]
-RSIDAY   = RDAY / (1. + RDAY / RSIYEA)          # Sideral day length       [s]
-OMEGA    = 2. * RPI / RSIDAY                    # Earth rotation parameter [s-1]
-
 ## Default names of dimensions
-UDIMS = {'x':'x', 'y':'y', 'z':'olevel', 't':'time_counter'}
-
-## All possibles name of dimensions in Nemo files
+UDIMS = {'x':'x', 'y':'y', 'z':'olevel', 't':'time_counter' }
+
+## All possible names of lon/lat variables in Nemo files
+LONNAME=['nav_lon', 'nav_lon_T', 'nav_lon_U', 'nav_lon_V', 'nav_lon_F', 'nav_lon_W',
+             'nav_lon_grid_T', 'nav_lon_grid_U', 'nav_lon_grid_V', 'nav_lon_grid_F', 'nav_lon_grid_W']
+LATNAME=['nav_lat', 'nav_lat_T', 'nav_lat_U', 'nav_lat_V', 'nav_lat_F', 'nav_lat_W',
+             'nav_lat_grid_T', 'nav_lat_grid_U', 'nav_lat_grid_V', 'nav_lat_grid_F', 'nav_lat_grid_W']
+
+## All possibles names 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', ]
+              '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', ]
+              '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
+TNAME = [ 't', 'T', 'tt', 'TT', 'time', 'time_counter', 'time_centered', 'time', 'TIME', 'TIME_COUNTER', 'TIME_CENTERED', ]
+
+## All possibles name of units of dimensions in NEMO files
 XUNIT = [ 'degrees_east', ]
 YUNIT = [ 'degrees_north', ]
@@ -119 +124 @@
 
 ## All possibles size of dimensions in Orca files
-XLENGTH = [ 180, 182, 360, 362, 1440 ]
-YLENGTH = [ 148, 149, 331, 332 ]
-ZLENGTH = [ 31, 75]
+XLENGTH   = [ 180, 182, 360, 362, 1440, 1442 ]
+YLENGTH   = [ 148, 149, 331, 332 ]
+ZLENGTH   = [ 31, 75]
+XYZLENGTH = [ [180,148,31], [182,149,31], [360,331,75], [362,332,75] ]
+
+## T, S array to plot TS diagrams
+Ttab = np.linspace (-2, 40, 100)
+Stab = np.linspace ( 0, 40, 100)
+
+Ttab = xr.DataArray ( Ttab, dims=('Temperature',), coords=(Ttab,) )
+Stab = xr.DataArray ( Stab, dims=('Salinity'   ,), coords=(Stab,) )
+
+Ttab.attrs.update ( {'unit':'degrees_celcius', 'long_name':'Temperature'} )
+Stab.attrs.update ( {'unit':'PSU',             'long_name':'Salinity'} )
+
+# nemo internal options
+import warnings
+from typing import TYPE_CHECKING, Literal, TypedDict
+
+Stack = list()
+
+if TYPE_CHECKING :
+    Options = Literal [ "Debug", "Trace", "Depth", "Stack" ]
+
+    class T_Options (TypedDict) :
+        Debug = bool
+        Trace = bool
+        Depth = int
+        Stack = list()
+
+OPTIONS = { 'Debug':False, 'Trace':False, 'Depth':-1, 'Stack':list() }
+
+class set_options :
+    """
+    Set options for nemo
+    """
+    def __init__ (self, **kwargs):
+        self.old = {}
+        for k, v in kwargs.items():
+            if k not in OPTIONS:
+                raise ValueError ( f"argument name {k!r} is not in the set of valid options {set(OPTIONS)!r}" )
+            self.old[k] = OPTIONS[k]
+        self._apply_update(kwargs)
+
+    def _apply_update (self, options_dict) : OPTIONS.update (options_dict)
+    def __enter__ (self) : return
+    def __exit__ (self, type, value, traceback) : self._apply_update (self.old)
+
+def get_options () -> dict :
+    """
+    Get options for nemo
+
+    See Also
+    ----------
+    set_options
+
+    """
+    return OPTIONS
+
+def return_stack () :
+    return Stack
+
+def PushStack (string:str) :
+    OPTIONS['Depth'] += 1
+    if OPTIONS['Trace'] : print ( '  '*OPTIONS['Depth'], '-->nemo:', string)
+    Stack.append (string)
+    return
+
+def PopStack (string:str) :
+    if OPTIONS['Trace'] : print ( '  '*OPTIONS['Depth'], '<--nemo:', string)
+    OPTIONS['Depth'] -= 1
+    Stack.pop ()
+    return
 
 ## ===========================================================================
@@ -129 +204 @@
     Returns type
     '''
+    PushStack ( f'__mmath__ ( ptab, {default=} )' )
+    
     mmath = default
-    if isinstance (ptab, xr.core.dataarray.DataArray) :
-        mmath = xr
-    if isinstance (ptab, np.ndarray) :
-        mmath = np
-    if isinstance (ptab, np.ma.MaskType) :
-        mmath = np.ma
-
+    if isinstance (ptab, xr.core.dataset.Dataset)     : mmath = 'dataset'
+    if isinstance (ptab, xr.core.dataarray.DataArray) : mmath = xr
+    if isinstance (ptab, np.ndarray)                  : mmath = np
+    if isinstance (ptab, np.ma.MaskType)              : mmath = np.ma
+
+    PopStack ( f'__math_ : {mmath}' )
     return mmath
 
-def __guess_nperio__ (jpj, jpi, nperio=None, out='nperio') :
+def __guess_nperio__ (jpj:int, jpi:int, nperio=None, out:str='nperio') :
     '''Tries to guess the value of nperio (periodicity parameter.
 
@@ -148 +224 @@
     nperio : periodicity parameter
     '''
+    PushStack ( f'__guess_nperio__ ( {jpj=} {jpi=} {nperio=}, {out=} )' )
     if nperio is None :
         nperio = __guess_config__ (jpj, jpi, nperio=None, out=out)
+    PopStack (  f'__guess_nperio__ : {nperio}' )
     return nperio
 
-def __guess_config__ (jpj, jpi, nperio=None, config=None, out='nperio') :
+def __guess_config__ (jpj:int, jpi:int, nperio:bool=None, config=None, out='nperio') :
     '''Tries to guess the value of nperio (periodicity parameter).
 
@@ -161 +239 @@
     nperio : periodicity parameter
     '''
-    print ( jpi, jpj)
+    PushStack ( f'__guessConfig__ ( {jpj=}, {jpi=}, {nperio=}, {config=}, {out=}')
+    
+    print (jpi, jpj)
     if nperio is None :
         ## Values for NEMO version < 4.2
@@ -195 +275 @@
                                 'nemo.py is not ready for this value' )
 
-    if out == 'nperio' :
+    if out == 'nperio'        :
+        PushStack ( f'__guessConfig__ : {nperio=}' )
         return nperio
-    if out == 'config' :
+    if out == 'config'        :
         return config
-    if out == 'perio'  :
+        PushStack ( f'__guessConfig__ : {config=}' )
+    if out == 'perio'         :
         return iperio, jperio, nfold, nftype
+        PushStack ( f'__guessConfig__ : {iperio=}, {jperio=}, {nfold=}, {nftype=}' )
     if out in ['full', 'all'] :
-        return {'nperio':nperio, 'iperio':iperio, 'jperio':jperio, 'nfold':nfold, 'nftype':nftype}
+        zdict = {'nperio':nperio, 'iperio':iperio, 'jperio':jperio, 'nfold':nfold, 'nftype':nftype}
+        PushStack ( f'__guessConfig__ : {zdict}' )
+        return zdict
 
 def __guess_point__ (ptab) :
@@ -215 +300 @@
     Credits : who is the original author ?
     '''
+    PushStack ( f'__guess_point__ ( ptab )')
 
     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 ('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 is None :
             raise AttributeError ('in nemo module : cd_type not found, and cannot by guessed')
         print ( f'Grid set as {gp} deduced from dims {ptab.dims}' )
+
+        PopStack ( f'__guess_point__ : {gp}' )
         return gp
     else :
+        PopStack ( f'__guess_point__ : cd_type not found, input is not an xarray data' )
         raise AttributeError  ('in nemo module : cd_type not found, input is not an xarray data')
 
-def get_shape ( ptab ) :
+def get_shape (ptab) :
     '''Get shape of ptab return a string with axes names
 
@@ -258 +333 @@
     etc ...
     '''
-
+    PushStack ( f'get_shape ( ptab )' )
     g_shape = ''
-    if __find_axis__ (ptab, 'x')[0] :
-        g_shape = 'X'
-    if __find_axis__ (ptab, 'y')[0] :
-        g_shape = 'Y' + g_shape
-    if __find_axis__ (ptab, 'z')[0] :
-        g_shape = 'Z' + g_shape
-    if __find_axis__ (ptab, 't')[0] :
-        g_shape = 'T' + g_shape
+    if __find_axis__ (ptab, 'x')[0] : g_shape = 'X'
+    if __find_axis__ (ptab, 'y')[0] : g_shape = 'Y' + g_shape
+    if __find_axis__ (ptab, 'z')[0] : g_shape = 'Z' + g_shape
+    if __find_axis__ (ptab, 't')[0] : g_shape = 'T' + g_shape
+        
+    PopStack ( f'get_shape : {g_shape}' )
     return g_shape
 
-def lbc_diag (nperio) :
+def lbc_diag (nperio:int) :
     '''Useful to switch between field with and without halo'''
+    PushStack ( f'lbc_diag ( {nperio=}' )
     lperio, aperio = nperio, False
-    if nperio == 4.2 :
-        lperio, aperio = 4, True
-    if nperio == 6.2 :
-        lperio, aperio = 6, True
+    if nperio == 4.2 : lperio, aperio = 4, True
+    if nperio == 6.2 : lperio, aperio = 6, True
+
+    PopStack ( f'lbc_diag : {lperio}, {aperio}' )
     return lperio, aperio
 
-def __find_axis__ (ptab, axis='z', back=True) :
+def __find_axis__ (ptab, axis='z', back:bool=True, verbose:bool=False) :
     '''Returns name and name of the requested axis'''
+    PushStack ( '__find_axis__ ( ptab, {axis=}, {back=}, {verbose=}' )
+    
     mmath = __mmath__ (ptab)
     ax, ix = None, None
@@ -286 +362 @@
     if axis in XNAME :
         ax_name, unit_list, length = XNAME, XUNIT, XLENGTH
+        if verbose : print ( f'Working on xaxis found by name : {axis=} : {XNAME=} {ax_name=} {unit_list=} {length=}' )
     if axis in YNAME :
         ax_name, unit_list, length = YNAME, YUNIT, YLENGTH
+        if verbose : print ( f'Working on yaxis found by name : {axis=} : {YNAME=} {ax_name=} {unit_list=} {length=}' )
     if axis in ZNAME :
         ax_name, unit_list, length = ZNAME, ZUNIT, ZLENGTH
+        if verbose : print ( f'Working on zaxis found by name : {axis=} : {ZNAME=} {ax_name=} {unit_list=} {length=}' )
     if axis in TNAME :
         ax_name, unit_list, length = TNAME, TUNIT, None
-
-    if mmath == xr :
+        if verbose : print ( f'Working on taxis found by name : {axis=} : {TNAME=} {ax_name=} {unit_list=} {length=}' )
+
+    if verbose :
+        print ( f'{ax_name=} {unit_list=} {length=}' )
+            
+    if mmath in [xr, 'dataset' ]: 
         # Try by name
-        for dim in ax_name :
-            if dim in ptab.dims :
-                ix, ax = ptab.dims.index (dim), dim
+        for i, dim in enumerate (ptab.dims) :
+            if verbose : print ( f'{i=} {dim=}' )
+            if dim in ax_name :
+                if verbose : print ( f'Rule 2 : {dim=} axis found by name : {XNAME=}' )
+                ix, ax = i, dim
 
         # If not found, try by axis attributes
         if not ix :
+            if verbose : print ( 'ix not found - 1' )
             for i, dim in enumerate (ptab.dims) :
+                if verbose : print ( f'{i=} {dim=}' )
                 if 'axis' in ptab.coords[dim].attrs.keys() :
                     l_axis = ptab.coords[dim].attrs['axis']
-                    if axis in ax_name and l_axis == 'X' :
+                    if verbose : print ( f'Rule 3 : Trying {i=} {dim=} {l_axis=}' )
+                    if l_axis in ax_name and l_axis == 'X' :
+                        if verbose : print ( f'Rule 3 : xaxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' )
                         ix, ax = (i, dim)
-                    if axis in ax_name and l_axis == 'Y' :
+                    if l_axis in ax_name and l_axis == 'Y' :
+                        if verbose : print ( f'Rule 3 : yaxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' )
                         ix, ax = (i, dim)
-                    if axis in ax_name and l_axis == 'Z' :
+                    if l_axis in ax_name and l_axis == 'Z' :
+                        if verbose : print ( f'Rule 3 : zaxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' )
                         ix, ax = (i, dim)
-                    if axis in ax_name and l_axis == 'T' :
+                    if l_axis in ax_name and l_axis == 'T' :
+                        if verbose : print ( f'Rule 3 : taxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' )
                         ix, ax = (i, dim)
 
@@ -319 +411 @@
                     for name in unit_list :
                         if name in ptab.coords[dim].attrs['units'] :
+                            if verbose : print ( f'Rule 4 : {name=} found by unit : {axis=} : {unit_list=} {i=} {dim=}' )
                             ix, ax = i, dim
 
     # If numpy array or dimension not found, try by length
-    if mmath != xr or not ix :
+    
+    if mmath not in [xr, 'dataset'] or not ix :
         if length :
-            l_shape = ptab.shape
+            if mmath in [xr, 'dataset'] :
+                l_shape =[ len (x) for x in ptab.dims ]
+            else :
+                l_shape = ptab.shape
             for nn in np.arange ( len(l_shape) ) :
                 if l_shape[nn] in length :
-                    ix = nn
+                    ix = nn ; ax = None
+                    if verbose : print ( f'Rule 5 : {ax_name=} axis found by length : {axis=} : {XNAME=} {ix=} {ax=}' )
 
     if ix and back :
-        ix -= len(ptab.shape)
-
+        if mmath in [xr, 'dataset'] :
+            ix -= len (ptab.dims)
+        else :
+            ix -= len (ptab.shape)
+
+    PopStack ( f'__find_axis__ : {ax}, {ix}' )
     return ax, ix
 
-def find_axis ( ptab, axis='z', back=True ) :
+def find_axis (ptab, axis:str='z', back:bool=True, verbose:bool=False) :
     '''Version of find_axis with no __'''
-    ix, xx = __find_axis__ (ptab, axis, back)
+    PushStack ( f'find_axis ( ptab, {axis=}, {back=}, {varbose=}' )
+    
+    ix, xx = __find_axis__ (ptab, axis, back, verbose)
+    
+    PopStack ( f'find_axis : {xx}, {ix}' )
     return xx, ix
 
-def fixed_lon (plon, center_lon=0.0) :
+def fixed_lon (plon, center_lon:float=0.0) :
     '''Returns corrected longitudes for nicer plots
 
@@ -348 +454 @@
     See https://gist.github.com/pelson/79cf31ef324774c97ae7
     '''
+    PushStack ( f'fixed_lon ( plon, {center_lon=}' )
     mmath = __mmath__ (plon)
 
@@ -359 +466 @@
 
     # Special case for eORCA025
-    if f_lon.shape [-1] == 1442 :
-        f_lon [..., -2, :] = f_lon [..., -3, :]
-    if f_lon.shape [-1] == 1440 :
-        f_lon [..., -1, :] = f_lon [..., -2, :]
-
-    if f_lon.min () > center_lon :
-        f_lon += -360.0
-    if f_lon.max () < center_lon :
-        f_lon +=  360.0
-
-    if f_lon.min () < center_lon-360.0 :
-        f_lon +=  360.0
-    if f_lon.max () > center_lon+360.0 :
-        f_lon += -360.0
-
+    if f_lon.shape [-1] == 1442 : f_lon [..., -2, :] = f_lon [..., -3, :]
+    if f_lon.shape [-1] == 1440 : f_lon [..., -1, :] = f_lon [..., -2, :]
+
+    if f_lon.min () > center_lon       : f_lon += -360.0
+    if f_lon.max () < center_lon       : f_lon +=  360.0
+
+    if f_lon.min () < center_lon-360.0 : f_lon +=  360.0
+    if f_lon.max () > center_lon+360.0 : f_lon += -360.0
+
+    PopStack ( f'fixed_lon : f_lon' )
     return f_lon
 
-def bounds_clolon ( pbounds_lon, plon, rad=False, deg=True) :
+def bounds_clolon (pbounds_lon, plon, rad:bool=False, deg:bool=True) :
     '''Choose closest to lon0 longitude, adding/substacting 360° if needed
     '''
-
-    if rad :
-        lon_range = 2.0*np.pi
-    if deg :
-        lon_range = 360.0
+    PushStack ( f'bounds_clolon ( pbounds_lon, plon, {rad=}, {deb=}' )
+    if rad : lon_range = 2.0*np.pi
+    if deg : lon_range = 360.0
     b_clolon = pbounds_lon.copy ()
 
@@ -392 +492 @@
                           b_clolon-lon_range,
                           b_clolon )
+
+    PopStack ( f'bounds_clolon : b_clolon' )
     return b_clolon
 
-def unify_dims ( dd, x='x', y='y', z='olevel', t='time_counter', verbose=False ) :
+def unify_dims ( dd: "DataArray or Dataset", x=UDIMS['x'], y=UDIMS['y'], z=UDIMS['z'], t=UDIMS['t'], verbose:bool=False ) -> "DataArray or Dataset" :
     '''Rename dimensions to unify them between NEMO versions
     '''
+    PushStack ( f'unify_dims ( dd, {x=}, {y=}, {z=}, {t=}, {verbose=}' )
     for xx in XNAME :
         if xx in dd.dims and xx != x :
-            if verbose :
-                print ( f"{xx} renamed to {x}" )
+            if verbose : print ( f"{xx} renamed to {x}" )
             dd = dd.rename ( {xx:x})
 
     for yy in YNAME :
         if yy in dd.dims and yy != y  :
-            if verbose :
-                print ( f"{yy} renamed to {y}" )
+            if verbose : print ( f"{yy} renamed to {y}" )
             dd = dd.rename ( {yy:y} )
 
     for zz in ZNAME :
         if zz in dd.dims and zz != z :
-            if verbose :
-                print ( f"{zz} renamed to {z}" )
+            if verbose : print ( f"{zz} renamed to {z}" )
             dd = dd.rename ( {zz:z} )
 
     for tt in TNAME  :
         if tt in dd.dims and tt != t :
-            if verbose :
-                print ( f"{tt} renamed to {t}" )
+            if verbose : print ( f"{tt} renamed to {t}" )
             dd = dd.rename ( {tt:t} )
 
+    PopStack ( f'unify_dims : dd' )
     return dd
 
 
 if SimpleImputer : 
-  def fill_empty (ptab, sval=np.nan, transpose=False) :
+    def fill_empty (ptab, sval:float=np.nan, transpose:bool=False) :
         '''Fill empty values
 
@@ -432 +532 @@
         values
         '''
+        PushStack ( f'fill_empty ( ptab, {sval=} {transpose=}' )
         mmath = __mmath__ (ptab)
 
@@ -446 +547 @@
             ztab.attrs.update (ptab.attrs)
 
+        PopStack ( f'fill_empty : ztab' )
         return ztab
 
+    def fill_latlon (plat, plon, sval=-1) :
+        '''Fill latitude/longitude 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
+        '''
+        PushStack ( f'fill_latlon ( plat, plon {sval=}' )
+
+        from sklearn.impute import SimpleImputer
+        mmath = __mmath__ (plon)
+        
+        imp = SimpleImputer (missing_values=sval, strategy='mean')
+        imp.fit (plon)
+        zlon = imp.transform (plon)
+        imp.fit (plat.T)
+        zlat = imp.transform (plat.T).T
+        
+        if mmath == xr :
+            zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords)
+            zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords)
+            zlon.attrs.update (plon.attrs)
+            zlat.attrs.update (plat.attrs)
+            
+        zlon = fixed_lon (zlon)
+
+        PopStack ( f'fixed_lon')
+        return zlat, zlon
+    
+    def fill_lonlat (plon, plat, 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
+        '''
+        PushStack ( f'fill_lonlat ( plon, plat {sval=}' )
+        from sklearn.impute import SimpleImputer
+        mmath = __mmath__ (plon)
+        
+        imp = SimpleImputer (missing_values=sval, strategy='mean')
+        imp.fit (plon)
+        zlon = imp.transform (plon)
+        imp.fit (plat.T)
+        zlat = imp.transform (plat.T).T
+        
+        if mmath == xr :
+            zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords)
+            zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords)
+            zlon.attrs.update (plon.attrs)
+            zlat.attrs.update (plat.attrs)
+            
+        zlon = fixed_lon (zlon)
+
+        PopStack ( f'fill_lonlat' )
+        return zlon, zlat
+
+    def fill_bounds_lonlat (pbounds_lon, pbounds_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
+        '''
+        PushStack ( f'fill_bounds_lonlat (pbounds_lon, pbounds_lat, {sval=}' )
+        mmath = __mmath__ (pbounds_lon)
+        
+        z_bounds_lon = np.empty_like (pbounds_lon)
+        z_bounds_lat = np.empty_like (pbounds_lat)
+        
+        imp = SimpleImputer (missing_values=sval, strategy='mean')
+        
+        for n in np.arange (4) :
+            imp.fit (pbounds_lon[:,:,n])
+            z_bounds_lon[:,:,n] = imp.transform (pbounds_lon[:,:,n])
+            imp.fit (pbounds_lat[:,:,n].T)
+            z_bounds_lat[:,:,n] = imp.transform (pbounds_lat[:,:,n].T).T
+            
+        if mmath == xr :
+            z_bounds_lon = xr.DataArray (pbounds_lon, dims=pbounds_lon.dims,
+                                            coords=pbounds_lon.coords)
+            z_bounds_lat = xr.DataArray (pbounds_lat, dims=pbounds_lat.dims,
+                                            coords=pbounds_lat.coords)
+            z_bounds_lon.attrs.update (pbounds_lat.attrs)
+            z_bounds_lat.attrs.update (pbounds_lat.attrs)
+
+        PopStack ( 'fill_bounds_lonlat' )
+        return z_bounds_lon, z_bounds_lat
     
 else : 
@@ -461 +651 @@
           values
         '''
+        PushStack ( f'fill_empty [void version] ( ptab, {sval=} {transpose=}' )
+
         print ( 'Error : module sklearn.impute.SimpleImputer not found' )
         print ( 'Can not call fill_empty' )
         print ( 'Call arguments where : ' )
         print ( f'{ptab.shape=} {sval=} {transpose=}' )
-  
-def fill_lonlat (plon, plat, 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__ (plon)
-
-    imp = SimpleImputer (missing_values=sval, strategy='mean')
-    imp.fit (plon)
-    zlon = imp.transform (plon)
-    imp.fit (plat.T)
-    zlat = imp.transform (plat.T).T
-
-    if mmath == xr :
-        zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords)
-        zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords)
-        zlon.attrs.update (plon.attrs)
-        zlat.attrs.update (plat.attrs)
-
-    zlon = fixed_lon (zlon)
-
-    return zlon, zlat
-
-def fill_bounds_lonlat (pbounds_lon, pbounds_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__ (pbounds_lon)
-
-    z_bounds_lon = np.empty ( pbounds_lon.shape )
-    z_bounds_lat = np.empty ( pbounds_lat.shape )
-
-    imp = SimpleImputer (missing_values=sval, strategy='mean')
-
-    for n in np.arange (4) :
-        imp.fit (pbounds_lon[:,:,n])
-        z_bounds_lon[:,:,n] = imp.transform (pbounds_lon[:,:,n])
-        imp.fit (pbounds_lat[:,:,n].T)
-        z_bounds_lat[:,:,n] = imp.transform (pbounds_lat[:,:,n].T).T
-
-    if mmath == xr :
-        z_bounds_lon = xr.DataArray (pbounds_lon, dims=pbounds_lon.dims,
-                                         coords=pbounds_lon.coords)
-        z_bounds_lat = xr.DataArray (pbounds_lat, dims=pbounds_lat.dims,
-                                         coords=pbounds_lat.coords)
-        z_bounds_lon.attrs.update (pbounds_lat.attrs)
-        z_bounds_lat.attrs.update (pbounds_lat.attrs)
-
-    return z_bounds_lon, z_bounds_lat
+        PopStack ( f'fill_empty [void version]' )
+
+
+    def fill_latlon (plat, plon, sval=-1) :
+        ''' Void version of fill_latlon, because module sklearn.impute.SimpleImputer is not available
+
+        Useful when NEMO has run with no wet points options :
+        some parts of the domain, with no ocean points, have no
+        lon/lat values
+        '''
+        PushStack ( f'fill_latlon [void_version] ( plat, plon {sval=}' )
+
+        print ( 'Error : module sklearn.impute.SimpleImputer not found' )
+        print ( 'Can not call fill_empty' )
+        print ( 'Call arguments where : ' )
+        print ( f'{ptab.shape=} {sval=} {transpose=}' )
+
+        PopStack ( f'fixed_lon [void version]')
+        
+    def fill_lonlat (plon, plat, sval=-1) :
+        ''''Void version of fill_lonlat, because module sklearn.impute.SimpleImputer is not available
+        
+        Useful when NEMO has run with no wet points options :
+        some parts of the domain, with no ocean points, have no
+        lon/lat values
+        '''
+        PushStack ( f'fill_lonlat [void Version] ( plon, plat {sval=}' )
+
+        print ( 'Error : module sklearn.impute.SimpleImputer not found' )
+        print ( 'Can not call fill_empty' )
+        print ( 'Call arguments where : ' )
+        print ( f'{ptab.shape=} {sval=} {transpose=}' )
+
+        PopStack ( f'fill_lonlat [void version]' )
+    
+    def fill_bounds_lonlat (pbounds_lon, pbounds_lat, sval=-1) :
+        ''''Void version of fill_bounds_lonlat, because module sklearn.impute.SimpleImputer is not available
+        
+        Useful when NEMO has run with no wet points options :
+        some parts of the domain, with no ocean points, as no
+        lon/lat values
+        '''
+        PushStack ( f'fill_bounds_lonlat [void version] (pbounds_lon, pbounds_lat, {sval=}' )
+
+        print ( 'Error : module sklearn.impute.SimpleImputer not found' )
+        print ( 'Can not call fill_empty' )
+        print ( 'Call arguments where : ' )
+        print ( f'{ptab.shape=} {sval=} {transpose=}' )
+
+        PopStack ( 'fill_bounds_lonlat [void version]' )
 
 def jeq (plat) :
@@ -527 +713 @@
     lat : latitudes of the grid. At least 2D.
     '''
+    PushStack ( f'jeq ( plat ) ' )
     mmath = __mmath__ (plat)
     jy = __find_axis__ (plat, 'y')[-1]
 
-    if mmath == xr :
-        jj = int ( np.mean ( np.argmin (np.abs (np.float64 (plat)),
-                                                axis=jy) ) )
-    else :
-        jj = np.argmin (np.abs (np.float64 (plat[...,:, 0])))
-
+    if mmath == xr : jj = int ( np.mean ( np.argmin (np.abs (np.float64 (plat)), axis=jy) ) )
+    else           : jj = np.argmin (np.abs (np.float64 (plat[...,:, 0])))
+
+    PopStack ( f'jeq : {jj}' )
     return jj
 
@@ -544 +729 @@
     plat (optionnal) : latitudes of the grid
     '''
+    PushStack ( f'lon1d ( plon, plat) ' ) 
     mmath = __mmath__ (plon)
     jpj, jpi  = plon.shape [-2:]
@@ -566 +752 @@
         lon_1d.attrs['long_name :']   = 'Longitude'
 
+    PopStack ( 'lon_1d' )
     return lon_1d
 
@@ -575 +762 @@
     diff [optional] : tolerance
     '''
+    PushStack ( f'latreg ( plat, {diff=}' )
     #mmath = __mmath__ (plat)
     if diff is None :
@@ -587 +775 @@
     lareg  = np.float64 (plat[...,jreg,:].mean(axis=-1))
 
+    PopStack ( f'latreg : {jreg}, lareg' )
     return jreg, lareg
 
@@ -594 +783 @@
     plat : latitudes of the grid (2D)
     '''
+    PushStack ( f'lat1d ( plat )' )
     mmath = __mmath__ (plat)
     iy = __find_axis__ (plat, 'y')[-1]
@@ -606 +796 @@
 
     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
+        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
+    else :  yrange = 90.-lat_reg # V or F grid
 
     if jpj-1 > jreg :
@@ -629 +816 @@
         lat_1d.attrs ['long_name :']   = 'Latitude'
 
+    PopStack ( f'lat1d' )
     return lat_1d
 
@@ -636 +824 @@
     plat, plon : latitudes and longitudes of the grid (2D)
     '''
-    return lat1d (plat),  lon1d (plon, plat)
+    PushStack ( f'latlon1d ( plat, plon) ' )
+    zla = lat1d (plat)
+    zlo = lon1d (plon, plat)
+
+    PopStack ( f'latlon1d' )
+    return zla, zlo
 
 def ff (plat) :
     '''Returns Coriolis factor
     '''
-    zff   = np.sin (RAD * plat) * OMEGA
+    PushStack ( 'ff ( plat )' )
+    zff   = np.sin (RAD * plat) * ROMEGA
+    PopStack ( 'ff' )
     return zff
 
@@ -647 +842 @@
     '''Return Beta factor (derivative of Coriolis factor)
     '''
-    zbeta = np.cos (RAD * plat) * OMEGA / RA
+    PushStack ( 'beta ( plat )' )
+    zbeta = np.cos (RAD * plat) * ROMEGA / RA
+    PopStack ( 'beta' )
     return zbeta
 
@@ -653 +850 @@
     '''Returns masked values outside a lat/lon box
     '''
+    PushStack ( f'mask_lonlat  (ptab, {x0=}, {x1=}, {y0=}, {y1=}, lon, lat, {sval=}' )
     mmath = __mmath__ (ptab)
     if mmath == xr :
@@ -660 +858 @@
     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)))
+                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, sval)
 
+    PopStack ( 'mask_lonlat' )
     return tab
 
@@ -684 +883 @@
 
     '''
+    PushStack ( f'extend ( ptab, {blon=}, {jplus=}, {jpi=}, {nperio=}' )
     mmath = __mmath__ (ptab)
 
-    if ptab.shape[-1] == 1 :
-        tabex = ptab
+    if ptab.shape[-1] == 1 : tabex = ptab
 
     else :
-        if jpi is None :
-            jpi = ptab.shape[-1]
-
-        if blon :
-            xplus = -360.0
-        else   :
-            xplus =    0.0
+        if not jpi : jpi = ptab.shape[-1]
+
+        if blon : xplus = -360.0
+        else    : xplus =    0.0
 
         if ptab.shape[-1] > jpi :
             tabex = ptab
         else :
-            if nperio in [ 0, 4.2 ] :
-                istart, le, la = 0, jpi+1, 0
-            if nperio == 1 :
-                istart, le, la = 0, jpi+1, 0
+            if nperio in [ 0, 4.2 ] : istart, le, la = 0, jpi+1, 0
+            if nperio == 1          : istart, le, la = 0, jpi+1, 0
             if nperio in [4, 6] : # OPA case with two halo points for periodicity
                 # Perfect, except at the pole that should be masked by lbc_plot
@@ -716 +910 @@
                 new_coords = []
                 for coord in ptab.dims :
-                    if coord == lon :
-                        new_coords.append ( np.arange( tabex.shape[-1]))
-                    else            :
-                        new_coords.append ( ptab.coords[coord].values)
+                    if coord == lon : new_coords.append ( np.arange( tabex.shape[-1]))
+                    else            : new_coords.append ( ptab.coords[coord].values)
                 tabex = xr.DataArray ( tabex, dims=ptab.dims,
                                            coords=new_coords )
@@ -727 +919 @@
                      ptab [..., istart+la:istart+la+jplus]          ),
                      axis=-1)
+
+    PopStack ( 'extend' )
     return tabex
 
-def orca2reg (dd, lat_name='nav_lat', lon_name='nav_lon',
-                  y_name='y', x_name='x') :
+def orca2reg (dd, lat_name=None, lon_name=None, y_name=None, x_name=None) :
     '''Assign an ORCA dataset on a regular grid.
 
@@ -741 +934 @@
       Returns : xarray dataset with rectangular grid. Incorrect above 20°N
     '''
+    PushStack ( f'orca2reg ( dd, {lat_name=}, {lon_name=}, {y_name=}, {x_name=}' )
+    if not x_name : x_name, ix = __find_axis__ (dd, axis='x')
+    if not y_name : y_name, jy = __find_axis__ (dd, axis='y')
+
+    if not lon_name :
+        for xn in LONNAME :
+            if xn in dd.variables : lon_name = xn
+    if not lat_name :
+        for yn in LATNAME :
+            if yn in dd.variables : lat_name = yn
+
+    print ( f'{dd.dims=}'   )
+    print ( f'{x_name=} {y_name=}' )
+    
     # Compute 1D longitude and latitude
-    (zlat, zlon) = latlon1d ( dd[lat_name], dd[lon_name])
-
+    ylat, xlon = fill_latlon (dd[lat_name], dd[lon_name], sval=-1)
+    (zlat, zlon) = latlon1d (ylat, xlon)
     zdd = dd
+    
     # Assign lon and lat as dimensions of the dataset
     if y_name in zdd.dims :
@@ -754 +962 @@
     # Force dimensions to be in the right order
     coord_order = ['lat', 'lon']
-    for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z',
+    for dim in [ 'olevel', 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z',
                  'time_counter', 'time', 'tbnds',
-                 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] :
-        if dim in zdd.dims :
-            coord_order.insert (0, dim)
+                 'bnds', 'axis_nbounds', 'nvertex', 'two2', 'two1', 'two', 'four',] :
+        if dim in zdd.dims : coord_order.insert (0, dim)
 
     zdd = zdd.transpose (*coord_order)
+
+    PopStack ( 'orca2reg' )
     return zdd
 
@@ -777 +986 @@
     See NEMO documentation for further details
     '''
+    PushStack ( f'lbc_init ( ptab, {nperio=}' )
     jpi, jpj = None, None
     ax, ix = __find_axis__ (ptab, 'x')
     ay, jy = __find_axis__ (ptab, 'y')
-    if ax :
-        jpi = ptab.shape[ix]
-    if ay :
-        jpj = ptab.shape[jy]
-
-    if nperio is None :
-        nperio = __guess_nperio__ (jpj, jpi, nperio)
+    if ax : jpi = ptab.shape[ix]
+    if ay : jpj = ptab.shape[jy]
+
+    if nperio is None : nperio = __guess_nperio__ (jpj, jpi, nperio)
 
     if nperio not in NPERIO_VALID_RANGE :
         raise AttributeError ( f'{nperio=} is not in the valid range {NPERIO_VALID_RANGE}' )
 
+    PopStack ( f'lbc_init : {jpj}, {jpi}, {nperio}' )
     return jpj, jpi, nperio
 
@@ -803 +1011 @@
     See NEMO documentation for further details
     '''
+    PushStack ( f'lbc_init ( ptab, {nperio=}, {cd_type=}, {psgn=}, {nemo_4u_bug=}' )
+
     jpi, nperio = lbc_init (ptab, nperio)[1:]
-    ax = __find_axis__ (ptab, 'x')[0]
-    ay = __find_axis__ (ptab, 'y')[0]
+    ax     = __find_axis__ (ptab, 'x')[0]
+    ay     = __find_axis__ (ptab, 'y')[0]
     psgn   = ptab.dtype.type (psgn)
     mmath  = __mmath__ (ptab)
 
-    if mmath == xr :
-        ztab = ptab.values.copy ()
-    else           :
-        ztab = ptab.copy ()
+    if mmath == xr : ztab = ptab.values.copy ()
+    else           : ztab = ptab.copy ()
 
     if ax :
@@ -858 +1066 @@
                 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 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
@@ -897 +1099 @@
         ztab.attrs = ptab.attrs
 
+    PopStack ( 'lbc' )
     return ztab
 
@@ -908 +1111 @@
     See NEMO documentation for further details
     '''
+    PushStack ( f'lbc_mask (ptab, {nperio=}, {cd_type=}, {sval=}' )
     jpi, nperio = lbc_init (ptab, nperio)[1:]
     ax = __find_axis__ (ptab, 'x')[0]
@@ -926 +1130 @@
             #> South (in which nperio cases ?)
             # --------------------------------
-            if nperio in [1, 3, 4, 5, 6] :
-                ztab [..., 0, :] = sval
+            if nperio in [1, 3, 4, 5, 6] : ztab [..., 0, :] = sval
 
             #
@@ -952 +1155 @@
                 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 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
@@ -979 +1176 @@
                     ztab [..., -2, jpi//2+1:-1] = sval
 
+    PopStack ( 'lbc_mask' )
     return ztab
 
@@ -994 +1192 @@
     See NEMO documentation for further details
     '''
+    PushStack ( f'lbc_plot (ptab, {nperio=}, {cd_type=}, {psgn=}, {sval=}' )
     jpi, nperio = lbc_init (ptab, nperio)[1:]
     ax = __find_axis__ (ptab, 'x')[0]
@@ -1012 +1211 @@
             #> Masks south
             # ------------
-            if nperio in [4, 6] :
-                ztab [..., 0, : ] = sval
+            if nperio in [4, 6] : ztab [..., 0, : ] = sval
 
             #
@@ -1037 +1235 @@
                 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 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
@@ -1062 +1254 @@
                     ztab [..., -2, jpi//2+1:-1] = sval
 
+    PopStack ( 'lbc_plot' )
     return ztab
 
 def lbc_add (ptab, nperio=None, cd_type=None, psgn=1) :
-    '''Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2
+    '''
+    Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2
 
     Periodicity and north fold halos has been removed in NEMO 4.2
@@ -1076 +1270 @@
     See NEMO documentation for further details
     '''
+    PushStack ( f'lbc_add ( ptab, {nperio=}, {cd_type=}, {psgn=} )' )
     mmath = __mmath__ (ptab)
     nperio = lbc_init (ptab, nperio)[-1]
@@ -1097 +1292 @@
                 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 ()
+                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)
@@ -1106 +1299 @@
                 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 'X' in lshape : ptab_ext [...,      1:-1] = ptab.copy ()
+                if 'Y' in lshape : ptab_ext [..., :-1      ] = ptab.copy ()
 
         if nperio == 4.2 :
@@ -1120 +1311 @@
             az = __find_axis__ (ptab, 'z')[0]
             at = __find_axis__ (ptab, 't')[0]
-            if az :
-                ptab_ext = ptab_ext.assign_coords ( {az:ptab.coords[az]} )
-            if at :
-                ptab_ext = ptab_ext.assign_coords ( {at:ptab.coords[at]} )
+            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)
 
+    PopStack ( 'lbc_add' )
     return ptab_ext
 
@@ -1141 +1331 @@
     See NEMO documentation for further details
     '''
+    PushStack ( f'lbc_del (ptab, {nperio=}, {cd_type=}, {psgn=}' )
     nperio = lbc_init (ptab, nperio)[-1]
     #lshape = get_shape (ptab)
@@ -1148 +1339 @@
     if nperio in [4.2, 6.2] :
         if ax or ay :
-            if ax and ay :
-                ztab = lbc (ptab[..., :-1, 1:-1],
-                            nperio=nperio, cd_type=cd_type, psgn=psgn)
+            if ax and ay : ztab = lbc (ptab[..., :-1, 1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn)
             else :
-                if ax :
-                    ztab = lbc (ptab[...,      1:-1],
-                                nperio=nperio, cd_type=cd_type, psgn=psgn)
-                if ay :
-                    ztab = lbc (ptab[..., -1],
-                                nperio=nperio, cd_type=cd_type, psgn=psgn)
-        else :
-            ztab = ptab
-    else :
-        ztab = ptab
-
+                if ax    : ztab = lbc (ptab[...,      1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn)
+                if ay    : ztab = lbc (ptab[..., -1]       , nperio=nperio, cd_type=cd_type, psgn=psgn)
+        else : ztab = ptab
+    else : ztab = ptab
+
+    PopStack ( 'lbc_del' )
     return ztab
 
@@ -1176 +1360 @@
     See NEMO documentation for further details
     '''
-
+    PushStack ( f'lbc_index ( {jj=}, {ii=}, {jpj=}, {jpi=}, {nperio=}, {cd_type=} )' )
     if nperio is None :
         nperio = __guess_nperio__ (jpj, jpi, nperio)
@@ -1186 +1370 @@
 
     mmath = __mmath__ (jj)
-    if mmath is None :
-        mmath=np
+    if not mmath : mmath=np
 
     #
@@ -1199 +1382 @@
     #
     def mod_ij (cond, jy_new, ix_new) :
+        PushStack ( f'mod_ij (cond, jy_new, ix_new)' )
         jy_r = mmath.where (cond, jy_new, jy)
         ix_r = mmath.where (cond, ix_new, ix)
+        PopStack ( 'mod_ij' )
         return jy_r, ix_r
     #
@@ -1251 +1436 @@
 
     ## 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 ()
-
+    jy += -1 ; ix += -1
+
+    if isinstance (jj, int) : jy = jy.item ()
+    if isinstance (ii, int) : ix = ix.item ()
+
+    PopStack ( f'lbc_index = {jy}, {ix}' )
     return jy, ix
 
-def find_ji (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False, out=None) :
+def find_ji (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False, drop_duplicates=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]
+    Usage: go find_ji  <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask]
     <grid latitudes><grid longitudes> 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)
+    Example : find_ji (40., -20., nav_lat, nav_lon, mask=1.0)
 
     Note : all longitudes and latitudes in degrees
 
-    Note : may work with 1D lon/lat (?)
-    '''
+    Note : may work with 1D lon_data/lat_data (?)
+    '''
+    PushStack ( f'find_ji ( lat_data, lon_data, lat_grid, lon_grid, mask, {verbose=}, {drop_duplicates=}, {out=} ) ')
     # Get grid dimensions
-    if len (lon_grid.shape) == 2 :
-        jpi = lon_grid.shape[-1]
-    else                         :
-        jpi = len(lon_grid)
+    if len (lon_grid.shape) == 2 : jpi = lon_grid.shape[-1]
+    else                         : jpi = len(lon_grid)
 
     #mmath = __mmath__ (lat_grid)
+
+    if verbose :
+        print ( f'{lat_data=}' )
+        print ( f'{lon_data=}' )
 
     # Compute distance from the point to all grid points (in RADian)
@@ -1288 +1474 @@
                + np.cos (RAD*lat_data) * np.cos (RAD*lat_grid) *
                  np.cos(RAD*(lon_data-lon_grid)) )
+
+    if verbose : print ( f'{type(arg)=} {arg.shape=}' )
+    
     # Send masked points to 'infinite'
     distance = np.arccos (arg) + 4.0*RPI*(1.0-mask)
 
-    # Truncates to alleviate some precision problem with some grids
+    if verbose : print ( f'{type(distance)=} {distance.shape=}' )
+
+    # Truncates to alleviate precision problem encountered 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 ())
+    if verbose : print ( f'{type(distance)=} {distance.shape=}' )
+        
+    # Compute index minimum of distance
+    try :
+        nn = len (lat_data)
+    except : 
+        nn = 0
+        jimin = distance.argmin ().astype(int)
+    else : 
+        jimin   = np.empty (nn, dtype=int )
+        for ji in np.arange (nn) :
+            jimin[ji] = distance[ji].argmin()
+    finally : 
+        if verbose :
+            print ( f'{type(jimin)=} {jimin.shape}' )
+            print ( f'{jimin=}' )
 
     # Compute 2D indices (Python/C flavor : starting at 0)
@@ -1304 +1506 @@
     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' \
-            +   f'- Grid:{lat_grid[jmin,imin]:4.1f}°N '   \
-            +   f'{lon_grid[jmin,imin]:4.1f}°E - Dist: {RA*distance[jmin,imin]:6.1f}km' \
-            +   f' {DAR*distance[jmin,imin]:5.2f}° ' \
-            +   f'- Azimuth: {RAD*azimuth:3.2f}RAD - {azimuth:5.1f}°' )
-
+        print ( f'{jmin=}' )
+        print ( f'{imin=}' )
+        
+    if drop_duplicates :
+        zz = np.vstack ( (np.array (jmin), np.array (imin)) )
+        zz = np.swapaxes (zz , 0, 1)
+        zz = np.unique ( zz, axis=0)
+        jmin = zz[:,-2]
+        imin = zz[:,-1]
+
+    PopStack ( 'find_ji' )
     if   out=='dict'                     :
-        return {'x':imin, 'y':jmin}
+        zdict = {UDIMS['y']:jmin, UDIMS['x']:imin}
+        return zdict
     elif out in ['array', 'numpy', 'np'] :
-        return np.array ( [jmin, imin] )
+        return np.array (jmin), np.array (imin)
     elif out in ['xarray', 'xr']         :
-        return xr.DataArray ( [jmin, imin] )
+        jmin = xr.DataArray (jmin, dims=('Num',))
+        imin = xr.DataArray (imin, dims=('Num',))
+        jmin.attrs.update ( {'long_name':'j-index' } )
+        imin.attrs.update ( {'long_name':'i-index' } )
+        jmin.name = 'j_index'
+        imin.name = 'j_index'
+        return jmin, imin
     elif out=='list'                     :
         return [jmin, imin]
@@ -1337 +1541 @@
     '''Returns curl of a vector field
     '''
+    PushStack ( f'curl ( tx, ty, e1f, e2f, {nperio=} )' )
     ax = __find_axis__ (tx, 'x')[0]
     ay = __find_axis__ (ty, 'y')[0]
@@ -1361 +1566 @@
     zcurl = lbc (zcurl, nperio=nperio, cd_type='F', psgn=1)
 
+    PopStack ( 'curl' )
     return zcurl
 
@@ -1366 +1572 @@
     '''Returns divergence of a vector field
     '''
+    PushStack ( f'div  (ux, uy, e1t, e2t, {nperio=}' )
     ax = __find_axis__ (ux, 'x')[0]
     ay = __find_axis__ (ux, 'y')[0]
@@ -1390 +1597 @@
     zdiv = lbc (zdiv, nperio=nperio, cd_type='T', psgn=1)
 
+    PopStack ( 'zdiv' )
     return zdiv
 
@@ -1399 +1607 @@
         glam, gphi : longitude and latitude of the points
     '''
-
+    PushStack ( 'geo2en (pxx, pyy, pzz, glam, gphi)' )
     gsinlon = np.sin (RAD * glam)
     gcoslon = np.cos (RAD * glam)
@@ -1408 +1616 @@
     ptn = - pxx * gcoslon * gsinlat  - pyy * gsinlon * gsinlat + pzz * gcoslat
 
+    PopStack ( geo2en )
     return pte, ptn
 
@@ -1417 +1626 @@
         glam, gphi : longitude and latitude of the points
     '''
-
+    PushStack ( 'en2geo ( pte, ptn, glam, gphi )' )
+    
     gsinlon = np.sin (RAD * glam)
     gcoslon = np.cos (RAD * glam)
@@ -1427 +1637 @@
     pzz =   ptn * gcoslat
 
+    PopStack ( 'en2geo' )
     return pxx, pyy, pzz
 
@@ -1434 +1645 @@
     if needed
     '''
+    PushStack ( f'clo_lon (lon, {lon0=}, {rad=}, {deg=} )' )
+    if rad and deg :
+        raise 
     mmath = __mmath__ (lon, np)
-    if rad :
-        lon_range = 2.*np.pi
-    if deg :
-        lon_range = 360.
+    if rad : lon_range = 2.*np.pi
+    if deg : lon_range = 360.
     c_lon = lon
-    c_lon = mmath.where (c_lon > lon0 + lon_range*0.5,
-                             c_lon-lon_range, clo_lon)
-    c_lon = mmath.where (c_lon < lon0 - lon_range*0.5,
-                             c_lon+lon_range, clo_lon)
-    c_lon = mmath.where (c_lon > lon0 + lon_range*0.5,
-                             c_lon-lon_range, clo_lon)
-    c_lon = mmath.where (c_lon < lon0 - lon_range*0.5,
-                             c_lon+lon_range, clo_lon)
-    if c_lon.shape == () :
-        c_lon = c_lon.item ()
+    c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, c_lon-lon_range, c_lon)
+    c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, c_lon+lon_range, c_lon)
+    c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, c_lon-lon_range, c_lon)
+    c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, c_lon+lon_range, c_lon)
+    if c_lon.shape == () : c_lon = c_lon.item ()
     if mmath == xr :
-        if lon.attrs :
-            c_lon.attrs.update ( lon.attrs )
+        if lon.attrs : c_lon.attrs.update (lon.attrs)
+
+    PopStack ( 'clo_lon' )
     return c_lon
 
@@ -1460 +1668 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpk = gdept_0.shape[0]
     kk = xr.DataArray(pk)
@@ -1474 +1684 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpk  = gdept_0.shape[0]
     if   isinstance (pz, xr.core.dataarray.DataArray ) :
-        zz   = xr.DataArray (pz.values, dims=('zz',))
+        zz   = xr.DataArray (pz.values , dims=('zz',))
     elif isinstance (pz,  np.ndarray) :
         zz   = xr.DataArray (pz.ravel(), dims=('zz',))
@@ -1497 +1709 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpk = gdept_0.shape[0]
     kk = xr.DataArray (pk)
@@ -1512 +1726 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpk = gdept_0.shape[0]
     if isinstance (pz, xr.core.dataarray.DataArray) :
@@ -1540 +1756 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     #print ('start depth2comp')
     if isinstance (pz, xr.core.dataarray.DataArray) :
@@ -1549 +1767 @@
     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
+    if type (pz) in [int, float] : return gz.item()
+    else                         : return gz
 
 def comp2depth (pz, depth0, fact ) :
@@ -1559 +1775 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     if isinstance (pz, xr.core.dataarray.DataArray) :
         zz   = pz.values
@@ -1576 +1794 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpi = plon_1d.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)*plon_1d[i0]+ xx*plon_1d[i1]
+    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)*plon_1d[i0]+ xx*plon_1d[i1]
     return gx.values
 
@@ -1590 +1810 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpi  = plon_1d.shape[0]
     if isinstance (px, xr.core.dataarray.DataArray) :
@@ -1614 +1836 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpj = plat_1d.shape[0]
     jj  = xr.DataArray (pj)
@@ -1628 +1852 @@
     Needed to use transforms in Matplotlib
     '''
+    # No stack here, to avoid problem in Matplotlib
+    
     jpj = plat_1d.shape[0]
     if isinstance (py, xr.core.dataarray.DataArray) :
@@ -1650 +1876 @@
     respect to east
     '''
+    PushStack ( f'angle_full ( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif, {nperio=} )') 
     mmath = __mmath__ (glamt)
 
@@ -1759 +1986 @@
         gcosf = gcosf.assign_coords ( glamf.coords )
 
+    PopStack ( 'angle_full' )
     return gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf
 
@@ -1765 +1993 @@
     respect to east
     '''
+    PushStack ( f'angle (glam, gphi, {nperio=}, {cd_type=} ) ')
     mmath = __mmath__ (glam)
 
@@ -1797 +2026 @@
         gcos = gcos.assign_coords ( glam.coords )
 
+    PopStack ( 'angle' )
     return gsin, gcos
 
-def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type ) :
+def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type='T' ) :
     '''Rotates the Repere: Change vector componantes between
     geographic grid --> stretched coordinates grid.
@@ -1805 +2035 @@
     All components are on the same grid (T, U, V or F)
     '''
-
+    PushStack ( f'rot_en2ij ( u_e, v_n, gsin, gcos, nperio, {cd_type=} )' )
     u_i = + u_e * gcos + v_n * gsin
     v_j = - u_e * gsin + v_n * gcos
@@ -1812 +2042 @@
     v_j = lbc (v_j, nperio=nperio, cd_type=cd_type, psgn=-1.0)
 
+    PopStack ( 'rot_en2ij' )
     return u_i, v_j
 
@@ -1820 +2051 @@
     All components are on the same grid (T, U, V or F)
     '''
+    PushStack ( f'rot_ij2en ( u_i, v_j, gsin, gcos, {nperio=}, {cd_type=} )' )
     u_e = + u_i * gcos - v_j * gsin
     v_n = + u_i * gsin + v_j * gcos
@@ -1826 +2058 @@
     v_n = lbc (v_n, nperio=nperio, cd_type=cd_type, psgn=1.0)
 
+    PopStack ( 'rot_ij2en' )
     return u_e, v_n
 
@@ -1837 +2070 @@
     Returns east-north components on the T grid point
     '''
+    PushStack ( f'rot_uv2en ( uo, vo, gsint, gcost, {nperio=}, {zdim=} )' )
     ut = u2t (uo, nperio=nperio, psgn=-1.0, zdim=zdim)
     vt = v2t (vo, nperio=nperio, psgn=-1.0, zdim=zdim)
@@ -1846 +2080 @@
     v_n = lbc (v_n, nperio=nperio, cd_type='T', psgn=1.0)
 
+    PopStack ( 'rot_uv2en' )
     return u_e, v_n
 
@@ -1857 +2092 @@
     Returns east-north components on the F grid point
     '''
+    PushStack ( f'rot_uv2enf ( uo, vo, gsint, gcost, {nperio=}, {zdim=} )' )
     uf = u2f (uo, nperio=nperio, psgn=-1.0, zdim=zdim)
     vf = v2f (vo, nperio=nperio, psgn=-1.0, zdim=zdim)
@@ -1866 +2102 @@
     v_n = lbc (v_n, nperio=nperio, cd_type='F', psgn= 1.0)
 
+    PopStack ( 'rot_uv2enf' )
     return u_e, v_n
 
@@ -1871 +2108 @@
     '''Interpolates an array from U grid to T grid (i-mean)
     '''
+    PushStack ( f'u2t (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (utab)
     utab_0 = mmath.where ( np.isnan(utab), 0., utab)
@@ -1897 +2135 @@
             if az != zdim :
                 ttab = ttab.rename( {az:zdim})
+                
+    PopStack ( 'u2t' )
     return ttab
 
@@ -1902 +2142 @@
     '''Interpolates an array from V grid to T grid (j-mean)
     '''
+    PushStack ( f'v2t (vtab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (vtab)
     #lperio, aperio = lbc_diag (nperio)
@@ -1927 +2168 @@
             if az != zdim :
                 ttab = ttab.rename( {az:zdim})
+                
+    PopStack ( 'v2t' )
     return ttab
 
@@ -1932 +2175 @@
     '''Interpolates an array from F grid to T grid (i- and j- means)
     '''
+    PushStack ( f'f2t (ftab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (ftab)
     ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab)
@@ -1937 +2181 @@
     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)
+
+    ttab = lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn)
+    
+    PopStack ( 'f2t' )
+    return ttab
 
 def t2u (ttab, nperio=None, psgn=1.0, zdim=None, action='ave') :
     '''Interpolates an array from T grid to U grid (i-mean)
     '''
+    PushStack ( f't2u (ttab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
+
     mmath = __mmath__ (ttab)
     ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
@@ -1966 +2216 @@
             if az != zdim :
                 utab = utab.rename( {az:zdim})
+
+    PopStack ( 't2u' )
     return utab
 
@@ -1971 +2223 @@
     '''Interpolates an array from T grid to V grid (j-mean)
     '''
+    PushStack ( f't2v (ttab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (ttab)
     ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
@@ -1995 +2248 @@
             if az != zdim :
                 vtab = vtab.rename( {az:zdim})
+ 
+    PopStack ( 't2v' )
     return vtab
 
@@ -2000 +2255 @@
     '''Interpolates an array from V grid to F grid (i-mean)
     '''
+    PushStack ( f'v2f (vtab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (vtab)
     vtab_0 = mmath.where ( np.isnan(vtab), 0., vtab)
@@ -2024 +2280 @@
             if az != zdim :
                 ftab = ftab.rename( {az:zdim})
-    return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)
+
+    ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)
+    
+    PopStack ( 'v2f' )
+    return ftab
 
 def u2f (utab, nperio=None, psgn=-1.0, zdim=None, action='ave') :
     '''Interpolates an array from U grid to F grid i-mean)
     '''
+    PushStack ( f'u2f (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (utab)
     utab_0 = mmath.where ( np.isnan(utab), 0., utab)
@@ -2053 +2314 @@
             if az != zdim :
                 ftab = ftab.rename( {az:zdim})
+    
+    PopStack ( 'u2f' )
     return ftab
 
@@ -2058 +2321 @@
     '''Interpolates an array on T grid to F grid (i- and j- means)
     '''
+    PushStack ( f't2f (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (ttab)
     ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
@@ -2064 +2328 @@
                      nperio=nperio, psgn=psgn, zdim=zdim, action=action)
 
-    return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)
+    ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 
+
+    PopStack ( 'v2f' )
+    return ftab
 
 def f2u (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') :
     '''Interpolates an array on F grid to U grid (j-mean)
     '''
+    PushStack ( f'f2u (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (ftab)
     ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab)
@@ -2091 +2359 @@
         if zdim and az and az != zdim :
             utab = utab.rename( {az:zdim})
+
+    PopStack ( 'f2u' )
     return utab
 
@@ -2096 +2366 @@
     '''Interpolates an array from F grid to V grid (i-mean)
     '''
+    PushStack ( f'f2v (ftab, {nperio=}, {psgn=}, {zdim=}, {action=} )' )
     mmath = __mmath__ (ftab)
     ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab)
@@ -2119 +2390 @@
             if az != zdim :
                 vtab = vtab.rename( {az:zdim})
+
+    PopStack ( 'f2v' )
     return vtab
 
@@ -2126 +2399 @@
     sval is the bottom value
     '''
+    PushStack ( f'w2t (wtab, {zcoord=}, {zdim=}, {sval=} )' )
     mmath = __mmath__ (wtab)
     wtab_0 = mmath.where ( np.isnan(wtab), 0., wtab)
@@ -2146 +2420 @@
         ttab[..., -1, :, :] = sval
 
+    PopStack ( 'w2t' )
     return ttab
 
@@ -2154 +2429 @@
     if extrap_surf==True, surface value is taken from 1st level value.
     '''
+    PushStack ( f't2w (utab, {zcoord=}, {zdim=}, {sval=}, {extrap_surf=} )' )
     mmath = __mmath__ (ttab)
     ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab)
@@ -2177 +2453 @@
         else :
             wtab = wtab.assign_coords ( {zdim:np.arange(ttab.shape[kz])+1.} )
+
+    PopStack ( 't2w' )
     return wtab
 
@@ -2186 +2464 @@
        nperio, cd_type : periodicity characteristics
     '''
-
+    PushStack ( f'fill (ptab, {perio=}, {cd_type=}, {npass=}, {sval=} ) ')
     mmath = __mmath__ (ptab)
 
@@ -2237 +2515 @@
         ztab = lbc_del (ztab, nperio=lperio)
 
+    PopStack ( 'fill' )
     return ztab
 
@@ -2256 +2535 @@
        modified eastward/nothward components to have correct polar projections in cartopy
     '''
+    PushStack (f' correct_uv (u, v, lat)')
     uv = np.sqrt (u*u + v*v)           # Original modulus
     zu = u
@@ -2262 +2542 @@
     uc = zu*uv/zz
     vc = zv*uv/zz      # Final corrected values
+
+    PopStack ( 'correct_uv' )
     return uc, vc
 
@@ -2267 +2549 @@
     '''Returns norm of a 2 components vector
     '''
-    return np.sqrt (u*u + v*v)
+    PushStack ( 'norm_uv (u, v)' )
+    zz = np.sqrt (u*u + v*v)
+    Poptack ( 'norm_uv' )
+    return zz
 
 def normalize_uv (u, v) :
     '''Normalizes 2 components vector
     '''
+    PushStack ( 'normalize_uv (u, v)' )
     uv = norm_uv (u, v)
-    return u/uv, v/uv
+    uu = u/uv
+    vv = v/uv
+    PopStack ( 'normalize_uv' )
+    return uu, vv
 
 def msf (vv, e1v_e3v, plat1d, depthw) :
@@ -2281 +2570 @@
     e1v_e3v : prodcut of scale factors e1v*e3v
     '''
-
+    PushStack ( 'msf (vv, e1v_e3v, plat1d, depthw)' )
     v_e1v_e3v = vv * e1v_e3v
     v_e1v_e3v.attrs = vv.attrs
@@ -2301 +2590 @@
     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.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=plat1d.attrs
 
+    PopStack ( 'msf' )
     return zomsf
 
@@ -2318 +2608 @@
          bsf0={'x':5, 'y':300} for eORCA1
     '''
+    PushStack ( f'bsf (uu, e2u_e3u, mask, {nperio=}, {bsf0=} )' )
+    
     u_e2u_e3u       = uu * e2u_e3u
     u_e2u_e3u.attrs = uu.attrs
@@ -2337 +2629 @@
     zbsf = lbc (zbsf, nperio=nperio, cd_type='F')
 
+    PopStack ( 'bsf' )
     return zbsf
 
@@ -2353 +2646 @@
 
         '''
+        PushStack ( f'namelist_read (ref, cfg, {out=}, {flat=}, {verbose=} )' )
         if ref :
             if isinstance (ref, str) :
@@ -2407 +2701 @@
                             print (nam, ':', value, ':', nml[nam][value])
 
+        PopStack ( 'namelist_read' )
         if out == 'dict' :
             return dict_namelist
@@ -2418 +2713 @@
         Read NEMO namelist(s) and return either a dictionnary or an xarray dataset
         '''
+        PushStack ( f'namelist_read [void version] (ref, cfg, {out=}, {flat=}, {verbose=} )' )
+
         print ( 'Error : module f90nml not found' )
         print ( 'Cannot call namelist_read' )
         print ( 'Call parameters where : ')
         print ( f'{err=} {ref=} {cfg=} {out=} {flat=} {verbose=}' )
+        PopStack ( 'namelist_read [void version]' )
 
 def fill_closed_seas (imask, nperio=None,  cd_type='T') :
@@ -2428 +2726 @@
     imask : mask, 1 on ocean, 0 on land
     '''
+    PushStack ( f'fill_closed_seas (imask, {nperio=}, {cd_type=} )' )
     from scipy import ndimage
 
@@ -2433 +2732 @@
     imask_filled = lbc ( imask_filled, nperio=nperio, cd_type=cd_type)
 
+    PopStack ( 'fill_closed_seas' )
     return imask_filled
 
@@ -2438 +2738 @@
 # Sea water state function parameters from NEMO code
 
-RDELTAS = 32.
+RDELTAS = 32.0
 R1_S0   = 0.875/35.16504
-R1_T0   = 1./40.
-R1_Z0   = 1.e-4
+R1_T0   = 1.0/40.
+R1_Z0   = 1.0e-4
 
 EOS000 =  8.0189615746e+02
@@ -2496 +2796 @@
 EOS013 =  3.7969820455e-01
 
-def rhop ( ptemp, psal ) :
+def rhop (ptemp, psal) :
     '''Returns potential density referenced to surface
 
     Computation from NEMO code
     '''
+    PushStack ( 'rhop ( ptemp, psal )' )
     zt      = ptemp * R1_T0                                  # Temperature (°C)
     zs      = np.sqrt ( np.abs( psal + RDELTAS ) * R1_S0 )   # Square root of salinity (PSS)
@@ -2513 +2814 @@
          + (((((EOS600*zs+ EOS500)*zs + EOS400)*zs + EOS300)*zs + EOS200)*zs + EOS100)*zs + EOS000 )
     #
+    PopStack ( 'rhop' )
     return prhop
 
@@ -2520 +2822 @@
     Computation from NEMO code
     '''
+    PushStack ( 'rho ( ptemp, psal) ' )
     zh      = pdep  * R1_Z0                                  # Depth (m)
     zt      = ptemp * R1_T0                                  # Temperature (°C)
@@ -2543 +2846 @@
          + ((((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 )
+                                         EOS200)*zs + EOS100)*zs + EOS000 )
     #
     prho  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
     #
+    PopStack ( 'rho' )
     return prho