Changeset 6899


Ignore:
Timestamp:
07/15/24 11:11:01 (10 hours ago)
Author:
omamce
Message:

MOSIAX, by OM : add ORCA2.4.2 case

Location:
TOOLS/MOSAIX
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/MOSAIX/Build_coordinates_mask.py

    r6190 r6899  
    8888    if nperio != 0 : 
    8989        if nperio != 6 : 
    90             print(f'Warning ! model = {model} and ocePerio = {nperio} instead of 6 !') 
     90            print (f'Warning ! model = {model} and ocePerio = {nperio} instead of 6 !') 
    9191    else : 
    9292        nperio = 6 
    … …  
    9595    if nperio != 0 : 
    9696        if nperio != 4 : 
    97             print(f'Attention ! model = {model} and ocePerio = {nperio} instead of 4 !') 
     97            print (f'Attention ! model = {model} and ocePerio = {nperio} instead of 4 !') 
    9898    else : 
    9999        nperio = 4 
    100100 
    101 print(f' model = {model}\n ocePerio = {nperio}\n bathy = {n_bathy}\n coord = {n_coord}\n fout  = {n_out}') 
    102 print(f' Switchs : nemo4Ubug = {nemoUbug}, straits = {straits}, coordPerio = {coordperio}, maskbathy = {maskbathynemo}') 
     101print (f' model = {model}\n ocePerio = {nperio}\n bathy = {n_bathy}\n coord = {n_coord}\n fout  = {n_out}') 
     102print (f' Switchs : nemo4Ubug = {nemoUbug}, straits = {straits}, coordPerio = {coordperio}, maskbathy = {maskbathynemo}') 
    103103 
    104104## 

  • TOOLS/MOSAIX/CalvingWeights.py

    r6666 r6899  
    7878# Adding arguments 
    7979parser.add_argument ('--oce'     , help='oce model name', type=str, default='eORCA1.2', 
    80                       choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 
     80                      choices=['ORCA2.3', 'ORCA2.4', 'ORCA2.4.2',  
     81                               'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 
    8182                               'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
    8283parser.add_argument ('--atm'     , type=str, default='ICO40', 
    … …  
    178179# Periodicity masking for NEMO 
    179180if nperio_dst == 0 : 
     181    if dst_Name in ['ORCA2.4.2',]                    : nperio_dst = 4.2 
    180182    if dst_Name in ['ORCA2.3', 'ORCA2.4']            : nperio_dst = 4 
    181183    if dst_Name in ['eORCA1.2', 'eORCA1.3']          : nperio_dst = 6 
    … …  
    356358        # 
    357359        #print ( 'Sum of remap_matrix : {:12.3e}'.format(np.sum(matrix_local)) ) 
    358         print ( 'Point in atmospheric grid : {index_src:4d} -- num_links: {num_links:6d}' ) 
     360        print ( f'Point in atmospheric grid : {index_src:4d} -- num_links: {num_links:6d}' ) 
    359361        print ('  ') 
    360362 

  • TOOLS/MOSAIX/CreateWeightsMask.bash

    r6666 r6899  
    7575echo ${Titre}"Defines model"${Norm} 
    7676# ================================= 
    77 CplModel=ORCA2.3xLMD9695 
     77#CplModel=ORCA2.3xLMD9695 
     78CplModel=ORCA2.4.2xLMD9695 
    7879#CplModel=ORCA2.3xICO30 
    7980#CplModel=ORCA2.3xICO40 
    … …  
    8889#CplModel=eORCA025.1xLMD256256 
    8990 
    90 #Version="v0" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" 
     91Version="v0" ; Comment="Test ORCA 2.4.2" 
    9192#Version="v1" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" 
    9293#Version="NoSearchRadius" ; Comment="For testing new routing" 
    93 Version="v3" ; Comment="Correction of ORCA masks to have a perfect conservation of run-off" 
     94#Version="v3" ; Comment="Correction of ORCA masks to have a perfect conservation of run-off" 
    9495 
    9596# If available, get model name from job name 
    … …  
    326327 
    327328case ${OCE} in # Periodicity type of ORCA grid 
     329    ( ORCA2.4.2            ) OcePerio=6.2 ;; # 4 ORCA 6 PALEORCA 
    328330    ( ORCA2.4              ) OcePerio=4.2 ;; # 4 ORCA 6 PALEORCA 
    329331    ( ORCA2*               ) OcePerio=4   ;; # 4 ORCA 6 PALEORCA 

  • TOOLS/MOSAIX/README.md

    r6235 r6899  
    11# MOSAIX 
     2 
     3## MOSAIX 
     4 
     5MOSAIX is a small set of fortran routines, bash scripts, python 
     6scripts and xml files allowing to generate coupling weights for the 
     7IPSL Earth System Model. It replaces the old MOSAIC 
     8software. It uses XIOS as a library to compute weights. In the present 
     9state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids. 
     10 
     11MOSAIX can generate second order conservative interpolation 
     12weights.  They are used in place of bilinear interpolation by 
     13[OASIS-MCT](https://portal.enes.org/oasis) 
     14 
     15As MOSAIX uses XIOS, it works in parallel, using MPI distributed memory. 
     16 
    217 
    318## SVN information 
    … …  
    3752to the following warning : 
    3853 
    39 ### Warning 
     54## Warning 
    4055 
    4156Warning, to install, configure, run, use any of MOSAIX software or to 
    … …  
    4762software by incorrectly or partially configured 
    4863 
    49 ## MOSAIX 
    50  
    51 MOSAIX is a small set of fortran routines, bash scripts, python 
    52 scripts and xml files allowing to generate coupling weights for the 
    53 IPSL Earth System Model. It is under development, and presently not 
    54 used by the production runs. It aims to replace the old MOSAIC 
    55 software. It uses XIOS as a library to compute weights. In the present 
    56 state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids. 
    57  
    58 MOSAIX can generate second order conservative interpolation 
    59 weights.  They are used in place of bilinear interpolation by 
    60 [OASIS-MCT](https://portal.enes.org/oasis) 
    61  
    62 As MOSAIX uses XIOS, it works in parallel, using MPI distributed memory. 
    6364 
    6465## Compatibility with MOSAIC 
    … …  
    8384 
    8485MOSAIX been developped and extensively used on Irene Skylake computer, 
    85 using the Intel compilation tools. MOSAIX is supposed to work on any 
    86 computer where XIOS works. 
     86using the Intel compilation tools. Irene AMD computer is also commonly used. MOSAIX is supposed to work on any computer where XIOS works. 
    8787 
    8888It has been tested on Mac OS X, using `gcc` as a compiler. It is 
    89 fragile, as regularly XIOS updates break the compatibility. 
     89fragile, as regularly XIOS or Mac OS updates break the compatibility. 
    9090 
    9191## TODO 
    … …  
    231231  -  _oorc_. ORCA grid at _T_ point. Masked on land, duplicated (from 
    232232   periodicity) point masked, with area equal 1. 
     233    
     234## Understanding ORCA periodicity parameter in MOSAIX (`OcePerio`) 
     235 
     236ORCA grids are global, with an east-west periodicity, and a special 
     237fold at the north pole. See https://www.nemo-ocean.eu for mode details. 
     238 
     239Depending on the NEMO code version, and ot the ORCA 
     240configuration, ones may have different situations : 
     241 
     242### North pole folding 
     243- Old versions of ORCA have a north pole folding around a T point. 
     244- Newer versions of ORCA have a north pole folding around a F point. 
     245 
     246### Halos 
     247- NEMO version < 4.2 have grid file that contains the halo points at 
     248  the east-west periodidicity and at the upper line. 
     249- Since version >= 4.2, grid files does not include halo points. 
     250 
     251### Values MOSAIX periodicity parameter `OcePerio` 
     252 
     253 
     254                  OcePerio         | NEMO<4.2     | NEMO>=4.2      | 
     255|----------------------------------|--------------|----------------| 
     256| Old ORCA grids (T-point folding) | `OcePerio=4` | `OcePerio=4.2` | 
     257| New ORCA grids (F-point folding) | `OcePerio=6` | `OcePerio=6.2` | 
     258 
    233259 
    234260## Generating input grid files 

  • TOOLS/MOSAIX/RunoffWeights.py

    r6666 r6899  
    101101# Adding arguments 
    102102parser.add_argument ('--oce'          , help='oce model name', type=str, default='eORCA1.2', 
    103                      choices=['ORCA2.3',  'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 
    104                               'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
     103                     choices=['ORCA2.3', 'ORCA2.4', 'ORCA2.4.2', 
     104                              'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 
     105                              'eORCA1.4', 'eORCA1.4.2', 
     106                              'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
    105107parser.add_argument ('--atm'          , help='atm model name', type=str, default='LMD9695'    ) 
    106108parser.add_argument ('--atmCoastWidth', type=int, default=1, 
    … …  
    231233#    if oce_jpi ==  362 : oce_perio = 6 # ORCA 1 
    232234#    if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 
    233 oce_preio = nemo.__guess_nperio__ (oce_jpj, oce_jpi, nperio=oce_perio) 
     235oce_perio = nemo.__guess_nperio__ (oce_jpj, oce_jpi, nperio=oce_perio) 
    234236 
    235237print ( f'Oce NPERIO parameter : {oce_perio}' ) 

  • TOOLS/MOSAIX/nemo.py

    r6666 r6899  
    1 # -*- coding: utf-8 -*- 
     17# -*- coding: utf-8 -*- 
    22## =========================================================================== 
    33## 
    … …  
    1818## 
    1919## =========================================================================== 
    20 '''Utilities to plot NEMO ORCA fields, 
     20''' 
     21Utilities to plot NEMO ORCA fields, 
    2122 
    2223Handles periodicity and other stuff 
    … …  
    3839import xarray as xr 
    3940 
    40 # Tries to import some moldules that are available in all Python installations 
    41 #  and used in only a few specific functions 
    4241try : 
    4342    from sklearn.impute import SimpleImputer 
    … …  
    6867REPSI = np.finfo (1.0).eps 
    6968 
     69RAAMO  = xr.DataArray (12)        ; RAAMO .name='RAAMO'  ; RAAMO.attrs.update  ({'units':"mth"    , 'long_name':"Number of months in one year" }) 
     70RJJHH  = xr.DataArray (24)        ; RJJHH .name='RJJHH'  ; RJJHH.attrs.update  ({'units':"h"      , 'long_name':"Number of hours in one day"} ) 
     71RHHMM  = xr.DataArray (60)        ; RHHMM .name='RHHMM'  ; RHHMM.attrs.update  ({'units':"min"    , 'long_name':"Number of minutes in one hour"} ) 
     72RMMSS  = xr.DataArray (60)        ; RMMSS .name='RMMSS'  ; RMMSS.attrs.update  ({'units':"s"      , 'long_name':"Number of seconds in one minute"} ) 
     73RA     = xr.DataArray (6371229.0) ; RA    .name='RA'     ; RA.attrs.update     ({'units':"m"      , 'long_name':"Earth radius"} ) 
     74GRAV   = xr.DataArray (9.80665)   ; GRAV  .name='GRAV'   ; GRAV.attrs.update   ({'units':"m/s2"   , 'long_name':"Gravity"} ) 
     75RT0    = xr.DataArray (273.15)    ; RT0   .name='RT0'    ; RT0.attrs.update    ({'units':"K"      , 'long_name':"Freezing point of fresh water"} ) 
     76RAU0   = xr.DataArray (1026.0)    ; RAU0  .name='RAU0'   ; RAU0.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of sea water"} ) 
     77SICE   = xr.DataArray (6.0)       ; SICE  .name='SICE'   ; SICE.attrs.update   ({'units':"psu"    , 'long_name':"Salinity of ice (for pisces)"} ) 
     78SOCE   = xr.DataArray (34.7)      ; SOCE  .name='SOCE'   ; SOCE.attrs.update   ({'units':"psu"    , 'long_name':"Salinity of sea (for pisces and isf)"} ) 
     79RLEVAP = xr.DataArray (2.5e+6)    ; RLEVAP.name='RLEVAP' ; RLEVAP.attrs.update ({'units':"J/K"    , 'long_name':"Latent heat of evaporation (water)"} ) 
     80VKARMN = xr.DataArray (0.4)       ; VKARMN.name='VKARMN' ; VKARMN.attrs.update ({'units':"-"      , 'long_name':"Von Karman constant"} ) 
     81STEFAN = xr.DataArray (5.67e-8)   ; STEFAN.name='STEFAN' ; STEFAN.attrs.update ({'units':"W/m2/K4", 'long_name':"Stefan-Boltzmann constant"} ) 
     82RHOS   = xr.DataArray (330.)      ; RHOS  .name='RHOS'   ; RHOS.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of snow"} ) 
     83RHOI   = xr.DataArray (917.)      ; RHOI  .name='RHOI'   ; RHOI.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of sea ice"} ) 
     84RHOW   = xr.DataArray (1000.)     ; RHOW  .name='RHOW'   ; RHOW.attrs.update   ({'units':"kg/m3"  , 'long_name':"Volumic mass of freshwater in melt ponds"} ) 
     85RCND_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"} ) 
     86RCPI   = xr.DataArray (2067.0)    ; RCPI  .name='RCPI'   ; RCPI.attrs.update   ({'units':"J/kg/K" , 'long_name':"Specific heat of fresh ice"} ) 
     87RLSUB  = xr.DataArray (2.834e+6)  ; RLSUB .name='RLSUB'  ; RLSUB.attrs.update  ({'units':"J/kg"   , 'long_name':"Pure ice latent heat of sublimation"} ) 
     88RLFUS  = xr.DataArray (0.334e+6)  ; RLFUS .name='RLFUS'  ; RLFUS.attrs.update  ({'units':"J/kg"   , 'long_name':"Latent heat of fusion of fresh ice"} ) 
     89RTMLT  = xr.DataArray (0.054)     ; RTMLT .name='RTMLT'  ; RTMLT.attrs.update  ({'units':"C/PSU"  , 'long_name':"Decrease of seawater meltpoint with salinity"} ) 
     90 
     91RDAY     = RJJHH * RHHMM * RMMSS                ; RDAY  .name='RDAY'   ; RDAY.attrs.update   ({'units':"s"      , 'long_name':"Day length"} ) 
     92RSIYEA   = 365.25 * RDAY * 2. * RPI / 6.283076  ; RSIYEA.name='RSIYEA' ; RSIYEA.attrs.update ({'units':"s"      , 'long_name':"Sideral year length"} ) 
     93RSIDAY   = RDAY / (1. + RDAY / RSIYEA)          ; RSIDAY.name='RSIDAY' ; RSIDAY.attrs.update ({'units':"s"      , 'long_name':"Sideral day length"} ) 
     94ROMEGA   = 2. * RPI / RSIDAY                    ; ROMEGA.name='ROMEGA' ; ROMEGA.attrs.update ({'units':"s-1"    , 'long_name':"Earth rotation parameter"} ) 
     95 
    7096NPERIO_VALID_RANGE = [0, 1, 4, 4.2, 5, 6, 6.2] 
    7197 
    72 RAAMO    =      12          # Number of months in one year 
    73 RJJHH    =      24          # Number of hours in one day 
    74 RHHMM    =      60          # Number of minutes in one hour 
    75 RMMSS    =      60          # Number of seconds in one minute 
    76 RA       = 6371229.0        # Earth radius                                  [m] 
    77 GRAV     =       9.80665    # Gravity                                       [m/s2] 
    78 RT0      =     273.15       # Freezing point of fresh water                 [Kelvin] 
    79 RAU0     =    1026.0        # Volumic mass of sea water                     [kg/m3] 
    80 SICE     =       6.0        # Salinity of ice (for pisces)                  [psu] 
    81 SOCE     =      34.7        # Salinity of sea (for pisces and isf)          [psu] 
    82 RLEVAP   =       2.5e+6     # Latent heat of evaporation (water)            [J/K] 
    83 VKARMN   =       0.4        # Von Karman constant 
    84 STEFAN   =       5.67e-8    # Stefan-Boltzmann constant                     [W/m2/K4] 
    85 RHOS     =     330.         # Volumic mass of snow                          [kg/m3] 
    86 RHOI     =     917.         # Volumic mass of sea ice                       [kg/m3] 
    87 RHOW     =    1000.         # Volumic mass of freshwater in melt ponds      [kg/m3] 
    88 RCND_I   =       2.034396   # Thermal conductivity of fresh ice             [W/m/K] 
    89 RCPI     =    2067.0        # Specific heat of fresh ice                    [J/kg/K] 
    90 RLSUB    =       2.834e+6   # Pure ice latent heat of sublimation           [J/kg] 
    91 RLFUS    =       0.334e+6   # Latent heat of fusion of fresh ice            [J/kg] 
    92 RTMLT    =       0.054      # Decrease of seawater meltpoint with salinity 
    93  
    94 RDAY     = RJJHH * RHHMM * RMMSS                # Day length               [s] 
    95 RSIYEA   = 365.25 * RDAY * 2. * RPI / 6.283076  # Sideral year length      [s] 
    96 RSIDAY   = RDAY / (1. + RDAY / RSIYEA)          # Sideral day length       [s] 
    97 OMEGA    = 2. * RPI / RSIDAY                    # Earth rotation parameter [s-1] 
    98  
    9998## Default names of dimensions 
    100 UDIMS = {'x':'x', 'y':'y', 'z':'olevel', 't':'time_counter'} 
    101  
    102 ## All possibles name of dimensions in Nemo files 
     99UDIMS = {'x':'x', 'y':'y', 'z':'olevel', 't':'time_counter' } 
     100 
     101## All possible names of lon/lat variables in Nemo files 
     102LONNAME=['nav_lon', 'nav_lon_T', 'nav_lon_U', 'nav_lon_V', 'nav_lon_F', 'nav_lon_W', 
     103             'nav_lon_grid_T', 'nav_lon_grid_U', 'nav_lon_grid_V', 'nav_lon_grid_F', 'nav_lon_grid_W'] 
     104LATNAME=['nav_lat', 'nav_lat_T', 'nav_lat_U', 'nav_lat_V', 'nav_lat_F', 'nav_lat_W', 
     105             'nav_lat_grid_T', 'nav_lat_grid_U', 'nav_lat_grid_V', 'nav_lat_grid_F', 'nav_lat_grid_W'] 
     106 
     107## All possibles names of dimensions in Nemo files 
    103108XNAME = [ 'x', 'X', 'X1', 'xx', 'XX', 
    104109              'x_grid_T', 'x_grid_U', 'x_grid_V', 'x_grid_F', 'x_grid_W', 
    105               'lon', 'nav_lon', 'longitude', 'X1', 'x_c', 'x_f', ] 
     110              'lon', 'nav_lon', 'longitude', 'X1', 'x_c', 'x_f'] 
    106111YNAME = [ 'y', 'Y', 'Y1', 'yy', 'YY', 
    107112              'y_grid_T', 'y_grid_U', 'y_grid_V', 'y_grid_F', 'y_grid_W', 
    108               'lat', 'nav_lat', 'latitude' , 'Y1', 'y_c', 'y_f', ] 
     113              'lat', 'nav_lat', 'latitude' , 'Y1', 'y_c', 'y_f'] 
    109114ZNAME = [ 'z', 'Z', 'Z1', 'zz', 'ZZ', 'depth', 'tdepth', 'udepth', 
    110115              'vdepth', 'wdepth', 'fdepth', 'deptht', 'depthu', 
    111116              'depthv', 'depthw', 'depthf', 'olevel', 'z_c', 'z_f', ] 
    112 TNAME = [ 't', 'T', 'tt', 'TT', 'time', 'time_counter', 'time_centered', ] 
    113  
    114 ## All possibles name of units of dimensions in Nemo files 
     117TNAME = [ 't', 'T', 'tt', 'TT', 'time', 'time_counter', 'time_centered', 'time', 'TIME', 'TIME_COUNTER', 'TIME_CENTERED', ] 
     118 
     119## All possibles name of units of dimensions in NEMO files 
    115120XUNIT = [ 'degrees_east', ] 
    116121YUNIT = [ 'degrees_north', ] 
    … …  
    119124 
    120125## All possibles size of dimensions in Orca files 
    121 XLENGTH = [ 180, 182, 360, 362, 1440 ] 
    122 YLENGTH = [ 148, 149, 331, 332 ] 
    123 ZLENGTH = [ 31, 75] 
     126XLENGTH   = [ 180, 182, 360, 362, 1440, 1442 ] 
     127YLENGTH   = [ 148, 149, 331, 332 ] 
     128ZLENGTH   = [ 31, 75] 
     129XYZLENGTH = [ [180,148,31], [182,149,31], [360,331,75], [362,332,75] ] 
     130 
     131## T, S array to plot TS diagrams 
     132Ttab = np.linspace (-2, 40, 100) 
     133Stab = np.linspace ( 0, 40, 100) 
     134 
     135Ttab = xr.DataArray ( Ttab, dims=('Temperature',), coords=(Ttab,) ) 
     136Stab = xr.DataArray ( Stab, dims=('Salinity'   ,), coords=(Stab,) ) 
     137 
     138Ttab.attrs.update ( {'unit':'degrees_celcius', 'long_name':'Temperature'} ) 
     139Stab.attrs.update ( {'unit':'PSU',             'long_name':'Salinity'} ) 
     140 
     141# nemo internal options 
     142import warnings 
     143from typing import TYPE_CHECKING, Literal, TypedDict 
     144 
     145Stack = list() 
     146 
     147if TYPE_CHECKING : 
     148    Options = Literal [ "Debug", "Trace", "Depth", "Stack" ] 
     149 
     150    class T_Options (TypedDict) : 
     151        Debug = bool 
     152        Trace = bool 
     153        Depth = int 
     154        Stack = list() 
     155 
     156OPTIONS = { 'Debug':False, 'Trace':False, 'Depth':-1, 'Stack':list() } 
     157 
     158class set_options : 
     159    """ 
     160    Set options for nemo 
     161    """ 
     162    def __init__ (self, **kwargs): 
     163        self.old = {} 
     164        for k, v in kwargs.items(): 
     165            if k not in OPTIONS: 
     166                raise ValueError ( f"argument name {k!r} is not in the set of valid options {set(OPTIONS)!r}" ) 
     167            self.old[k] = OPTIONS[k] 
     168        self._apply_update(kwargs) 
     169 
     170    def _apply_update (self, options_dict) : OPTIONS.update (options_dict) 
     171    def __enter__ (self) : return 
     172    def __exit__ (self, type, value, traceback) : self._apply_update (self.old) 
     173 
     174def get_options () -> dict : 
     175    """ 
     176    Get options for nemo 
     177 
     178    See Also 
     179    ---------- 
     180    set_options 
     181 
     182    """ 
     183    return OPTIONS 
     184 
     185def return_stack () : 
     186    return Stack 
     187 
     188def PushStack (string:str) : 
     189    OPTIONS['Depth'] += 1 
     190    if OPTIONS['Trace'] : print ( '  '*OPTIONS['Depth'], '-->nemo:', string) 
     191    Stack.append (string) 
     192    return 
     193 
     194def PopStack (string:str) : 
     195    if OPTIONS['Trace'] : print ( '  '*OPTIONS['Depth'], '<--nemo:', string) 
     196    OPTIONS['Depth'] -= 1 
     197    Stack.pop () 
     198    return 
    124199 
    125200## =========================================================================== 
    … …  
    129204    Returns type 
    130205    ''' 
     206    PushStack ( f'__mmath__ ( ptab, {default=} )' ) 
     207     
    131208    mmath = default 
    132     if isinstance (ptab, xr.core.dataarray.DataArray) : 
    133         mmath = xr 
    134     if isinstance (ptab, np.ndarray) : 
    135         mmath = np 
    136     if isinstance (ptab, np.ma.MaskType) : 
    137         mmath = np.ma 
    138  
     209    if isinstance (ptab, xr.core.dataset.Dataset)     : mmath = 'dataset' 
     210    if isinstance (ptab, xr.core.dataarray.DataArray) : mmath = xr 
     211    if isinstance (ptab, np.ndarray)                  : mmath = np 
     212    if isinstance (ptab, np.ma.MaskType)              : mmath = np.ma 
     213 
     214    PopStack ( f'__math_ : {mmath}' ) 
    139215    return mmath 
    140216 
    141 def __guess_nperio__ (jpj, jpi, nperio=None, out='nperio') : 
     217def __guess_nperio__ (jpj:int, jpi:int, nperio=None, out:str='nperio') : 
    142218    '''Tries to guess the value of nperio (periodicity parameter. 
    143219 
    … …  
    148224    nperio : periodicity parameter 
    149225    ''' 
     226    PushStack ( f'__guess_nperio__ ( {jpj=} {jpi=} {nperio=}, {out=} )' ) 
    150227    if nperio is None : 
    151228        nperio = __guess_config__ (jpj, jpi, nperio=None, out=out) 
     229    PopStack (  f'__guess_nperio__ : {nperio}' ) 
    152230    return nperio 
    153231 
    154 def __guess_config__ (jpj, jpi, nperio=None, config=None, out='nperio') : 
     232def __guess_config__ (jpj:int, jpi:int, nperio:bool=None, config=None, out='nperio') : 
    155233    '''Tries to guess the value of nperio (periodicity parameter). 
    156234 
    … …  
    161239    nperio : periodicity parameter 
    162240    ''' 
    163     print ( jpi, jpj) 
     241    PushStack ( f'__guessConfig__ ( {jpj=}, {jpi=}, {nperio=}, {config=}, {out=}') 
     242     
     243    print (jpi, jpj) 
    164244    if nperio is None : 
    165245        ## Values for NEMO version < 4.2 
    … …  
    195275                                'nemo.py is not ready for this value' ) 
    196276 
    197     if out == 'nperio' : 
     277    if out == 'nperio'        : 
     278        PushStack ( f'__guessConfig__ : {nperio=}' ) 
    198279        return nperio 
    199     if out == 'config' : 
     280    if out == 'config'        : 
    200281        return config 
    201     if out == 'perio'  : 
     282        PushStack ( f'__guessConfig__ : {config=}' ) 
     283    if out == 'perio'         : 
    202284        return iperio, jperio, nfold, nftype 
     285        PushStack ( f'__guessConfig__ : {iperio=}, {jperio=}, {nfold=}, {nftype=}' ) 
    203286    if out in ['full', 'all'] : 
    204         return {'nperio':nperio, 'iperio':iperio, 'jperio':jperio, 'nfold':nfold, 'nftype':nftype} 
     287        zdict = {'nperio':nperio, 'iperio':iperio, 'jperio':jperio, 'nfold':nfold, 'nftype':nftype} 
     288        PushStack ( f'__guessConfig__ : {zdict}' ) 
     289        return zdict 
    205290 
    206291def __guess_point__ (ptab) : 
    … …  
    215300    Credits : who is the original author ? 
    216301    ''' 
     302    PushStack ( f'__guess_point__ ( ptab )') 
    217303 
    218304    gp = None 
    219305    mmath = __mmath__ (ptab) 
    220306    if mmath == xr : 
    221         if ('x_c' in ptab.dims and 'y_c' in ptab.dims ) : 
    222             gp = 'T' 
    223         if ('x_f' in ptab.dims and 'y_c' in ptab.dims ) : 
    224             gp = 'U' 
    225         if ('x_c' in ptab.dims and 'y_f' in ptab.dims ) : 
    226             gp = 'V' 
    227         if ('x_f' in ptab.dims and 'y_f' in ptab.dims ) : 
    228             gp = 'F' 
    229         if ('x_c' in ptab.dims and 'y_c' in ptab.dims 
    230               and 'z_c' in ptab.dims )                : 
    231             gp = 'T' 
    232         if ('x_c' in ptab.dims and 'y_c' in ptab.dims 
    233                 and 'z_f' in ptab.dims )                : 
    234             gp = 'W' 
    235         if ('x_f' in ptab.dims and 'y_c' in ptab.dims 
    236                 and 'z_f' in ptab.dims )                : 
    237             gp = 'U' 
    238         if ('x_c' in ptab.dims and 'y_f' in ptab.dims 
    239                 and 'z_f' in ptab.dims ) : 
    240             gp = 'V' 
    241         if ('x_f' in ptab.dims and 'y_f' in ptab.dims 
    242                 and 'z_f' in ptab.dims ) : 
    243             gp = 'F' 
     307        if ('x_c' in ptab.dims and 'y_c' in ptab.dims )                        : gp = 'T' 
     308        if ('x_f' in ptab.dims and 'y_c' in ptab.dims )                        : gp = 'U' 
     309        if ('x_c' in ptab.dims and 'y_f' in ptab.dims )                        : gp = 'V' 
     310        if ('x_f' in ptab.dims and 'y_f' in ptab.dims )                        : gp = 'F' 
     311        if ('x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_c' in ptab.dims ) : gp = 'T' 
     312        if ('x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims ) : gp = 'W' 
     313        if ('x_f' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims ) : gp = 'U' 
     314        if ('x_c' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims ) : gp = 'V' 
     315        if ('x_f' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims ) : gp = 'F' 
    244316 
    245317        if gp is None : 
    246318            raise AttributeError ('in nemo module : cd_type not found, and cannot by guessed') 
    247319        print ( f'Grid set as {gp} deduced from dims {ptab.dims}' ) 
     320 
     321        PopStack ( f'__guess_point__ : {gp}' ) 
    248322        return gp 
    249323    else : 
     324        PopStack ( f'__guess_point__ : cd_type not found, input is not an xarray data' ) 
    250325        raise AttributeError  ('in nemo module : cd_type not found, input is not an xarray data') 
    251326 
    252 def get_shape ( ptab ) : 
     327def get_shape (ptab) : 
    253328    '''Get shape of ptab return a string with axes names 
    254329 
    … …  
    258333    etc ... 
    259334    ''' 
    260  
     335    PushStack ( f'get_shape ( ptab )' ) 
    261336    g_shape = '' 
    262     if __find_axis__ (ptab, 'x')[0] : 
    263         g_shape = 'X' 
    264     if __find_axis__ (ptab, 'y')[0] : 
    265         g_shape = 'Y' + g_shape 
    266     if __find_axis__ (ptab, 'z')[0] : 
    267         g_shape = 'Z' + g_shape 
    268     if __find_axis__ (ptab, 't')[0] : 
    269         g_shape = 'T' + g_shape 
     337    if __find_axis__ (ptab, 'x')[0] : g_shape = 'X' 
     338    if __find_axis__ (ptab, 'y')[0] : g_shape = 'Y' + g_shape 
     339    if __find_axis__ (ptab, 'z')[0] : g_shape = 'Z' + g_shape 
     340    if __find_axis__ (ptab, 't')[0] : g_shape = 'T' + g_shape 
     341         
     342    PopStack ( f'get_shape : {g_shape}' ) 
    270343    return g_shape 
    271344 
    272 def lbc_diag (nperio) : 
     345def lbc_diag (nperio:int) : 
    273346    '''Useful to switch between field with and without halo''' 
     347    PushStack ( f'lbc_diag ( {nperio=}' ) 
    274348    lperio, aperio = nperio, False 
    275     if nperio == 4.2 : 
    276         lperio, aperio = 4, True 
    277     if nperio == 6.2 : 
    278         lperio, aperio = 6, True 
     349    if nperio == 4.2 : lperio, aperio = 4, True 
     350    if nperio == 6.2 : lperio, aperio = 6, True 
     351 
     352    PopStack ( f'lbc_diag : {lperio}, {aperio}' ) 
    279353    return lperio, aperio 
    280354 
    281 def __find_axis__ (ptab, axis='z', back=True) : 
     355def __find_axis__ (ptab, axis='z', back:bool=True, verbose:bool=False) : 
    282356    '''Returns name and name of the requested axis''' 
     357    PushStack ( '__find_axis__ ( ptab, {axis=}, {back=}, {verbose=}' ) 
     358     
    283359    mmath = __mmath__ (ptab) 
    284360    ax, ix = None, None 
    … …  
    286362    if axis in XNAME : 
    287363        ax_name, unit_list, length = XNAME, XUNIT, XLENGTH 
     364        if verbose : print ( f'Working on xaxis found by name : {axis=} : {XNAME=} {ax_name=} {unit_list=} {length=}' ) 
    288365    if axis in YNAME : 
    289366        ax_name, unit_list, length = YNAME, YUNIT, YLENGTH 
     367        if verbose : print ( f'Working on yaxis found by name : {axis=} : {YNAME=} {ax_name=} {unit_list=} {length=}' ) 
    290368    if axis in ZNAME : 
    291369        ax_name, unit_list, length = ZNAME, ZUNIT, ZLENGTH 
     370        if verbose : print ( f'Working on zaxis found by name : {axis=} : {ZNAME=} {ax_name=} {unit_list=} {length=}' ) 
    292371    if axis in TNAME : 
    293372        ax_name, unit_list, length = TNAME, TUNIT, None 
    294  
    295     if mmath == xr : 
     373        if verbose : print ( f'Working on taxis found by name : {axis=} : {TNAME=} {ax_name=} {unit_list=} {length=}' ) 
     374 
     375    if verbose : 
     376        print ( f'{ax_name=} {unit_list=} {length=}' ) 
     377             
     378    if mmath in [xr, 'dataset' ]:  
    296379        # Try by name 
    297         for dim in ax_name : 
    298             if dim in ptab.dims : 
    299                 ix, ax = ptab.dims.index (dim), dim 
     380        for i, dim in enumerate (ptab.dims) : 
     381            if verbose : print ( f'{i=} {dim=}' ) 
     382            if dim in ax_name : 
     383                if verbose : print ( f'Rule 2 : {dim=} axis found by name : {XNAME=}' ) 
     384                ix, ax = i, dim 
    300385 
    301386        # If not found, try by axis attributes 
    302387        if not ix : 
     388            if verbose : print ( 'ix not found - 1' ) 
    303389            for i, dim in enumerate (ptab.dims) : 
     390                if verbose : print ( f'{i=} {dim=}' ) 
    304391                if 'axis' in ptab.coords[dim].attrs.keys() : 
    305392                    l_axis = ptab.coords[dim].attrs['axis'] 
    306                     if axis in ax_name and l_axis == 'X' : 
     393                    if verbose : print ( f'Rule 3 : Trying {i=} {dim=} {l_axis=}' ) 
     394                    if l_axis in ax_name and l_axis == 'X' : 
     395                        if verbose : print ( f'Rule 3 : xaxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' ) 
    307396                        ix, ax = (i, dim) 
    308                     if axis in ax_name and l_axis == 'Y' : 
     397                    if l_axis in ax_name and l_axis == 'Y' : 
     398                        if verbose : print ( f'Rule 3 : yaxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' ) 
    309399                        ix, ax = (i, dim) 
    310                     if axis in ax_name and l_axis == 'Z' : 
     400                    if l_axis in ax_name and l_axis == 'Z' : 
     401                        if verbose : print ( f'Rule 3 : zaxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' ) 
    311402                        ix, ax = (i, dim) 
    312                     if axis in ax_name and l_axis == 'T' : 
     403                    if l_axis in ax_name and l_axis == 'T' : 
     404                        if verbose : print ( f'Rule 3 : taxis found by name : {ax=} {l_axis=} {axis=} : {ax_name=} {l_axis=} {i=} {dim=}' ) 
    313405                        ix, ax = (i, dim) 
    314406 
    … …  
    319411                    for name in unit_list : 
    320412                        if name in ptab.coords[dim].attrs['units'] : 
     413                            if verbose : print ( f'Rule 4 : {name=} found by unit : {axis=} : {unit_list=} {i=} {dim=}' ) 
    321414                            ix, ax = i, dim 
    322415 
    323416    # If numpy array or dimension not found, try by length 
    324     if mmath != xr or not ix : 
     417     
     418    if mmath not in [xr, 'dataset'] or not ix : 
    325419        if length : 
    326             l_shape = ptab.shape 
     420            if mmath in [xr, 'dataset'] : 
     421                l_shape =[ len (x) for x in ptab.dims ] 
     422            else : 
     423                l_shape = ptab.shape 
    327424            for nn in np.arange ( len(l_shape) ) : 
    328425                if l_shape[nn] in length : 
    329                     ix = nn 
     426                    ix = nn ; ax = None 
     427                    if verbose : print ( f'Rule 5 : {ax_name=} axis found by length : {axis=} : {XNAME=} {ix=} {ax=}' ) 
    330428 
    331429    if ix and back : 
    332         ix -= len(ptab.shape) 
    333  
     430        if mmath in [xr, 'dataset'] : 
     431            ix -= len (ptab.dims) 
     432        else : 
     433            ix -= len (ptab.shape) 
     434 
     435    PopStack ( f'__find_axis__ : {ax}, {ix}' ) 
    334436    return ax, ix 
    335437 
    336 def find_axis ( ptab, axis='z', back=True ) : 
     438def find_axis (ptab, axis:str='z', back:bool=True, verbose:bool=False) : 
    337439    '''Version of find_axis with no __''' 
    338     ix, xx = __find_axis__ (ptab, axis, back) 
     440    PushStack ( f'find_axis ( ptab, {axis=}, {back=}, {varbose=}' ) 
     441     
     442    ix, xx = __find_axis__ (ptab, axis, back, verbose) 
     443     
     444    PopStack ( f'find_axis : {xx}, {ix}' ) 
    339445    return xx, ix 
    340446 
    341 def fixed_lon (plon, center_lon=0.0) : 
     447def fixed_lon (plon, center_lon:float=0.0) : 
    342448    '''Returns corrected longitudes for nicer plots 
    343449 
    … …  
    348454    See https://gist.github.com/pelson/79cf31ef324774c97ae7 
    349455    ''' 
     456    PushStack ( f'fixed_lon ( plon, {center_lon=}' ) 
    350457    mmath = __mmath__ (plon) 
    351458 
    … …  
    359466 
    360467    # Special case for eORCA025 
    361     if f_lon.shape [-1] == 1442 : 
    362         f_lon [..., -2, :] = f_lon [..., -3, :] 
    363     if f_lon.shape [-1] == 1440 : 
    364         f_lon [..., -1, :] = f_lon [..., -2, :] 
    365  
    366     if f_lon.min () > center_lon : 
    367         f_lon += -360.0 
    368     if f_lon.max () < center_lon : 
    369         f_lon +=  360.0 
    370  
    371     if f_lon.min () < center_lon-360.0 : 
    372         f_lon +=  360.0 
    373     if f_lon.max () > center_lon+360.0 : 
    374         f_lon += -360.0 
    375  
     468    if f_lon.shape [-1] == 1442 : f_lon [..., -2, :] = f_lon [..., -3, :] 
     469    if f_lon.shape [-1] == 1440 : f_lon [..., -1, :] = f_lon [..., -2, :] 
     470 
     471    if f_lon.min () > center_lon       : f_lon += -360.0 
     472    if f_lon.max () < center_lon       : f_lon +=  360.0 
     473 
     474    if f_lon.min () < center_lon-360.0 : f_lon +=  360.0 
     475    if f_lon.max () > center_lon+360.0 : f_lon += -360.0 
     476 
     477    PopStack ( f'fixed_lon : f_lon' ) 
    376478    return f_lon 
    377479 
    378 def bounds_clolon ( pbounds_lon, plon, rad=False, deg=True) : 
     480def bounds_clolon (pbounds_lon, plon, rad:bool=False, deg:bool=True) : 
    379481    '''Choose closest to lon0 longitude, adding/substacting 360° if needed 
    380482    ''' 
    381  
    382     if rad : 
    383         lon_range = 2.0*np.pi 
    384     if deg : 
    385         lon_range = 360.0 
     483    PushStack ( f'bounds_clolon ( pbounds_lon, plon, {rad=}, {deb=}' ) 
     484    if rad : lon_range = 2.0*np.pi 
     485    if deg : lon_range = 360.0 
    386486    b_clolon = pbounds_lon.copy () 
    387487 
    … …  
    392492                          b_clolon-lon_range, 
    393493                          b_clolon ) 
     494 
     495    PopStack ( f'bounds_clolon : b_clolon' ) 
    394496    return b_clolon 
    395497 
    396 def unify_dims ( dd, x='x', y='y', z='olevel', t='time_counter', verbose=False ) : 
     498def unify_dims ( dd: "DataArray or Dataset", x=UDIMS['x'], y=UDIMS['y'], z=UDIMS['z'], t=UDIMS['t'], verbose:bool=False ) -> "DataArray or Dataset" : 
    397499    '''Rename dimensions to unify them between NEMO versions 
    398500    ''' 
     501    PushStack ( f'unify_dims ( dd, {x=}, {y=}, {z=}, {t=}, {verbose=}' ) 
    399502    for xx in XNAME : 
    400503        if xx in dd.dims and xx != x : 
    401             if verbose : 
    402                 print ( f"{xx} renamed to {x}" ) 
     504            if verbose : print ( f"{xx} renamed to {x}" ) 
    403505            dd = dd.rename ( {xx:x}) 
    404506 
    405507    for yy in YNAME : 
    406508        if yy in dd.dims and yy != y  : 
    407             if verbose : 
    408                 print ( f"{yy} renamed to {y}" ) 
     509            if verbose : print ( f"{yy} renamed to {y}" ) 
    409510            dd = dd.rename ( {yy:y} ) 
    410511 
    411512    for zz in ZNAME : 
    412513        if zz in dd.dims and zz != z : 
    413             if verbose : 
    414                 print ( f"{zz} renamed to {z}" ) 
     514            if verbose : print ( f"{zz} renamed to {z}" ) 
    415515            dd = dd.rename ( {zz:z} ) 
    416516 
    417517    for tt in TNAME  : 
    418518        if tt in dd.dims and tt != t : 
    419             if verbose : 
    420                 print ( f"{tt} renamed to {t}" ) 
     519            if verbose : print ( f"{tt} renamed to {t}" ) 
    421520            dd = dd.rename ( {tt:t} ) 
    422521 
     522    PopStack ( f'unify_dims : dd' ) 
    423523    return dd 
    424524 
    425525 
    426526if SimpleImputer :  
    427   def fill_empty (ptab, sval=np.nan, transpose=False) : 
     527    def fill_empty (ptab, sval:float=np.nan, transpose:bool=False) : 
    428528        '''Fill empty values 
    429529 
    … …  
    432532        values 
    433533        ''' 
     534        PushStack ( f'fill_empty ( ptab, {sval=} {transpose=}' ) 
    434535        mmath = __mmath__ (ptab) 
    435536 
    … …  
    446547            ztab.attrs.update (ptab.attrs) 
    447548 
     549        PopStack ( f'fill_empty : ztab' ) 
    448550        return ztab 
    449551 
     552    def fill_latlon (plat, plon, sval=-1) : 
     553        '''Fill latitude/longitude values 
     554 
     555        Useful when NEMO has run with no wet points options : 
     556        some parts of the domain, with no ocean points, have no 
     557        lon/lat values 
     558        ''' 
     559        PushStack ( f'fill_latlon ( plat, plon {sval=}' ) 
     560 
     561        from sklearn.impute import SimpleImputer 
     562        mmath = __mmath__ (plon) 
     563         
     564        imp = SimpleImputer (missing_values=sval, strategy='mean') 
     565        imp.fit (plon) 
     566        zlon = imp.transform (plon) 
     567        imp.fit (plat.T) 
     568        zlat = imp.transform (plat.T).T 
     569         
     570        if mmath == xr : 
     571            zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords) 
     572            zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords) 
     573            zlon.attrs.update (plon.attrs) 
     574            zlat.attrs.update (plat.attrs) 
     575             
     576        zlon = fixed_lon (zlon) 
     577 
     578        PopStack ( f'fixed_lon') 
     579        return zlat, zlon 
     580     
     581    def fill_lonlat (plon, plat, sval=-1) : 
     582        '''Fill longitude/latitude values 
     583         
     584        Useful when NEMO has run with no wet points options : 
     585        some parts of the domain, with no ocean points, have no 
     586        lon/lat values 
     587        ''' 
     588        PushStack ( f'fill_lonlat ( plon, plat {sval=}' ) 
     589        from sklearn.impute import SimpleImputer 
     590        mmath = __mmath__ (plon) 
     591         
     592        imp = SimpleImputer (missing_values=sval, strategy='mean') 
     593        imp.fit (plon) 
     594        zlon = imp.transform (plon) 
     595        imp.fit (plat.T) 
     596        zlat = imp.transform (plat.T).T 
     597         
     598        if mmath == xr : 
     599            zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords) 
     600            zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords) 
     601            zlon.attrs.update (plon.attrs) 
     602            zlat.attrs.update (plat.attrs) 
     603             
     604        zlon = fixed_lon (zlon) 
     605 
     606        PopStack ( f'fill_lonlat' ) 
     607        return zlon, zlat 
     608 
     609    def fill_bounds_lonlat (pbounds_lon, pbounds_lat, sval=-1) : 
     610        '''Fill longitude/latitude bounds values 
     611         
     612        Useful when NEMO has run with no wet points options : 
     613        some parts of the domain, with no ocean points, as no 
     614        lon/lat values 
     615        ''' 
     616        PushStack ( f'fill_bounds_lonlat (pbounds_lon, pbounds_lat, {sval=}' ) 
     617        mmath = __mmath__ (pbounds_lon) 
     618         
     619        z_bounds_lon = np.empty_like (pbounds_lon) 
     620        z_bounds_lat = np.empty_like (pbounds_lat) 
     621         
     622        imp = SimpleImputer (missing_values=sval, strategy='mean') 
     623         
     624        for n in np.arange (4) : 
     625            imp.fit (pbounds_lon[:,:,n]) 
     626            z_bounds_lon[:,:,n] = imp.transform (pbounds_lon[:,:,n]) 
     627            imp.fit (pbounds_lat[:,:,n].T) 
     628            z_bounds_lat[:,:,n] = imp.transform (pbounds_lat[:,:,n].T).T 
     629             
     630        if mmath == xr : 
     631            z_bounds_lon = xr.DataArray (pbounds_lon, dims=pbounds_lon.dims, 
     632                                            coords=pbounds_lon.coords) 
     633            z_bounds_lat = xr.DataArray (pbounds_lat, dims=pbounds_lat.dims, 
     634                                            coords=pbounds_lat.coords) 
     635            z_bounds_lon.attrs.update (pbounds_lat.attrs) 
     636            z_bounds_lat.attrs.update (pbounds_lat.attrs) 
     637 
     638        PopStack ( 'fill_bounds_lonlat' ) 
     639        return z_bounds_lon, z_bounds_lat 
    450640     
    451641else :  
    … …  
    461651          values 
    462652        ''' 
     653        PushStack ( f'fill_empty [void version] ( ptab, {sval=} {transpose=}' ) 
     654 
    463655        print ( 'Error : module sklearn.impute.SimpleImputer not found' ) 
    464656        print ( 'Can not call fill_empty' ) 
    465657        print ( 'Call arguments where : ' ) 
    466658        print ( f'{ptab.shape=} {sval=} {transpose=}' ) 
    467    
    468 def fill_lonlat (plon, plat, sval=-1) : 
    469     '''Fill longitude/latitude values 
    470  
    471     Useful when NEMO has run with no wet points options : 
    472     some parts of the domain, with no ocean points, have no 
    473     lon/lat values 
    474     ''' 
    475     from sklearn.impute import SimpleImputer 
    476     mmath = __mmath__ (plon) 
    477  
    478     imp = SimpleImputer (missing_values=sval, strategy='mean') 
    479     imp.fit (plon) 
    480     zlon = imp.transform (plon) 
    481     imp.fit (plat.T) 
    482     zlat = imp.transform (plat.T).T 
    483  
    484     if mmath == xr : 
    485         zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords) 
    486         zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords) 
    487         zlon.attrs.update (plon.attrs) 
    488         zlat.attrs.update (plat.attrs) 
    489  
    490     zlon = fixed_lon (zlon) 
    491  
    492     return zlon, zlat 
    493  
    494 def fill_bounds_lonlat (pbounds_lon, pbounds_lat, sval=-1) : 
    495     '''Fill longitude/latitude bounds values 
    496  
    497     Useful when NEMO has run with no wet points options : 
    498     some parts of the domain, with no ocean points, as no 
    499     lon/lat values 
    500     ''' 
    501     mmath = __mmath__ (pbounds_lon) 
    502  
    503     z_bounds_lon = np.empty ( pbounds_lon.shape ) 
    504     z_bounds_lat = np.empty ( pbounds_lat.shape ) 
    505  
    506     imp = SimpleImputer (missing_values=sval, strategy='mean') 
    507  
    508     for n in np.arange (4) : 
    509         imp.fit (pbounds_lon[:,:,n]) 
    510         z_bounds_lon[:,:,n] = imp.transform (pbounds_lon[:,:,n]) 
    511         imp.fit (pbounds_lat[:,:,n].T) 
    512         z_bounds_lat[:,:,n] = imp.transform (pbounds_lat[:,:,n].T).T 
    513  
    514     if mmath == xr : 
    515         z_bounds_lon = xr.DataArray (pbounds_lon, dims=pbounds_lon.dims, 
    516                                          coords=pbounds_lon.coords) 
    517         z_bounds_lat = xr.DataArray (pbounds_lat, dims=pbounds_lat.dims, 
    518                                          coords=pbounds_lat.coords) 
    519         z_bounds_lon.attrs.update (pbounds_lat.attrs) 
    520         z_bounds_lat.attrs.update (pbounds_lat.attrs) 
    521  
    522     return z_bounds_lon, z_bounds_lat 
     659        PopStack ( f'fill_empty [void version]' ) 
     660 
     661 
     662    def fill_latlon (plat, plon, sval=-1) : 
     663        ''' Void version of fill_latlon, because module sklearn.impute.SimpleImputer is not available 
     664 
     665        Useful when NEMO has run with no wet points options : 
     666        some parts of the domain, with no ocean points, have no 
     667        lon/lat values 
     668        ''' 
     669        PushStack ( f'fill_latlon [void_version] ( plat, plon {sval=}' ) 
     670 
     671        print ( 'Error : module sklearn.impute.SimpleImputer not found' ) 
     672        print ( 'Can not call fill_empty' ) 
     673        print ( 'Call arguments where : ' ) 
     674        print ( f'{ptab.shape=} {sval=} {transpose=}' ) 
     675 
     676        PopStack ( f'fixed_lon [void version]') 
     677         
     678    def fill_lonlat (plon, plat, sval=-1) : 
     679        ''''Void version of fill_lonlat, because module sklearn.impute.SimpleImputer is not available 
     680         
     681        Useful when NEMO has run with no wet points options : 
     682        some parts of the domain, with no ocean points, have no 
     683        lon/lat values 
     684        ''' 
     685        PushStack ( f'fill_lonlat [void Version] ( plon, plat {sval=}' ) 
     686 
     687        print ( 'Error : module sklearn.impute.SimpleImputer not found' ) 
     688        print ( 'Can not call fill_empty' ) 
     689        print ( 'Call arguments where : ' ) 
     690        print ( f'{ptab.shape=} {sval=} {transpose=}' ) 
     691 
     692        PopStack ( f'fill_lonlat [void version]' ) 
     693     
     694    def fill_bounds_lonlat (pbounds_lon, pbounds_lat, sval=-1) : 
     695        ''''Void version of fill_bounds_lonlat, because module sklearn.impute.SimpleImputer is not available 
     696         
     697        Useful when NEMO has run with no wet points options : 
     698        some parts of the domain, with no ocean points, as no 
     699        lon/lat values 
     700        ''' 
     701        PushStack ( f'fill_bounds_lonlat [void version] (pbounds_lon, pbounds_lat, {sval=}' ) 
     702 
     703        print ( 'Error : module sklearn.impute.SimpleImputer not found' ) 
     704        print ( 'Can not call fill_empty' ) 
     705        print ( 'Call arguments where : ' ) 
     706        print ( f'{ptab.shape=} {sval=} {transpose=}' ) 
     707 
     708        PopStack ( 'fill_bounds_lonlat [void version]' ) 
    523709 
    524710def jeq (plat) : 
    … …  
    527713    lat : latitudes of the grid. At least 2D. 
    528714    ''' 
     715    PushStack ( f'jeq ( plat ) ' ) 
    529716    mmath = __mmath__ (plat) 
    530717    jy = __find_axis__ (plat, 'y')[-1] 
    531718 
    532     if mmath == xr : 
    533         jj = int ( np.mean ( np.argmin (np.abs (np.float64 (plat)), 
    534                                                 axis=jy) ) ) 
    535     else : 
    536         jj = np.argmin (np.abs (np.float64 (plat[...,:, 0]))) 
    537  
     719    if mmath == xr : jj = int ( np.mean ( np.argmin (np.abs (np.float64 (plat)), axis=jy) ) ) 
     720    else           : jj = np.argmin (np.abs (np.float64 (plat[...,:, 0]))) 
     721 
     722    PopStack ( f'jeq : {jj}' ) 
    538723    return jj 
    539724 
    … …  
    544729    plat (optionnal) : latitudes of the grid 
    545730    ''' 
     731    PushStack ( f'lon1d ( plon, plat) ' )  
    546732    mmath = __mmath__ (plon) 
    547733    jpj, jpi  = plon.shape [-2:] 
    … …  
    566752        lon_1d.attrs['long_name :']   = 'Longitude' 
    567753 
     754    PopStack ( 'lon_1d' ) 
    568755    return lon_1d 
    569756 
    … …  
    575762    diff [optional] : tolerance 
    576763    ''' 
     764    PushStack ( f'latreg ( plat, {diff=}' ) 
    577765    #mmath = __mmath__ (plat) 
    578766    if diff is None : 
    … …  
    587775    lareg  = np.float64 (plat[...,jreg,:].mean(axis=-1)) 
    588776 
     777    PopStack ( f'latreg : {jreg}, lareg' ) 
    589778    return jreg, lareg 
    590779 
    … …  
    594783    plat : latitudes of the grid (2D) 
    595784    ''' 
     785    PushStack ( f'lat1d ( plat )' ) 
    596786    mmath = __mmath__ (plat) 
    597787    iy = __find_axis__ (plat, 'y')[-1] 
    … …  
    606796 
    607797    if np.abs (lat_eq) < dy/100. : # T, U or W grid 
    608         if jpj-1 > jreg : 
    609             dys = (90.-lat_reg) / (jpj-jreg-1)*0.5 
    610         else            : 
    611             dys = (90.-lat_reg) / 2.0 
     798        if jpj-1 > jreg : dys = (90.-lat_reg) / (jpj-jreg-1)*0.5 
     799        else            : dys = (90.-lat_reg) / 2.0 
    612800        yrange = 90.-dys-lat_reg 
    613     else                           :  # V or F grid 
    614         yrange = 90.-lat_reg 
     801    else :  yrange = 90.-lat_reg # V or F grid 
    615802 
    616803    if jpj-1 > jreg : 
    … …  
    629816        lat_1d.attrs ['long_name :']   = 'Latitude' 
    630817 
     818    PopStack ( f'lat1d' ) 
    631819    return lat_1d 
    632820 
    … …  
    636824    plat, plon : latitudes and longitudes of the grid (2D) 
    637825    ''' 
    638     return lat1d (plat),  lon1d (plon, plat) 
     826    PushStack ( f'latlon1d ( plat, plon) ' ) 
     827    zla = lat1d (plat) 
     828    zlo = lon1d (plon, plat) 
     829 
     830    PopStack ( f'latlon1d' ) 
     831    return zla, zlo 
    639832 
    640833def ff (plat) : 
    641834    '''Returns Coriolis factor 
    642835    ''' 
    643     zff   = np.sin (RAD * plat) * OMEGA 
     836    PushStack ( 'ff ( plat )' ) 
     837    zff   = np.sin (RAD * plat) * ROMEGA 
     838    PopStack ( 'ff' ) 
    644839    return zff 
    645840 
    … …  
    647842    '''Return Beta factor (derivative of Coriolis factor) 
    648843    ''' 
    649     zbeta = np.cos (RAD * plat) * OMEGA / RA 
     844    PushStack ( 'beta ( plat )' ) 
     845    zbeta = np.cos (RAD * plat) * ROMEGA / RA 
     846    PopStack ( 'beta' ) 
    650847    return zbeta 
    651848 
    … …  
    653850    '''Returns masked values outside a lat/lon box 
    654851    ''' 
     852    PushStack ( f'mask_lonlat  (ptab, {x0=}, {x1=}, {y0=}, {y1=}, lon, lat, {sval=}' ) 
    655853    mmath = __mmath__ (ptab) 
    656854    if mmath == xr : 
    … …  
    660858    mask = np.logical_and (np.logical_and(lat>y0, lat<y1), 
    661859            np.logical_or (np.logical_or ( 
    662                 np.logical_and(lon>x0, lon<x1), 
    663                 np.logical_and(lon+360>x0, lon+360<x1)), 
    664                 np.logical_and(lon-360>x0, lon-360<x1))) 
     860                np.logical_and (lon>x0, lon<x1), 
     861                np.logical_and (lon+360>x0, lon+360<x1)), 
     862                np.logical_and (lon-360>x0, lon-360<x1))) 
    665863    tab = mmath.where (mask, ptab, sval) 
    666864 
     865    PopStack ( 'mask_lonlat' ) 
    667866    return tab 
    668867 
    … …  
    684883 
    685884    ''' 
     885    PushStack ( f'extend ( ptab, {blon=}, {jplus=}, {jpi=}, {nperio=}' ) 
    686886    mmath = __mmath__ (ptab) 
    687887 
    688     if ptab.shape[-1] == 1 : 
    689         tabex = ptab 
     888    if ptab.shape[-1] == 1 : tabex = ptab 
    690889 
    691890    else : 
    692         if jpi is None : 
    693             jpi = ptab.shape[-1] 
    694  
    695         if blon : 
    696             xplus = -360.0 
    697         else   : 
    698             xplus =    0.0 
     891        if not jpi : jpi = ptab.shape[-1] 
     892 
     893        if blon : xplus = -360.0 
     894        else    : xplus =    0.0 
    699895 
    700896        if ptab.shape[-1] > jpi : 
    701897            tabex = ptab 
    702898        else : 
    703             if nperio in [ 0, 4.2 ] : 
    704                 istart, le, la = 0, jpi+1, 0 
    705             if nperio == 1 : 
    706                 istart, le, la = 0, jpi+1, 0 
     899            if nperio in [ 0, 4.2 ] : istart, le, la = 0, jpi+1, 0 
     900            if nperio == 1          : istart, le, la = 0, jpi+1, 0 
    707901            if nperio in [4, 6] : # OPA case with two halo points for periodicity 
    708902                # Perfect, except at the pole that should be masked by lbc_plot 
    … …  
    716910                new_coords = [] 
    717911                for coord in ptab.dims : 
    718                     if coord == lon : 
    719                         new_coords.append ( np.arange( tabex.shape[-1])) 
    720                     else            : 
    721                         new_coords.append ( ptab.coords[coord].values) 
     912                    if coord == lon : new_coords.append ( np.arange( tabex.shape[-1])) 
     913                    else            : new_coords.append ( ptab.coords[coord].values) 
    722914                tabex = xr.DataArray ( tabex, dims=ptab.dims, 
    723915                                           coords=new_coords ) 
    … …  
    727919                     ptab [..., istart+la:istart+la+jplus]          ), 
    728920                     axis=-1) 
     921 
     922    PopStack ( 'extend' ) 
    729923    return tabex 
    730924 
    731 def orca2reg (dd, lat_name='nav_lat', lon_name='nav_lon', 
    732                   y_name='y', x_name='x') : 
     925def orca2reg (dd, lat_name=None, lon_name=None, y_name=None, x_name=None) : 
    733926    '''Assign an ORCA dataset on a regular grid. 
    734927 
    … …  
    741934      Returns : xarray dataset with rectangular grid. Incorrect above 20°N 
    742935    ''' 
     936    PushStack ( f'orca2reg ( dd, {lat_name=}, {lon_name=}, {y_name=}, {x_name=}' ) 
     937    if not x_name : x_name, ix = __find_axis__ (dd, axis='x') 
     938    if not y_name : y_name, jy = __find_axis__ (dd, axis='y') 
     939 
     940    if not lon_name : 
     941        for xn in LONNAME : 
     942            if xn in dd.variables : lon_name = xn 
     943    if not lat_name : 
     944        for yn in LATNAME : 
     945            if yn in dd.variables : lat_name = yn 
     946 
     947    print ( f'{dd.dims=}'   ) 
     948    print ( f'{x_name=} {y_name=}' ) 
     949     
    743950    # Compute 1D longitude and latitude 
    744     (zlat, zlon) = latlon1d ( dd[lat_name], dd[lon_name]) 
    745  
     951    ylat, xlon = fill_latlon (dd[lat_name], dd[lon_name], sval=-1) 
     952    (zlat, zlon) = latlon1d (ylat, xlon) 
    746953    zdd = dd 
     954     
    747955    # Assign lon and lat as dimensions of the dataset 
    748956    if y_name in zdd.dims : 
    … …  
    754962    # Force dimensions to be in the right order 
    755963    coord_order = ['lat', 'lon'] 
    756     for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z', 
     964    for dim in [ 'olevel', 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z', 
    757965                 'time_counter', 'time', 'tbnds', 
    758                  'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] : 
    759         if dim in zdd.dims : 
    760             coord_order.insert (0, dim) 
     966                 'bnds', 'axis_nbounds', 'nvertex', 'two2', 'two1', 'two', 'four',] : 
     967        if dim in zdd.dims : coord_order.insert (0, dim) 
    761968 
    762969    zdd = zdd.transpose (*coord_order) 
     970 
     971    PopStack ( 'orca2reg' ) 
    763972    return zdd 
    764973 
    … …  
    777986    See NEMO documentation for further details 
    778987    ''' 
     988    PushStack ( f'lbc_init ( ptab, {nperio=}' ) 
    779989    jpi, jpj = None, None 
    780990    ax, ix = __find_axis__ (ptab, 'x') 
    781991    ay, jy = __find_axis__ (ptab, 'y') 
    782     if ax : 
    783         jpi = ptab.shape[ix] 
    784     if ay : 
    785         jpj = ptab.shape[jy] 
    786  
    787     if nperio is None : 
    788         nperio = __guess_nperio__ (jpj, jpi, nperio) 
     992    if ax : jpi = ptab.shape[ix] 
     993    if ay : jpj = ptab.shape[jy] 
     994 
     995    if nperio is None : nperio = __guess_nperio__ (jpj, jpi, nperio) 
    789996 
    790997    if nperio not in NPERIO_VALID_RANGE : 
    791998        raise AttributeError ( f'{nperio=} is not in the valid range {NPERIO_VALID_RANGE}' ) 
    792999 
     1000    PopStack ( f'lbc_init : {jpj}, {jpi}, {nperio}' ) 
    7931001    return jpj, jpi, nperio 
    7941002 
    … …  
    8031011    See NEMO documentation for further details 
    8041012    ''' 
     1013    PushStack ( f'lbc_init ( ptab, {nperio=}, {cd_type=}, {psgn=}, {nemo_4u_bug=}' ) 
     1014 
    8051015    jpi, nperio = lbc_init (ptab, nperio)[1:] 
    806     ax = __find_axis__ (ptab, 'x')[0] 
    807     ay = __find_axis__ (ptab, 'y')[0] 
     1016    ax     = __find_axis__ (ptab, 'x')[0] 
     1017    ay     = __find_axis__ (ptab, 'y')[0] 
    8081018    psgn   = ptab.dtype.type (psgn) 
    8091019    mmath  = __mmath__ (ptab) 
    8101020 
    811     if mmath == xr : 
    812         ztab = ptab.values.copy () 
    813     else           : 
    814         ztab = ptab.copy () 
     1021    if mmath == xr : ztab = ptab.values.copy () 
     1022    else           : ztab = ptab.copy () 
    8151023 
    8161024    if ax : 
    … …  
    8581066                if cd_type in [ 'T', 'W' ] : # T-, W-point 
    8591067                    ztab [..., -1, jpi//2:  ] = psgn * ztab [..., -1, jpi//2:0:-1  ] 
    860  
    861                 if cd_type == 'U' : 
    862                     ztab [..., -1, jpi//2-1:-1] = psgn * ztab [..., -1, jpi//2:0:-1] 
    863  
    864                 if cd_type == 'V' : 
    865                     ztab [..., -1, 1:       ] = psgn * ztab [..., -2, jpi-1:0:-1   ] 
    866  
    867                 if cd_type == 'F' : 
    868                     ztab [..., -1, 0:-1     ] = psgn * ztab [..., -2, -1:0:-1      ] 
     1068                if cd_type == 'U' : ztab [..., -1, jpi//2-1:-1] = psgn * ztab [..., -1, jpi//2:0:-1] 
     1069                if cd_type == 'V' : ztab [..., -1, 1:       ] = psgn * ztab [..., -2, jpi-1:0:-1   ] 
     1070                if cd_type == 'F' : ztab [..., -1, 0:-1     ] = psgn * ztab [..., -2, -1:0:-1      ] 
    8691071 
    8701072            if nperio in [5, 6] :            #  North fold F-point pivot 
    … …  
    8971099        ztab.attrs = ptab.attrs 
    8981100 
     1101    PopStack ( 'lbc' ) 
    8991102    return ztab 
    9001103 
    … …  
    9081111    See NEMO documentation for further details 
    9091112    ''' 
     1113    PushStack ( f'lbc_mask (ptab, {nperio=}, {cd_type=}, {sval=}' ) 
    9101114    jpi, nperio = lbc_init (ptab, nperio)[1:] 
    9111115    ax = __find_axis__ (ptab, 'x')[0] 
    … …  
    9261130            #> South (in which nperio cases ?) 
    9271131            # -------------------------------- 
    928             if nperio in [1, 3, 4, 5, 6] : 
    929                 ztab [..., 0, :] = sval 
     1132            if nperio in [1, 3, 4, 5, 6] : ztab [..., 0, :] = sval 
    9301133 
    9311134            # 
    … …  
    9521155                if cd_type in [ 'T', 'W' ] : # T-, W-point 
    9531156                    ztab [..., -1, jpi//2  :  ] = sval 
    954  
    955                 if cd_type == 'U' : 
    956                     ztab [..., -1, jpi//2-1:-1] = sval 
    957  
    958                 if cd_type == 'V' : 
    959                     ztab [..., -1, 1:       ] = sval 
    960  
    961                 if cd_type == 'F' : 
    962                     ztab [..., -1, 0:-1     ] = sval 
     1157                if cd_type == 'U' : ztab [..., -1, jpi//2-1:-1] = sval 
     1158                if cd_type == 'V' : ztab [..., -1, 1:         ] = sval 
     1159                if cd_type == 'F' : ztab [..., -1, 0:-1       ] = sval 
    9631160 
    9641161            if nperio in [5, 6] :            #  North fold F-point pivot 
    … …  
    9791176                    ztab [..., -2, jpi//2+1:-1] = sval 
    9801177 
     1178    PopStack ( 'lbc_mask' ) 
    9811179    return ztab 
    9821180 
    … …  
    9941192    See NEMO documentation for further details 
    9951193    ''' 
     1194    PushStack ( f'lbc_plot (ptab, {nperio=}, {cd_type=}, {psgn=}, {sval=}' ) 
    9961195    jpi, nperio = lbc_init (ptab, nperio)[1:] 
    9971196    ax = __find_axis__ (ptab, 'x')[0] 
    … …  
    10121211            #> Masks south 
    10131212            # ------------ 
    1014             if nperio in [4, 6] : 
    1015                 ztab [..., 0, : ] = sval 
     1213            if nperio in [4, 6] : ztab [..., 0, : ] = sval 
    10161214 
    10171215            # 
    … …  
    10371235                if cd_type in [ 'T', 'W' ] : # T-, W-point 
    10381236                    ztab [..., -1, jpi//2:  ] = sval 
    1039  
    1040                 if cd_type == 'U' : 
    1041                     ztab [..., -1, jpi//2-1:-1] = sval 
    1042  
    1043                 if cd_type == 'V' : 
    1044                     ztab [..., -1, 1:       ] = sval 
    1045  
    1046                 if cd_type == 'F' : 
    1047                     ztab [..., -1, 0:-1     ] = sval 
     1237                if cd_type == 'U' : ztab [..., -1, jpi//2-1:-1] = sval 
     1238                if cd_type == 'V' : ztab [..., -1, 1:       ] = sval 
     1239                if cd_type == 'F' : ztab [..., -1, 0:-1     ] = sval 
    10481240 
    10491241            if nperio in [5, 6] :            #  North fold F-point pivot 
    … …  
    10621254                    ztab [..., -2, jpi//2+1:-1] = sval 
    10631255 
     1256    PopStack ( 'lbc_plot' ) 
    10641257    return ztab 
    10651258 
    10661259def lbc_add (ptab, nperio=None, cd_type=None, psgn=1) : 
    1067     '''Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2 
     1260    ''' 
     1261    Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2 
    10681262 
    10691263    Periodicity and north fold halos has been removed in NEMO 4.2 
    … …  
    10761270    See NEMO documentation for further details 
    10771271    ''' 
     1272    PushStack ( f'lbc_add ( ptab, {nperio=}, {cd_type=}, {psgn=} )' ) 
    10781273    mmath = __mmath__ (ptab) 
    10791274    nperio = lbc_init (ptab, nperio)[-1] 
    … …  
    10971292                ptab_ext.values[..., :-1, 1:-1] = ptab.values.copy () 
    10981293            else : 
    1099                 if 'X' in lshape     : 
    1100                     ptab_ext.values[...,      1:-1] = ptab.values.copy () 
    1101                 if 'Y' in lshape     : 
    1102                     ptab_ext.values[..., :-1      ] = ptab.values.copy () 
     1294                if 'X' in lshape : ptab_ext.values[...,      1:-1] = ptab.values.copy () 
     1295                if 'Y' in lshape : ptab_ext.values[..., :-1      ] = ptab.values.copy () 
    11031296        else           : 
    11041297            ptab_ext =               np.zeros (ext_shape) 
    … …  
    11061299                ptab_ext       [..., :-1, 1:-1] = ptab.copy () 
    11071300            else : 
    1108                 if 'X' in lshape     : 
    1109                     ptab_ext       [...,      1:-1] = ptab.copy () 
    1110                 if 'Y' in lshape     : 
    1111                     ptab_ext       [..., :-1      ] = ptab.copy () 
     1301                if 'X' in lshape : ptab_ext [...,      1:-1] = ptab.copy () 
     1302                if 'Y' in lshape : ptab_ext [..., :-1      ] = ptab.copy () 
    11121303 
    11131304        if nperio == 4.2 : 
    … …  
    11201311            az = __find_axis__ (ptab, 'z')[0] 
    11211312            at = __find_axis__ (ptab, 't')[0] 
    1122             if az : 
    1123                 ptab_ext = ptab_ext.assign_coords ( {az:ptab.coords[az]} ) 
    1124             if at : 
    1125                 ptab_ext = ptab_ext.assign_coords ( {at:ptab.coords[at]} ) 
     1313            if az : ptab_ext = ptab_ext.assign_coords ( {az:ptab.coords[az]} ) 
     1314            if at : ptab_ext = ptab_ext.assign_coords ( {at:ptab.coords[at]} ) 
    11261315 
    11271316    else : ptab_ext = lbc (ptab, nperio=nperio, cd_type=cd_type, psgn=psgn) 
    11281317 
     1318    PopStack ( 'lbc_add' ) 
    11291319    return ptab_ext 
    11301320 
    … …  
    11411331    See NEMO documentation for further details 
    11421332    ''' 
     1333    PushStack ( f'lbc_del (ptab, {nperio=}, {cd_type=}, {psgn=}' ) 
    11431334    nperio = lbc_init (ptab, nperio)[-1] 
    11441335    #lshape = get_shape (ptab) 
    … …  
    11481339    if nperio in [4.2, 6.2] : 
    11491340        if ax or ay : 
    1150             if ax and ay : 
    1151                 ztab = lbc (ptab[..., :-1, 1:-1], 
    1152                             nperio=nperio, cd_type=cd_type, psgn=psgn) 
     1341            if ax and ay : ztab = lbc (ptab[..., :-1, 1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn) 
    11531342            else : 
    1154                 if ax : 
    1155                     ztab = lbc (ptab[...,      1:-1], 
    1156                                 nperio=nperio, cd_type=cd_type, psgn=psgn) 
    1157                 if ay : 
    1158                     ztab = lbc (ptab[..., -1], 
    1159                                 nperio=nperio, cd_type=cd_type, psgn=psgn) 
    1160         else : 
    1161             ztab = ptab 
    1162     else : 
    1163         ztab = ptab 
    1164  
     1343                if ax    : ztab = lbc (ptab[...,      1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn) 
     1344                if ay    : ztab = lbc (ptab[..., -1]       , nperio=nperio, cd_type=cd_type, psgn=psgn) 
     1345        else : ztab = ptab 
     1346    else : ztab = ptab 
     1347 
     1348    PopStack ( 'lbc_del' ) 
    11651349    return ztab 
    11661350 
    … …  
    11761360    See NEMO documentation for further details 
    11771361    ''' 
    1178  
     1362    PushStack ( f'lbc_index ( {jj=}, {ii=}, {jpj=}, {jpi=}, {nperio=}, {cd_type=} )' ) 
    11791363    if nperio is None : 
    11801364        nperio = __guess_nperio__ (jpj, jpi, nperio) 
    … …  
    11861370 
    11871371    mmath = __mmath__ (jj) 
    1188     if mmath is None : 
    1189         mmath=np 
     1372    if not mmath : mmath=np 
    11901373 
    11911374    # 
    … …  
    11991382    # 
    12001383    def mod_ij (cond, jy_new, ix_new) : 
     1384        PushStack ( f'mod_ij (cond, jy_new, ix_new)' ) 
    12011385        jy_r = mmath.where (cond, jy_new, jy) 
    12021386        ix_r = mmath.where (cond, ix_new, ix) 
     1387        PopStack ( 'mod_ij' ) 
    12031388        return jy_r, ix_r 
    12041389    # 
    … …  
    12511436 
    12521437    ## Restore convention to Python/C : indexes start at 0 
    1253     jy += -1 
    1254     ix += -1 
    1255  
    1256     if isinstance (jj, int) : 
    1257         jy = jy.item () 
    1258     if isinstance (ii, int) : 
    1259         ix = ix.item () 
    1260  
     1438    jy += -1 ; ix += -1 
     1439 
     1440    if isinstance (jj, int) : jy = jy.item () 
     1441    if isinstance (ii, int) : ix = ix.item () 
     1442 
     1443    PopStack ( f'lbc_index = {jy}, {ix}' ) 
    12611444    return jy, ix 
    12621445 
    1263 def find_ji (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False, out=None) : 
     1446def find_ji (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False, drop_duplicates=False, out=None) : 
    12641447    ''' 
    12651448    Description: seeks J,I indices of the grid point which is the closest 
    12661449       of a given point 
    12671450 
    1268     Usage: go FindJI  <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask] 
     1451    Usage: go find_ji  <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask] 
    12691452    <grid latitudes><grid longitudes> are 2D fields on J/I (Y/X) dimensions 
    12701453    mask : if given, seek only non masked grid points (i.e with mask=1) 
    12711454 
    1272     Example : findIJ (40, -20, nav_lat, nav_lon, mask=1.0) 
     1455    Example : find_ji (40., -20., nav_lat, nav_lon, mask=1.0) 
    12731456 
    12741457    Note : all longitudes and latitudes in degrees 
    12751458 
    1276     Note : may work with 1D lon/lat (?) 
    1277     ''' 
     1459    Note : may work with 1D lon_data/lat_data (?) 
     1460    ''' 
     1461    PushStack ( f'find_ji ( lat_data, lon_data, lat_grid, lon_grid, mask, {verbose=}, {drop_duplicates=}, {out=} ) ') 
    12781462    # Get grid dimensions 
    1279     if len (lon_grid.shape) == 2 : 
    1280         jpi = lon_grid.shape[-1] 
    1281     else                         : 
    1282         jpi = len(lon_grid) 
     1463    if len (lon_grid.shape) == 2 : jpi = lon_grid.shape[-1] 
     1464    else                         : jpi = len(lon_grid) 
    12831465 
    12841466    #mmath = __mmath__ (lat_grid) 
     1467 
     1468    if verbose : 
     1469        print ( f'{lat_data=}' ) 
     1470        print ( f'{lon_data=}' ) 
    12851471 
    12861472    # Compute distance from the point to all grid points (in RADian) 
    … …  
    12881474               + np.cos (RAD*lat_data) * np.cos (RAD*lat_grid) * 
    12891475                 np.cos(RAD*(lon_data-lon_grid)) ) 
     1476 
     1477    if verbose : print ( f'{type(arg)=} {arg.shape=}' ) 
     1478     
    12901479    # Send masked points to 'infinite' 
    12911480    distance = np.arccos (arg) + 4.0*RPI*(1.0-mask) 
    12921481 
    1293     # Truncates to alleviate some precision problem with some grids 
     1482    if verbose : print ( f'{type(distance)=} {distance.shape=}' ) 
     1483 
     1484    # Truncates to alleviate precision problem encountered with some grids 
    12941485    prec = int (1E7) 
    12951486    distance = (distance*prec).astype(int) / prec 
    1296  
    1297     # Compute minimum of distance, and index of minimum 
    1298     # 
    1299     #distance_min = distance.min    () 
    1300     jimin        = int (distance.argmin ()) 
     1487    if verbose : print ( f'{type(distance)=} {distance.shape=}' ) 
     1488         
     1489    # Compute index minimum of distance 
     1490    try : 
     1491        nn = len (lat_data) 
     1492    except :  
     1493        nn = 0 
     1494        jimin = distance.argmin ().astype(int) 
     1495    else :  
     1496        jimin   = np.empty (nn, dtype=int ) 
     1497        for ji in np.arange (nn) : 
     1498            jimin[ji] = distance[ji].argmin() 
     1499    finally :  
     1500        if verbose : 
     1501            print ( f'{type(jimin)=} {jimin.shape}' ) 
     1502            print ( f'{jimin=}' ) 
    13011503 
    13021504    # Compute 2D indices (Python/C flavor : starting at 0) 
    … …  
    13041506    imin = jimin - jmin*jpi 
    13051507 
    1306     # Result 
    13071508    if verbose : 
    1308         # Compute distance achieved 
    1309         #mindist = distance [jmin, imin] 
    1310  
    1311         # Compute azimuth 
    1312         dlon = lon_data-lon_grid[jmin,imin] 
    1313         arg  = np.sin (RAD*dlon) / ( 
    1314                   np.cos(RAD*lat_data)*np.tan(RAD*lat_grid[jmin,imin]) 
    1315                 - np.sin(RAD*lat_data)*np.cos(RAD*dlon)) 
    1316         azimuth = DAR*np.arctan (arg) 
    1317         print ( f'I={imin:d} J={jmin:d} - Data:{lat_data:5.1f}°N {lon_data:5.1f}°E' \ 
    1318             +   f'- Grid:{lat_grid[jmin,imin]:4.1f}°N '   \ 
    1319             +   f'{lon_grid[jmin,imin]:4.1f}°E - Dist: {RA*distance[jmin,imin]:6.1f}km' \ 
    1320             +   f' {DAR*distance[jmin,imin]:5.2f}° ' \ 
    1321             +   f'- Azimuth: {RAD*azimuth:3.2f}RAD - {azimuth:5.1f}°' ) 
    1322  
     1509        print ( f'{jmin=}' ) 
     1510        print ( f'{imin=}' ) 
     1511         
     1512    if drop_duplicates : 
     1513        zz = np.vstack ( (np.array (jmin), np.array (imin)) ) 
     1514        zz = np.swapaxes (zz , 0, 1) 
     1515        zz = np.unique ( zz, axis=0) 
     1516        jmin = zz[:,-2] 
     1517        imin = zz[:,-1] 
     1518 
     1519    PopStack ( 'find_ji' ) 
    13231520    if   out=='dict'                     : 
    1324         return {'x':imin, 'y':jmin} 
     1521        zdict = {UDIMS['y']:jmin, UDIMS['x']:imin} 
     1522        return zdict 
    13251523    elif out in ['array', 'numpy', 'np'] : 
    1326         return np.array ( [jmin, imin] ) 
     1524        return np.array (jmin), np.array (imin) 
    13271525    elif out in ['xarray', 'xr']         : 
    1328         return xr.DataArray ( [jmin, imin] ) 
     1526        jmin = xr.DataArray (jmin, dims=('Num',)) 
     1527        imin = xr.DataArray (imin, dims=('Num',)) 
     1528        jmin.attrs.update ( {'long_name':'j-index' } ) 
     1529        imin.attrs.update ( {'long_name':'i-index' } ) 
     1530        jmin.name = 'j_index' 
     1531        imin.name = 'j_index' 
     1532        return jmin, imin 
    13291533    elif out=='list'                     : 
    13301534        return [jmin, imin] 
    … …  
    13371541    '''Returns curl of a vector field 
    13381542    ''' 
     1543    PushStack ( f'curl ( tx, ty, e1f, e2f, {nperio=} )' ) 
    13391544    ax = __find_axis__ (tx, 'x')[0] 
    13401545    ay = __find_axis__ (ty, 'y')[0] 
    … …  
    13611566    zcurl = lbc (zcurl, nperio=nperio, cd_type='F', psgn=1) 
    13621567 
     1568    PopStack ( 'curl' ) 
    13631569    return zcurl 
    13641570 
    … …  
    13661572    '''Returns divergence of a vector field 
    13671573    ''' 
     1574    PushStack ( f'div  (ux, uy, e1t, e2t, {nperio=}' ) 
    13681575    ax = __find_axis__ (ux, 'x')[0] 
    13691576    ay = __find_axis__ (ux, 'y')[0] 
    … …  
    13901597    zdiv = lbc (zdiv, nperio=nperio, cd_type='T', psgn=1) 
    13911598 
     1599    PopStack ( 'zdiv' ) 
    13921600    return zdiv 
    13931601 
    … …  
    13991607        glam, gphi : longitude and latitude of the points 
    14001608    ''' 
    1401  
     1609    PushStack ( 'geo2en (pxx, pyy, pzz, glam, gphi)' ) 
    14021610    gsinlon = np.sin (RAD * glam) 
    14031611    gcoslon = np.cos (RAD * glam) 
    … …  
    14081616    ptn = - pxx * gcoslon * gsinlat  - pyy * gsinlon * gsinlat + pzz * gcoslat 
    14091617 
     1618    PopStack ( geo2en ) 
    14101619    return pte, ptn 
    14111620 
    … …  
    14171626        glam, gphi : longitude and latitude of the points 
    14181627    ''' 
    1419  
     1628    PushStack ( 'en2geo ( pte, ptn, glam, gphi )' ) 
     1629     
    14201630    gsinlon = np.sin (RAD * glam) 
    14211631    gcoslon = np.cos (RAD * glam) 
    … …  
    14271637    pzz =   ptn * gcoslat 
    14281638 
     1639    PopStack ( 'en2geo' ) 
    14291640    return pxx, pyy, pzz 
    14301641 
    … …  
    14341645    if needed 
    14351646    ''' 
     1647    PushStack ( f'clo_lon (lon, {lon0=}, {rad=}, {deg=} )' ) 
     1648    if rad and deg : 
     1649        raise  
    14361650    mmath = __mmath__ (lon, np) 
    1437     if rad : 
    1438         lon_range = 2.*np.pi 
    1439     if deg : 
    1440         lon_range = 360. 
     1651    if rad : lon_range = 2.*np.pi 
     1652    if deg : lon_range = 360. 
    14411653    c_lon = lon 
    1442     c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, 
    1443                              c_lon-lon_range, clo_lon) 
    1444     c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, 
    1445                              c_lon+lon_range, clo_lon) 
    1446     c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, 
    1447                              c_lon-lon_range, clo_lon) 
    1448     c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, 
    1449                              c_lon+lon_range, clo_lon) 
    1450     if c_lon.shape == () : 
    1451         c_lon = c_lon.item () 
     1654    c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, c_lon-lon_range, c_lon) 
     1655    c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, c_lon+lon_range, c_lon) 
     1656    c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, c_lon-lon_range, c_lon) 
     1657    c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, c_lon+lon_range, c_lon) 
     1658    if c_lon.shape == () : c_lon = c_lon.item () 
    14521659    if mmath == xr : 
    1453         if lon.attrs : 
    1454             c_lon.attrs.update ( lon.attrs ) 
     1660        if lon.attrs : c_lon.attrs.update (lon.attrs) 
     1661 
     1662    PopStack ( 'clo_lon' ) 
    14551663    return c_lon 
    14561664 
    … …  
    14601668    Needed to use transforms in Matplotlib 
    14611669    ''' 
     1670    # No stack here, to avoid problem in Matplotlib 
     1671     
    14621672    jpk = gdept_0.shape[0] 
    14631673    kk = xr.DataArray(pk) 
    … …  
    14741684    Needed to use transforms in Matplotlib 
    14751685    ''' 
     1686    # No stack here, to avoid problem in Matplotlib 
     1687     
    14761688    jpk  = gdept_0.shape[0] 
    14771689    if   isinstance (pz, xr.core.dataarray.DataArray ) : 
    1478         zz   = xr.DataArray (pz.values, dims=('zz',)) 
     1690        zz   = xr.DataArray (pz.values , dims=('zz',)) 
    14791691    elif isinstance (pz,  np.ndarray) : 
    14801692        zz   = xr.DataArray (pz.ravel(), dims=('zz',)) 
    … …  
    14971709    Needed to use transforms in Matplotlib 
    14981710    ''' 
     1711    # No stack here, to avoid problem in Matplotlib 
     1712     
    14991713    jpk = gdept_0.shape[0] 
    15001714    kk = xr.DataArray (pk) 
    … …  
    15121726    Needed to use transforms in Matplotlib 
    15131727    ''' 
     1728    # No stack here, to avoid problem in Matplotlib 
     1729     
    15141730    jpk = gdept_0.shape[0] 
    15151731    if isinstance (pz, xr.core.dataarray.DataArray) : 
    … …  
    15401756    Needed to use transforms in Matplotlib 
    15411757    ''' 
     1758    # No stack here, to avoid problem in Matplotlib 
     1759     
    15421760    #print ('start depth2comp') 
    15431761    if isinstance (pz, xr.core.dataarray.DataArray) : 
    … …  
    15491767    gz = np.where ( zz>depth0, (zz-depth0)*fact+depth0, zz) 
    15501768    #print ( f'depth2comp : {gz=}' ) 
    1551     if type (pz) in [int, float] : 
    1552         return gz.item() 
    1553     else : 
    1554         return gz 
     1769    if type (pz) in [int, float] : return gz.item() 
     1770    else                         : return gz 
    15551771 
    15561772def comp2depth (pz, depth0, fact ) : 
    … …  
    15591775    Needed to use transforms in Matplotlib 
    15601776    ''' 
     1777    # No stack here, to avoid problem in Matplotlib 
     1778     
    15611779    if isinstance (pz, xr.core.dataarray.DataArray) : 
    15621780        zz   = pz.values 
    … …  
    15761794    Needed to use transforms in Matplotlib 
    15771795    ''' 
     1796    # No stack here, to avoid problem in Matplotlib 
     1797     
    15781798    jpi = plon_1d.shape[0] 
    1579     ii = xr.DataArray (pi) 
    1580     i =  np.maximum (0, np.minimum (jpi-1, ii    )) 
    1581     i0 = np.floor (i).astype (int) 
    1582     i1 = np.maximum (0, np.minimum (jpi-1,  i0+1)) 
    1583     xx = i - i0 
    1584     gx = (1.0-xx)*plon_1d[i0]+ xx*plon_1d[i1] 
     1799    ii  = xr.DataArray (pi) 
     1800    i   =  np.maximum (0, np.minimum (jpi-1, ii    )) 
     1801    i0  = np.floor (i).astype (int) 
     1802    i1  = np.maximum (0, np.minimum (jpi-1,  i0+1)) 
     1803    xx  = i - i0 
     1804    gx  = (1.0-xx)*plon_1d[i0]+ xx*plon_1d[i1] 
    15851805    return gx.values 
    15861806 
    … …  
    15901810    Needed to use transforms in Matplotlib 
    15911811    ''' 
     1812    # No stack here, to avoid problem in Matplotlib 
     1813     
    15921814    jpi  = plon_1d.shape[0] 
    15931815    if isinstance (px, xr.core.dataarray.DataArray) : 
    … …  
    16141836    Needed to use transforms in Matplotlib 
    16151837    ''' 
     1838    # No stack here, to avoid problem in Matplotlib 
     1839     
    16161840    jpj = plat_1d.shape[0] 
    16171841    jj  = xr.DataArray (pj) 
    … …  
    16281852    Needed to use transforms in Matplotlib 
    16291853    ''' 
     1854    # No stack here, to avoid problem in Matplotlib 
     1855     
    16301856    jpj = plat_1d.shape[0] 
    16311857    if isinstance (py, xr.core.dataarray.DataArray) : 
    … …  
    16501876    respect to east 
    16511877    ''' 
     1878    PushStack ( f'angle_full ( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif, {nperio=} )')  
    16521879    mmath = __mmath__ (glamt) 
    16531880 
    … …  
    17591986        gcosf = gcosf.assign_coords ( glamf.coords ) 
    17601987 
     1988    PopStack ( 'angle_full' ) 
    17611989    return gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 
    17621990 
    … …  
    17651993    respect to east 
    17661994    ''' 
     1995    PushStack ( f'angle (glam, gphi, {nperio=}, {cd_type=} ) ') 
    17671996    mmath = __mmath__ (glam) 
    17681997 
    … …  
    17972026        gcos = gcos.assign_coords ( glam.coords ) 
    17982027 
     2028    PopStack ( 'angle' ) 
    17992029    return gsin, gcos 
    18002030 
    1801 def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type ) : 
     2031def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type='T' ) : 
    18022032    '''Rotates the Repere: Change vector componantes between 
    18032033    geographic grid --> stretched coordinates grid. 
    … …  
    18052035    All components are on the same grid (T, U, V or F) 
    18062036    ''' 
    1807  
     2037    PushStack ( f'rot_en2ij ( u_e, v_n, gsin, gcos, nperio, {cd_type=} )' ) 
    18082038    u_i = + u_e * gcos + v_n * gsin 
    18092039    v_j = - u_e * gsin + v_n * gcos 
    … …  
    18122042    v_j = lbc (v_j, nperio=nperio, cd_type=cd_type, psgn=-1.0) 
    18132043 
     2044    PopStack ( 'rot_en2ij' ) 
    18142045    return u_i, v_j 
    18152046 
    … …  
    18202051    All components are on the same grid (T, U, V or F) 
    18212052    ''' 
     2053    PushStack ( f'rot_ij2en ( u_i, v_j, gsin, gcos, {nperio=}, {cd_type=} )' ) 
    18222054    u_e = + u_i * gcos - v_j * gsin 
    18232055    v_n = + u_i * gsin + v_j * gcos 
    … …  
    18262058    v_n = lbc (v_n, nperio=nperio, cd_type=cd_type, psgn=1.0) 
    18272059 
     2060    PopStack ( 'rot_ij2en' ) 
    18282061    return u_e, v_n 
    18292062 
    … …  
    18372070    Returns east-north components on the T grid point 
    18382071    ''' 
     2072    PushStack ( f'rot_uv2en ( uo, vo, gsint, gcost, {nperio=}, {zdim=} )' ) 
    18392073    ut = u2t (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 
    18402074    vt = v2t (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 
    … …  
    18462080    v_n = lbc (v_n, nperio=nperio, cd_type='T', psgn=1.0) 
    18472081 
     2082    PopStack ( 'rot_uv2en' ) 
    18482083    return u_e, v_n 
    18492084 
    … …  
    18572092    Returns east-north components on the F grid point 
    18582093    ''' 
     2094    PushStack ( f'rot_uv2enf ( uo, vo, gsint, gcost, {nperio=}, {zdim=} )' ) 
    18592095    uf = u2f (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 
    18602096    vf = v2f (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 
    … …  
    18662102    v_n = lbc (v_n, nperio=nperio, cd_type='F', psgn= 1.0) 
    18672103 
     2104    PopStack ( 'rot_uv2enf' ) 
    18682105    return u_e, v_n 
    18692106 
    … …  
    18712108    '''Interpolates an array from U grid to T grid (i-mean) 
    18722109    ''' 
     2110    PushStack ( f'u2t (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    18732111    mmath = __mmath__ (utab) 
    18742112    utab_0 = mmath.where ( np.isnan(utab), 0., utab) 
    … …  
    18972135            if az != zdim : 
    18982136                ttab = ttab.rename( {az:zdim}) 
     2137                 
     2138    PopStack ( 'u2t' ) 
    18992139    return ttab 
    19002140 
    … …  
    19022142    '''Interpolates an array from V grid to T grid (j-mean) 
    19032143    ''' 
     2144    PushStack ( f'v2t (vtab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    19042145    mmath = __mmath__ (vtab) 
    19052146    #lperio, aperio = lbc_diag (nperio) 
    … …  
    19272168            if az != zdim : 
    19282169                ttab = ttab.rename( {az:zdim}) 
     2170                 
     2171    PopStack ( 'v2t' ) 
    19292172    return ttab 
    19302173 
    … …  
    19322175    '''Interpolates an array from F grid to T grid (i- and j- means) 
    19332176    ''' 
     2177    PushStack ( f'f2t (ftab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    19342178    mmath = __mmath__ (ftab) 
    19352179    ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab) 
    … …  
    19372181    ttab = v2t (f2v (ftab_0, nperio=nperio, psgn=psgn, zdim=zdim, action=action), 
    19382182                     nperio=nperio, psgn=psgn, zdim=zdim, action=action) 
    1939     return lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn) 
     2183 
     2184    ttab = lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn) 
     2185     
     2186    PopStack ( 'f2t' ) 
     2187    return ttab 
    19402188 
    19412189def t2u (ttab, nperio=None, psgn=1.0, zdim=None, action='ave') : 
    19422190    '''Interpolates an array from T grid to U grid (i-mean) 
    19432191    ''' 
     2192    PushStack ( f't2u (ttab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
     2193 
    19442194    mmath = __mmath__ (ttab) 
    19452195    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 
    … …  
    19662216            if az != zdim : 
    19672217                utab = utab.rename( {az:zdim}) 
     2218 
     2219    PopStack ( 't2u' ) 
    19682220    return utab 
    19692221 
    … …  
    19712223    '''Interpolates an array from T grid to V grid (j-mean) 
    19722224    ''' 
     2225    PushStack ( f't2v (ttab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    19732226    mmath = __mmath__ (ttab) 
    19742227    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 
    … …  
    19952248            if az != zdim : 
    19962249                vtab = vtab.rename( {az:zdim}) 
     2250  
     2251    PopStack ( 't2v' ) 
    19972252    return vtab 
    19982253 
    … …  
    20002255    '''Interpolates an array from V grid to F grid (i-mean) 
    20012256    ''' 
     2257    PushStack ( f'v2f (vtab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    20022258    mmath = __mmath__ (vtab) 
    20032259    vtab_0 = mmath.where ( np.isnan(vtab), 0., vtab) 
    … …  
    20242280            if az != zdim : 
    20252281                ftab = ftab.rename( {az:zdim}) 
    2026     return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 
     2282 
     2283    ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 
     2284     
     2285    PopStack ( 'v2f' ) 
     2286    return ftab 
    20272287 
    20282288def u2f (utab, nperio=None, psgn=-1.0, zdim=None, action='ave') : 
    20292289    '''Interpolates an array from U grid to F grid i-mean) 
    20302290    ''' 
     2291    PushStack ( f'u2f (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    20312292    mmath = __mmath__ (utab) 
    20322293    utab_0 = mmath.where ( np.isnan(utab), 0., utab) 
    … …  
    20532314            if az != zdim : 
    20542315                ftab = ftab.rename( {az:zdim}) 
     2316     
     2317    PopStack ( 'u2f' ) 
    20552318    return ftab 
    20562319 
    … …  
    20582321    '''Interpolates an array on T grid to F grid (i- and j- means) 
    20592322    ''' 
     2323    PushStack ( f't2f (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    20602324    mmath = __mmath__ (ttab) 
    20612325    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 
    … …  
    20642328                     nperio=nperio, psgn=psgn, zdim=zdim, action=action) 
    20652329 
    2066     return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 
     2330    ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn)  
     2331 
     2332    PopStack ( 'v2f' ) 
     2333    return ftab 
    20672334 
    20682335def f2u (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') : 
    20692336    '''Interpolates an array on F grid to U grid (j-mean) 
    20702337    ''' 
     2338    PushStack ( f'f2u (utab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    20712339    mmath = __mmath__ (ftab) 
    20722340    ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab) 
    … …  
    20912359        if zdim and az and az != zdim : 
    20922360            utab = utab.rename( {az:zdim}) 
     2361 
     2362    PopStack ( 'f2u' ) 
    20932363    return utab 
    20942364 
    … …  
    20962366    '''Interpolates an array from F grid to V grid (i-mean) 
    20972367    ''' 
     2368    PushStack ( f'f2v (ftab, {nperio=}, {psgn=}, {zdim=}, {action=} )' ) 
    20982369    mmath = __mmath__ (ftab) 
    20992370    ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab) 
    … …  
    21192390            if az != zdim : 
    21202391                vtab = vtab.rename( {az:zdim}) 
     2392 
     2393    PopStack ( 'f2v' ) 
    21212394    return vtab 
    21222395 
    … …  
    21262399    sval is the bottom value 
    21272400    ''' 
     2401    PushStack ( f'w2t (wtab, {zcoord=}, {zdim=}, {sval=} )' ) 
    21282402    mmath = __mmath__ (wtab) 
    21292403    wtab_0 = mmath.where ( np.isnan(wtab), 0., wtab) 
    … …  
    21462420        ttab[..., -1, :, :] = sval 
    21472421 
     2422    PopStack ( 'w2t' ) 
    21482423    return ttab 
    21492424 
    … …  
    21542429    if extrap_surf==True, surface value is taken from 1st level value. 
    21552430    ''' 
     2431    PushStack ( f't2w (utab, {zcoord=}, {zdim=}, {sval=}, {extrap_surf=} )' ) 
    21562432    mmath = __mmath__ (ttab) 
    21572433    ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 
    … …  
    21772453        else : 
    21782454            wtab = wtab.assign_coords ( {zdim:np.arange(ttab.shape[kz])+1.} ) 
     2455 
     2456    PopStack ( 't2w' ) 
    21792457    return wtab 
    21802458 
    … …  
    21862464       nperio, cd_type : periodicity characteristics 
    21872465    ''' 
    2188  
     2466    PushStack ( f'fill (ptab, {perio=}, {cd_type=}, {npass=}, {sval=} ) ') 
    21892467    mmath = __mmath__ (ptab) 
    21902468 
    … …  
    22372515        ztab = lbc_del (ztab, nperio=lperio) 
    22382516 
     2517    PopStack ( 'fill' ) 
    22392518    return ztab 
    22402519 
    … …  
    22562535       modified eastward/nothward components to have correct polar projections in cartopy 
    22572536    ''' 
     2537    PushStack (f' correct_uv (u, v, lat)') 
    22582538    uv = np.sqrt (u*u + v*v)           # Original modulus 
    22592539    zu = u 
    … …  
    22622542    uc = zu*uv/zz 
    22632543    vc = zv*uv/zz      # Final corrected values 
     2544 
     2545    PopStack ( 'correct_uv' ) 
    22642546    return uc, vc 
    22652547 
    … …  
    22672549    '''Returns norm of a 2 components vector 
    22682550    ''' 
    2269     return np.sqrt (u*u + v*v) 
     2551    PushStack ( 'norm_uv (u, v)' ) 
     2552    zz = np.sqrt (u*u + v*v) 
     2553    Poptack ( 'norm_uv' ) 
     2554    return zz 
    22702555 
    22712556def normalize_uv (u, v) : 
    22722557    '''Normalizes 2 components vector 
    22732558    ''' 
     2559    PushStack ( 'normalize_uv (u, v)' ) 
    22742560    uv = norm_uv (u, v) 
    2275     return u/uv, v/uv 
     2561    uu = u/uv 
     2562    vv = v/uv 
     2563    PopStack ( 'normalize_uv' ) 
     2564    return uu, vv 
    22762565 
    22772566def msf (vv, e1v_e3v, plat1d, depthw) : 
    … …  
    22812570    e1v_e3v : prodcut of scale factors e1v*e3v 
    22822571    ''' 
    2283  
     2572    PushStack ( 'msf (vv, e1v_e3v, plat1d, depthw)' ) 
    22842573    v_e1v_e3v = vv * e1v_e3v 
    22852574    v_e1v_e3v.attrs = vv.attrs 
    … …  
    23012590    if az != new_az : 
    23022591        zomsf = zomsf.rename ( {az:new_az} ) 
    2303     zomsf.attrs['standard_name'] = 'Meridional stream function' 
    2304     zomsf.attrs['long_name'] = 'Meridional stream function' 
    2305     zomsf.attrs['units'] = 'Sv' 
    2306     zomsf[new_az].attrs  = depthw.attrs 
     2592    zomsf.attrs ['standard_name'] = 'Meridional stream function' 
     2593    zomsf.attrs ['long_name'] = 'Meridional stream function' 
     2594    zomsf.attrs ['units'] = 'Sv' 
     2595    zomsf [new_az].attrs  = depthw.attrs 
    23072596    zomsf.lat.attrs=plat1d.attrs 
    23082597 
     2598    PopStack ( 'msf' ) 
    23092599    return zomsf 
    23102600 
    … …  
    23182608         bsf0={'x':5, 'y':300} for eORCA1 
    23192609    ''' 
     2610    PushStack ( f'bsf (uu, e2u_e3u, mask, {nperio=}, {bsf0=} )' ) 
     2611     
    23202612    u_e2u_e3u       = uu * e2u_e3u 
    23212613    u_e2u_e3u.attrs = uu.attrs 
    … …  
    23372629    zbsf = lbc (zbsf, nperio=nperio, cd_type='F') 
    23382630 
     2631    PopStack ( 'bsf' ) 
    23392632    return zbsf 
    23402633 
    … …  
    23532646 
    23542647        ''' 
     2648        PushStack ( f'namelist_read (ref, cfg, {out=}, {flat=}, {verbose=} )' ) 
    23552649        if ref : 
    23562650            if isinstance (ref, str) : 
    … …  
    24072701                            print (nam, ':', value, ':', nml[nam][value]) 
    24082702 
     2703        PopStack ( 'namelist_read' ) 
    24092704        if out == 'dict' : 
    24102705            return dict_namelist 
    … …  
    24182713        Read NEMO namelist(s) and return either a dictionnary or an xarray dataset 
    24192714        ''' 
     2715        PushStack ( f'namelist_read [void version] (ref, cfg, {out=}, {flat=}, {verbose=} )' ) 
     2716 
    24202717        print ( 'Error : module f90nml not found' ) 
    24212718        print ( 'Cannot call namelist_read' ) 
    24222719        print ( 'Call parameters where : ') 
    24232720        print ( f'{err=} {ref=} {cfg=} {out=} {flat=} {verbose=}' ) 
     2721        PopStack ( 'namelist_read [void version]' ) 
    24242722 
    24252723def fill_closed_seas (imask, nperio=None,  cd_type='T') : 
    … …  
    24282726    imask : mask, 1 on ocean, 0 on land 
    24292727    ''' 
     2728    PushStack ( f'fill_closed_seas (imask, {nperio=}, {cd_type=} )' ) 
    24302729    from scipy import ndimage 
    24312730 
    … …  
    24332732    imask_filled = lbc ( imask_filled, nperio=nperio, cd_type=cd_type) 
    24342733 
     2734    PopStack ( 'fill_closed_seas' ) 
    24352735    return imask_filled 
    24362736 
    … …  
    24382738# Sea water state function parameters from NEMO code 
    24392739 
    2440 RDELTAS = 32. 
     2740RDELTAS = 32.0 
    24412741R1_S0   = 0.875/35.16504 
    2442 R1_T0   = 1./40. 
    2443 R1_Z0   = 1.e-4 
     2742R1_T0   = 1.0/40. 
     2743R1_Z0   = 1.0e-4 
    24442744 
    24452745EOS000 =  8.0189615746e+02 
    … …  
    24962796EOS013 =  3.7969820455e-01 
    24972797 
    2498 def rhop ( ptemp, psal ) : 
     2798def rhop (ptemp, psal) : 
    24992799    '''Returns potential density referenced to surface 
    25002800 
    25012801    Computation from NEMO code 
    25022802    ''' 
     2803    PushStack ( 'rhop ( ptemp, psal )' ) 
    25032804    zt      = ptemp * R1_T0                                  # Temperature (°C) 
    25042805    zs      = np.sqrt ( np.abs( psal + RDELTAS ) * R1_S0 )   # Square root of salinity (PSS) 
    … …  
    25132814         + (((((EOS600*zs+ EOS500)*zs + EOS400)*zs + EOS300)*zs + EOS200)*zs + EOS100)*zs + EOS000 ) 
    25142815    # 
     2816    PopStack ( 'rhop' ) 
    25152817    return prhop 
    25162818 
    … …  
    25202822    Computation from NEMO code 
    25212823    ''' 
     2824    PushStack ( 'rho ( ptemp, psal) ' ) 
    25222825    zh      = pdep  * R1_Z0                                  # Depth (m) 
    25232826    zt      = ptemp * R1_T0                                  # Temperature (°C) 
    … …  
    25432846         + ((((EOS510*zs  + EOS410)*zs + EOS310)*zs + EOS210)*zs + EOS110)*zs + EOS010)*zt 
    25442847         + (((((EOS600*zs + EOS500)*zs + EOS400)*zs + EOS300)*zs + 
    2545                                        EOS200)*zs + EOS100)*zs + EOS000 ) 
     2848                                         EOS200)*zs + EOS100)*zs + EOS000 ) 
    25462849    # 
    25472850    prho  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    25482851    # 
     2852    PopStack ( 'rho' ) 
    25492853    return prho 
    25502854 
Note: See TracChangeset for help on using the changeset viewer.