Index: /TOOLS/MOSAIX/Build_coordinates_mask.py
===================================================================
--- /TOOLS/MOSAIX/Build_coordinates_mask.py	(revision 5946)
+++ /TOOLS/MOSAIX/Build_coordinates_mask.py	(revision 5947)
@@ -10,141 +10,78 @@
 import xarray as xr
 import nemo
-import datetime, os, platform
+import datetime, os, platform, argparse
+
+## SVN information
+__Author__   = "$Author: $"
+__Date__     = "$Date: $"
+__Revision__ = "$Revision: $"
+__Id__       = "$Id: $"
+__HeadURL__  = "$HeadURL: $"
+__SVN_Date__ = "$SVN_Date: $"
+###
+
+### ===== Handling command line parameters ==================================================
+# Creating a parser
+parser = argparse.ArgumentParser (
+    description = """Read coordinates and bathymetry to build masks, grid bounds and areas for MOSAIX""",
+    epilog='-------- End of the help message --------')
+
+# Adding arguments
+parser.add_argument ('--model'   , help='oce model name', type=str, default='eORCA1.2', choices=('paleORCA', 'ORCA2.3', 'ORCA2.4', 'eORCA1.2', 'eORCA025', 'eORCA025.1') )
+parser.add_argument ('--bathy'   , help='bathymetry file name', type=str, default='bathy_meter.nc' )
+parser.add_argument ('--coord'   , help='coordinates file name', type=str, default='coordinates.nc' )
+parser.add_argument ('--fout'    , help='output file name (given as an input to MOSAIX)', type=str, default='coordinates_mask.nc' )
+parser.add_argument ('--ocePerio', help='periodicity of ocean grid', type=int, default=0, choices=(1, 2, 3, 4, 5, 6) )
+parser.add_argument ('--nemo4Ubug',help='reproduce NEMO lbc_lnk bug for U grid and periodicity 4', type=bool, default=False, choices=(True, False) )
+parser.add_argument ('--straits' , help='modify grid metrics for selected strait points (depends on model name)', type=bool, default=False, choices=(True, False) )
+parser.add_argument ('--coordPerio', help='impose periodicity of coordinates and areas through lbc', type=bool, default=False, choices=(True, False) )
+parser.add_argument ('--maskbathy', help='use the same formula as in NEMO to compute mask_T from the bathymetry file. Does not always give the same result as NEMO ???', type=bool, default=False, choices=(True, False) )
+
+# Parse command line
+myargs = parser.parse_args()
+
+# ocean model name
+model = myargs.model
+
+# bathymetry input file
+n_bathy = myargs.bathy
+# coordinates input file
+n_coord = myargs.coord
+# coordinates and mask output file
+n_out = myargs.fout
 
 # reproduce periodicity lbclnk bug for U and nperio = 4
-nemoUbug = False
+nemoUbug = myargs.nemo4Ubug
 # change metrics (and bathymetry) for straits
-straits = False
+straits = myargs.straits
 # periodicity imposed on coordinates, metrics and areas
-coordperio = False
+coordperio = myargs.coordPerio
 # function used to compute mask_T from bathymetry
-maskbathynemo = False
-
-# type of grid treated : ORCA 1, 2, paleorca # to be improved
-
-#model   = 'eORCA1.2'
-#model   = 'ORCA2.3'
-#model   = 'ORCA2.4'
-model   = 'paleORCA'
-#model   = 'ORCA.05'
-
-if model == 'eORCA1.2' :
-    n_coord = 'eORCA1.2coordinates.nc' #'eORCA_R1_coordinates_v1.0.nc' 
-    n_bathy = 'eORCA1.2bathy_meter.nc' #'eORCA_R1_bathy_meter_v2.2.nc'
-    #n_coord = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_coordinates_v1.0.nc'
-    #n_bathy = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_bathy_meter_v2.2.nc'
-    nperio  = 6
-    n_out   = 'eORCA1.2_coordinates_mask_test.nc'
-
-if model == 'ORCA2.3' :
-    nperio  = 4
-    n_coord = 'coordinates_ORCA2.nc' # 'ORCA2.3_coordinates.nc' # 'coordinates_ORCA2_LIM3_PISCES.nc' #
-    n_bathy = 'bathy_meter_ORCA2.nc' # 'ORCA2.3_bathy_meter.nc' # 'bathy_meter_ORCA2_LIM3_PISCES.nc' #
-    n_out   = 'My_coordinates_xios_ORCA2.nc'
-    
-if model == 'ORCA2.4' :
-    nperio  = 4
-    n_coord = 'ORCA2.4_coordinates.nc'
-    n_bathy = 'ORCA2.4_bathy_meter.nc'
-    n_out   = 'ORCA2.4_coordinates_mask_test.nc'
-
-if model == 'paleORCA' :
-    n_coord = 'coordinates_paleo.nc'
-    n_bathy = 'bathy_meter_paleo.nc'
-    nperio = 6
-    n_out = 'paleo_coordinates_mask_test.nc'
-
-#---------------------------------------------------------------
-# function to generate grid bounds from NEMO coordinates.nc file
-def set_bounds (cdgrd) :
-    '''
-    Constructs lon/lat bounds
-    Bounds are numerated counter clockwise, from bottom left
-    See NEMO file OPA_SRC/IOM/iom.F90, ROUTINE set_grid_bounds, for more details
-    '''
-    # Define offset of coordinate representing bottom-left corner
-    if cdgrd in ['T', 'W'] : 
-        icnr = -1 ; jcnr = -1
-        corner_lon = f_coord ['glamf'].copy() ; corner_lat = f_coord ['gphif'].copy()
-        center_lon = f_coord ['glamt'].copy() ; center_lat = f_coord ['gphit'].copy()
-    if cdgrd == 'U'        : 
-        icnr =  0 ; jcnr = -1
-        corner_lon = f_coord ['glamv'].copy() ; corner_lat = f_coord ['gphiv'].copy()
-        center_lon = f_coord ['glamu'].copy() ; center_lat = f_coord ['gphiu'].copy()
-    if cdgrd == 'V'        : 
-        icnr = -1 ; jcnr =  0
-        corner_lon = f_coord ['glamu'].copy() ; corner_lat = f_coord ['gphiu'].copy()
-        center_lon = f_coord ['glamv'].copy() ; center_lat = f_coord ['gphiv'].copy()
-    if cdgrd == 'F'        : 
-        icnr = -1 ; jcnr = -1
-        corner_lon = f_coord ['glamt'].copy() ; corner_lat = f_coord ['gphit'].copy()
-        center_lon = f_coord ['glamf'].copy() ; center_lat = f_coord ['gphif'].copy()
-      
-    #y_grid, x_grid = corner_lon.shape ;
-    nvertex = 4
-    dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]
-    
-    bounds_lon = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
-    bounds_lat = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
-    
-    idx = [(jcnr,icnr), (jcnr,icnr+1), (jcnr+1,icnr+1), (jcnr+1,icnr)]
-    
-    # Compute cell vertices that can be defined, 
-    # and complete with periodicity
-    for nn in range (nvertex) :
-        tmp = np.roll (corner_lon, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
-        bounds_lon[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
-        tmp = np.roll (corner_lat, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
-        bounds_lat[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
-        bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
-        bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
-    
-    # Zero-size cells at closed boundaries if cell points provided,
-    # otherwise they are closed cells with unrealistic bounds
-    if not (nperio == 1 or nperio == 4 or nperio == 6) :
-        for nn in range (nvertex) : 
-            bounds_lon[:,0,nn]   = center_lon[:,0]   # (West or jpni = 1), closed E-W
-            bounds_lat[:,0,nn]   = center_lat[:,0]
-            bounds_lon[:,jpi,nn] = center_lon[:,jpi] # (East or jpni = 1), closed E-W
-            bounds_lat[:,jpi,nn] = center_lat[:,jpi]
-    if nperio != 2 :
-        for nn in range (nvertex) :
-            bounds_lon[0,:,nn] = center_lon[0,:] # (South or jpnj = 1), not symmetric
-            bounds_lat[0,:,nn] = center_lat[0,:]
-    if nperio < 3 :
-        for nn in range (nvertex) :
-            bounds_lon[jpj,:,nn] = center_lon[jpj,:] # (North or jpnj = 1), no north fold
-            bounds_lat[jpj,:,nn] = center_lat[jpj,:]
-    
-    # Rotate cells at the north fold
-    if nperio >= 3 :
-        # Working array for location of northfold
-        z_fld = nemo.lbc (np.ones ((jpj, jpi)), nperio=nperio, cd_type=cdgrd, psgn=-1., nemo_4U_bug=nemoUbug)
-        z_fld = np.repeat((z_fld == -1.0)[...,np.newaxis],4,axis=2)
-        # circular shift of 2 indices in bounds third index
-        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
-        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
-        bounds_lon[:,:,:]  = np.where (z_fld, bounds_lon_tmp[:,:,:] , bounds_lon[:,:,:] )
-        bounds_lat[:,:,:]  = np.where (z_fld, bounds_lat_tmp[:,:,:] , bounds_lat[:,:,:] )
-    
-    # Invert cells at the symmetric equator
-    if nperio == 2 :
-        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
-        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
-        bounds_lon [0,:,:] = bounds_lon [0,:,:]
-        bounds_lat [0,:,:] = bounds_lat [0,:,:]
-    
-    #bounds_lon.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
-    #bounds_lat.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
-    bounds_lon.attrs['units']       = 'degrees_east'
-    bounds_lat.attrs['units']       = 'degrees_north'
-    bounds_lon.name = 'bounds_lon_grid_' + cdgrd
-    bounds_lat.name = 'bounds_lat_grid_' + cdgrd
-    # remove _FillValue
-    bounds_lon.encoding['_FillValue'] = None
-    bounds_lat.encoding['_FillValue'] = None
-
-    return bounds_lon, bounds_lat
-#------------------------------------------------------------
+maskbathynemo = myargs.maskbathy
+
+# type of grid periodicity
+nperio = myargs.ocePerio
+
+# check of periodicity with type of grid 
+if model in ('eORCA1.2', 'paleORCA', 'ORCA025', 'eORCA025', 'eORCA025.1') :
+    if nperio != 0 :
+        if nperio != 6 :
+            print(f'Warning ! model = {model} and ocePerio = {nperio} instead of 6 !')
+    else :
+        nperio = 6
+
+if model in ('ORCA2.3', 'ORCA2.4') :
+    if nperio != 0 :
+        if nperio != 4 :
+            print(f'Attention ! model = {model} and ocePerio = {nperio} instead of 4 !')
+    else :
+        nperio = 4
+
+print(f' model = {model}\n ocePerio = {nperio}\n bathy = {n_bathy}\n coord = {n_coord}\n fout  = {n_out}')
+print(f' Switchs : nemo4Ubug = {nemoUbug}, straits = {straits}, coordPerio = {coordperio}, maskbathy = {maskbathynemo}')
+
+##
+#!!! bathymetry and coordinates files
+##
 
 # open input files while removing the time dimension
@@ -160,4 +97,5 @@
     print ('failed to suppress time')
 
+# rename latitude, longitude and grid variables in bathymetry
 Bathymetry = f_bathy['Bathymetry'].copy()
 
@@ -174,4 +112,5 @@
 Bathymetry = Bathymetry.rename ({'y':'y_grid_T', 'x':'x_grid_T'})
 
+# modify Bathymetry in straits for select grids (might change computed mask_T)
 if straits :
     # Open straits for usual grids
@@ -182,32 +121,20 @@
         Bathymetry[87,159] = 137.
 
-# Determine periodicity from grid type
-jpj, jpi = Bathymetry.shape
-try : 
-    if nperio != None :
-        print ('nperio specified : ', nperio)
-except :
-    nperio = None
-    if jpi ==  182 : nperio = 4 # ORCA2. \!/ We choose legacy orca2
-    if jpi ==  362 : nperio = 6 # ORCA1.
-    if jpi == 1442 : nperio = 6 # ORCA025.
-    #
-    if nperio == None :
-        sys.exit ('nperio not found, and cannot by guessed' )
-    else :
-        print ('nperio set as {:d} (deduced from data size {:d} {:d})'.format (nperio, jpj, jpi))
-
+# impose periodicity of bathymetry
 Bathymetry = nemo.lbc (Bathymetry, nperio=nperio, cd_type='T')
-#Suppress ocean points at southernmost position only if nperio in {3, 4, 5, 6}
+# suppress ocean points at southernmost position only if nperio in {3, 4, 5, 6}
 if nperio in (3, 4, 5, 6) :
-    Bathymetry[0,:] = 0.0 
-
-# Creation of mask_T
-
+    Bathymetry[0,:] = 0.0
+
+##
+#!!! Create masks from bathymetry
+##
+
+# Creation of mask_T following choosed option maskbathy
 if maskbathynemo :
     # Use same formula as domzgr. 
-    mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='int8')
+    mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='f4')
 else :
-    mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='int8')
+    mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='f4')
 
 # Creation of U, V, W, F masks from mask_T
@@ -217,143 +144,69 @@
 mask_W = mask_T
 
-mask_U = nemo.lbc (mask_U, nperio=nperio, cd_type='U', nemo_4U_bug=nemoUbug).astype (dtype='int8')
-mask_V = nemo.lbc (mask_V, nperio=nperio, cd_type='V').astype (dtype='int8')
-mask_F = nemo.lbc (mask_F, nperio=nperio, cd_type='F').astype (dtype='int8')
-mask_W = nemo.lbc (mask_W, nperio=nperio, cd_type='W').astype (dtype='int8')
-
-mask_U = mask_U.rename ( {'y_grid_T'      : 'y_grid_U'      , 'x_grid_T'      : 'x_grid_U', 
-                          'nav_lat_grid_T': 'nav_lat_grid_U', 'nav_lon_grid_T': 'nav_lon_grid_U'} )
-mask_V = mask_V.rename ( {'y_grid_T'      : 'y_grid_V'      , 'x_grid_T'      : 'x_grid_V',
-                          'nav_lat_grid_T': 'nav_lat_grid_V', 'nav_lon_grid_T': 'nav_lon_grid_V'} )
-mask_W = mask_W.rename ( {'y_grid_T'      : 'y_grid_W'      , 'x_grid_T'      : 'x_grid_W',
-                          'nav_lat_grid_T': 'nav_lat_grid_W', 'nav_lon_grid_T': 'nav_lon_grid_W'} )
-mask_F = mask_F.rename ( {'y_grid_T'      : 'y_grid_F'      , 'x_grid_T'      : 'x_grid_F',
-                          'nav_lat_grid_T': 'nav_lat_grid_F', 'nav_lon_grid_T': 'nav_lon_grid_F'} )
-
-mask_T.name = 'mask_T'
-mask_U.name = 'mask_U'
-mask_V.name = 'mask_V'
-mask_F.name = 'mask_F'
-mask_W.name = 'mask_W'
-
-# create masks without duplicate points : maskutil_[TUVWF]
+# loop on TUVFW to modify coordinates and attributes of mask_[TUVFW] and maskutil_[TUVFW]
 for cd_type in ['T', 'U', 'V', 'F', 'W'] :
     MaskName = 'mask_'     + cd_type
     UtilName = 'maskutil_' + cd_type
+    # impose periodicity of chosen grid model on masks
+    locals()[MaskName] = nemo.lbc (locals()[MaskName], nperio=nperio, cd_type=cd_type, nemo_4U_bug=nemoUbug).astype (dtype='f4')
+    # rename masks coordinates
+    if cd_type != 'T' :
+        locals()[MaskName] = locals()[MaskName].rename \
+                    ( {'y_grid_T'      :       'y_grid_'+cd_type, 'x_grid_T'      :       'x_grid_'+cd_type,
+                       'nav_lat_grid_T': 'nav_lat_grid_'+cd_type, 'nav_lon_grid_T': 'nav_lon_grid_'+cd_type} )
+
+    # create masks without duplicate points : maskutil_[TUVWF]
     locals()[UtilName] = nemo.lbc_mask (locals()[MaskName].copy(), nperio=nperio, cd_type=cd_type)
+ 
+    #set name attribute of mask dataset
+    locals()[MaskName].name = MaskName
     locals()[UtilName].name = UtilName
 
-# add masks attributes
-mask_T.attrs    ['cell_measures'] = 'area: area_grid_T'
-mask_U.attrs    ['cell_measures'] = 'area: area_grid_U'
-mask_V.attrs    ['cell_measures'] = 'area: area_grid_V'
-mask_W.attrs    ['cell_measures'] = 'area: area_grid_W'
-mask_F.attrs    ['cell_measures'] = 'area: area_grid_F'
-
-maskutil_T.attrs['cell_measures'] = 'area: area_grid_T'
-maskutil_U.attrs['cell_measures'] = 'area: area_grid_U'
-maskutil_V.attrs['cell_measures'] = 'area: area_grid_V'
-maskutil_W.attrs['cell_measures'] = 'area: area_grid_W'
-maskutil_F.attrs['cell_measures'] = 'area: area_grid_F'
-
-# create grid coordinates from NEMO coordinates.nc file
-if coordperio :
-    nav_lon_grid_T = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='T')
-    nav_lat_grid_T = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='T')
-    nav_lon_grid_U = nemo.lbc (f_coord ['glamu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
-    nav_lat_grid_U = nemo.lbc (f_coord ['gphiu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
-    nav_lon_grid_V = nemo.lbc (f_coord ['glamv'].copy(), nperio=nperio, cd_type='V')
-    nav_lat_grid_V = nemo.lbc (f_coord ['gphiv'].copy(), nperio=nperio, cd_type='V')
-    nav_lon_grid_W = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='W')
-    nav_lat_grid_W = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='W')
-    nav_lon_grid_F = nemo.lbc (f_coord ['glamf'].copy(), nperio=nperio, cd_type='F')
-    nav_lat_grid_F = nemo.lbc (f_coord ['gphif'].copy(), nperio=nperio, cd_type='F')
-else :
-    nav_lon_grid_T = f_coord ['glamt'].copy()
-    nav_lat_grid_T = f_coord ['gphit'].copy()
-    nav_lon_grid_U = f_coord ['glamu'].copy()
-    nav_lat_grid_U = f_coord ['gphiu'].copy()
-    nav_lon_grid_V = f_coord ['glamv'].copy()
-    nav_lat_grid_V = f_coord ['gphiv'].copy()
-    nav_lon_grid_W = f_coord ['glamt'].copy()
-    nav_lat_grid_W = f_coord ['gphit'].copy()
-    nav_lon_grid_F = f_coord ['glamf'].copy()
-    nav_lat_grid_F = f_coord ['gphif'].copy()
-
-nav_lon_grid_T = nav_lon_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )
-nav_lat_grid_T = nav_lat_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )
-nav_lon_grid_U = nav_lon_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )
-nav_lat_grid_U = nav_lat_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )
-nav_lon_grid_V = nav_lon_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )
-nav_lat_grid_V = nav_lat_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )
-nav_lon_grid_W = nav_lon_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )
-nav_lat_grid_W = nav_lat_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )
-nav_lon_grid_F = nav_lon_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )
-nav_lat_grid_F = nav_lat_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )
+    # remove _FillVallue from mskutil
+    locals()[MaskName].encoding['_FillValue'] = None
+    locals()[UtilName].encoding['_FillValue'] = None
+
+    # add masks attributes
+    locals()[MaskName].attrs['cell_measures'] = 'area: area_grid_'+cd_type
+    locals()[UtilName].attrs['cell_measures'] = 'area: area_grid_'+cd_type
+
+##
+#!!! create grid coordinates from NEMO coordinates.nc file
+##
+
+angle = { 'lon' : 'glam', 'lat' : 'gphi' }
+gridv = { 'T' : 't', 'U' : 'u', 'V' : 'v', 'F' : 'f', 'W' : 't' }
+for cd_type in ['T', 'U', 'V', 'F', 'W'] :
+    for dir_type in ['lon', 'lat'] :
+        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
+        dir_name = angle[dir_type]+gridv[cd_type]
+# impose or not periodicity on coordinates read from file
+        if coordperio :
+            locals()[coord_name] = nemo.lbc (f_coord [dir_name].copy(), nperio=nperio, cd_type=cd_type, nemo_4U_bug=nemoUbug)
+        else :
+            locals()[coord_name] = f_coord [dir_name].copy()
+
+        locals()[coord_name] = locals()[coord_name].rename( {'y':'y_grid_'+cd_type, 'x':'x_grid_'+cd_type} )
 
 # remove _FillValue and missing_value
-for cd_type in ['T', 'U', 'V', 'F', 'W'] :
-    for dir_type in ['lat', 'lon'] :
-        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
         locals()[coord_name].encoding['_FillValue'] = None
         locals()[coord_name].encoding['missing_value'] = None
 
-nav_lon_grid_T.name = 'nav_lon_grid_T'
-nav_lat_grid_T.name = 'nav_lat_grid_T'
-nav_lon_grid_U.name = 'nav_lon_grid_U'
-nav_lat_grid_U.name = 'nav_lat_grid_U'
-nav_lon_grid_V.name = 'nav_lon_grid_V'
-nav_lat_grid_V.name = 'nav_lat_grid_V'
-nav_lon_grid_F.name = 'nav_lon_grid_F'
-nav_lat_grid_F.name = 'nav_lat_grid_F'
-nav_lon_grid_W.name = 'nav_lon_grid_W'
-nav_lat_grid_W.name = 'nav_lat_grid_W'
+# define name attribute
+        locals()[coord_name].name = coord_name
 
 # add coordinates attributes
-nav_lon_grid_T.attrs['standard_name'] = 'longitude'
-nav_lon_grid_T.attrs['long_name']     = 'Longitude'
-nav_lon_grid_T.attrs['units']         = 'degrees_east'
-nav_lat_grid_T.attrs['standard_name'] = 'latitude'
-nav_lat_grid_T.attrs['long_name']     = 'Latitude'
-nav_lat_grid_T.attrs['units']         = 'degrees_north'
-
-nav_lon_grid_U.attrs['standard_name'] = 'longitude'
-nav_lon_grid_U.attrs['long_name']     = 'Longitude'
-nav_lon_grid_U.attrs['units']         = 'degrees_east'
-nav_lat_grid_U.attrs['standard_name'] = 'latitude'
-nav_lat_grid_U.attrs['long_name']     = 'Latitude'
-nav_lat_grid_U.attrs['units']         = 'degrees_north'
-
-nav_lon_grid_V.attrs['standard_name'] = 'longitude'
-nav_lon_grid_V.attrs['long_name']     = 'Longitude'
-nav_lon_grid_V.attrs['units']         = 'degrees_east'
-nav_lat_grid_V.attrs['standard_name'] = 'latitude'
-nav_lat_grid_V.attrs['long_name']     = 'Latitude'
-nav_lat_grid_V.attrs['units']         = 'degrees_north'
-
-nav_lon_grid_W.attrs['standard_name'] = 'longitude'
-nav_lon_grid_W.attrs['long_name']     = 'Longitude'
-nav_lon_grid_W.attrs['units']         = 'degrees_east'
-nav_lat_grid_W.attrs['standard_name'] = 'latitude'
-nav_lat_grid_W.attrs['long_name']     = 'Latitude'
-nav_lat_grid_W.attrs['units']         = 'degrees_north'
-
-nav_lon_grid_F.attrs['standard_name'] = 'longitude'
-nav_lon_grid_F.attrs['long_name']     = 'Longitude'
-nav_lon_grid_F.attrs['units']         = 'degrees_east'
-nav_lat_grid_F.attrs['standard_name'] = 'latitude'
-nav_lat_grid_F.attrs['long_name']     = 'Latitude'
-nav_lat_grid_F.attrs['units']         = 'degrees_north'
-
-nav_lon_grid_T.attrs['bounds'] = 'bounds_lon_grid_T'
-nav_lat_grid_T.attrs['bounds'] = 'bounds_lat_grid_T'
-nav_lon_grid_U.attrs['bounds'] = 'bounds_lon_grid_U'
-nav_lat_grid_U.attrs['bounds'] = 'bounds_lat_grid_U'
-nav_lon_grid_V.attrs['bounds'] = 'bounds_lon_grid_V'
-nav_lat_grid_V.attrs['bounds'] = 'bounds_lat_grid_V'
-nav_lon_grid_W.attrs['bounds'] = 'bounds_lon_grid_W'
-nav_lat_grid_W.attrs['bounds'] = 'bounds_lat_grid_W'
-nav_lon_grid_F.attrs['bounds'] = 'bounds_lon_grid_F'
-nav_lat_grid_F.attrs['bounds'] = 'bounds_lat_grid_F'
+        locals()[coord_name].attrs['bounds'] ='bounds_' + dir_type + '_grid_' + cd_type
+
+    locals()['nav_lon_grid_'+cd_type].attrs['standard_name'] = 'longitude'
+    locals()['nav_lon_grid_'+cd_type].attrs['long_name']     = 'Longitude'
+    locals()['nav_lon_grid_'+cd_type].attrs['units']         = 'degrees_east'
+    locals()['nav_lat_grid_'+cd_type].attrs['standard_name'] = 'latitude'
+    locals()['nav_lat_grid_'+cd_type].attrs['long_name']     = 'Latitude'
+    locals()['nav_lat_grid_'+cd_type].attrs['units']         = 'degrees_north'
+
+##
+#!!! compute areas of cells at coordinate points
+##
 
 # create areas variables for grid cells from NEMO coordinates.nc file
@@ -367,6 +220,4 @@
         # some point are undefined but you need to have e1 and e2 .NE. 0
         locals()[coordName]=xr.where(locals()[coordName] == 0.0, 1.0e2, locals()[coordName])
-
-#Compute cells areas
 
 # Correct areas for straits
@@ -417,45 +268,122 @@
         e1v[275,661] =  25.e3
 
-if coordperio :
-    area_grid_T = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='T')
-    area_grid_U = nemo.lbc (e1u*e2u, nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
-    area_grid_V = nemo.lbc (e1v*e2v, nperio=nperio, cd_type='V')
-    area_grid_W = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='W')
-    area_grid_F = nemo.lbc (e1f*e2f, nperio=nperio, cd_type='F')
-else :
-    area_grid_T = e1t*e2t
-    area_grid_U = e1u*e2u
-    area_grid_V = e1v*e2v
-    area_grid_W = e1t*e2t
-    area_grid_F = e1f*e2f
-
-area_grid_T.name = 'area_grid_T'
-area_grid_U.name = 'area_grid_U'
-area_grid_V.name = 'area_grid_V'
-area_grid_F.name = 'area_grid_F'
-area_grid_W.name = 'area_grid_W'
-
-area_grid_T = area_grid_T.rename ({'y':'y_grid_T', 'x':'x_grid_T'})
-area_grid_U = area_grid_U.rename ({'y':'y_grid_U', 'x':'x_grid_U'})
-area_grid_V = area_grid_V.rename ({'y':'y_grid_V', 'x':'x_grid_V'})
-area_grid_W = area_grid_W.rename ({'y':'y_grid_W', 'x':'x_grid_W'})
-area_grid_F = area_grid_F.rename ({'y':'y_grid_F', 'x':'x_grid_F'})
-
-area_grid_T.attrs['standard_name'] = 'cell_area'
-area_grid_T.attrs['units']         = 'm2'
-area_grid_U.attrs['standard_name'] = 'cell_area'
-area_grid_U.attrs['units']         = 'm2'
-area_grid_V.attrs['standard_name'] = 'cell_area'
-area_grid_V.attrs['units']         = 'm2'
-area_grid_W.attrs['standard_name'] = 'cell_area'
-area_grid_W.attrs['units']         = 'm2'
-area_grid_F.attrs['standard_name'] = 'cell_area'
-area_grid_F.attrs['units']         = 'm2'
+# compute cells areas
 
 for cd_type in ['T', 'U', 'V', 'F', 'W'] :
     areaName = 'area_grid_' + cd_type
+    if coordperio :
+        locals()[areaName] = nemo.lbc (locals()['e1'+gridv[cd_type]]*locals()['e2'+gridv[cd_type]],
+                                                           nperio=nperio, cd_type=cd_type, nemo_4U_bug=nemoUbug)
+    else :
+        locals()[areaName] = locals()['e1'+gridv[cd_type]]*locals()['e2'+gridv[cd_type]]
+
+    # rename indices
+    locals()[areaName] = locals()[areaName].rename ({'y':'y_grid_'+cd_type, 'x':'x_grid_'+cd_type})
+
+    # add attributes
+    locals()[areaName].name = areaName
+    locals()[areaName].attrs['standard_name'] = 'cell_area'
+    locals()[areaName].attrs['units']         = 'm2'
+
+    # remove fill values
     locals()[areaName].encoding['_FillValue'] = None
 
-# compute grid cells bounds
+##
+#!!! compute grid bounds !!!
+##
+
+#---------------------------------------------------------------
+# function to generate grid bounds from NEMO coordinates.nc file
+def set_bounds (cdgrd) :
+    '''
+    Constructs lon/lat bounds
+    Bounds are numerated counter clockwise, from bottom left
+    See NEMO file OPA_SRC/IOM/iom.F90, ROUTINE set_grid_bounds, for more details
+    '''
+    # Define offset of coordinate representing bottom-left corner
+    if cdgrd in ['T', 'W'] : 
+        icnr = -1 ; jcnr = -1
+        corner_lon = f_coord ['glamf'].copy() ; corner_lat = f_coord ['gphif'].copy()
+        center_lon = f_coord ['glamt'].copy() ; center_lat = f_coord ['gphit'].copy()
+    if cdgrd == 'U'        : 
+        icnr =  0 ; jcnr = -1
+        corner_lon = f_coord ['glamv'].copy() ; corner_lat = f_coord ['gphiv'].copy()
+        center_lon = f_coord ['glamu'].copy() ; center_lat = f_coord ['gphiu'].copy()
+    if cdgrd == 'V'        : 
+        icnr = -1 ; jcnr =  0
+        corner_lon = f_coord ['glamu'].copy() ; corner_lat = f_coord ['gphiu'].copy()
+        center_lon = f_coord ['glamv'].copy() ; center_lat = f_coord ['gphiv'].copy()
+    if cdgrd == 'F'        : 
+        icnr = -1 ; jcnr = -1
+        corner_lon = f_coord ['glamt'].copy() ; corner_lat = f_coord ['gphit'].copy()
+        center_lon = f_coord ['glamf'].copy() ; center_lat = f_coord ['gphif'].copy()
+      
+    jpj, jpi = corner_lon.shape ;
+    nvertex = 4
+    dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]
+    
+    bounds_lon = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
+    bounds_lat = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
+    
+    idx = [(jcnr,icnr), (jcnr,icnr+1), (jcnr+1,icnr+1), (jcnr+1,icnr)]
+    
+    # Compute cell vertices that can be defined, 
+    # and complete with periodicity
+    for nn in range (nvertex) :
+        tmp = np.roll (corner_lon, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
+        bounds_lon[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
+        tmp = np.roll (corner_lat, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
+        bounds_lat[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
+        bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
+        bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
+    
+    # Zero-size cells at closed boundaries if cell points provided,
+    # otherwise they are closed cells with unrealistic bounds
+    if not (nperio == 1 or nperio == 4 or nperio == 6) :
+        for nn in range (nvertex) : 
+            bounds_lon[:,0,nn]   = center_lon[:,0]   # (West or jpni = 1), closed E-W
+            bounds_lat[:,0,nn]   = center_lat[:,0]
+            bounds_lon[:,jpi,nn] = center_lon[:,jpi] # (East or jpni = 1), closed E-W
+            bounds_lat[:,jpi,nn] = center_lat[:,jpi]
+    if nperio != 2 :
+        for nn in range (nvertex) :
+            bounds_lon[0,:,nn] = center_lon[0,:] # (South or jpnj = 1), not symmetric
+            bounds_lat[0,:,nn] = center_lat[0,:]
+    if nperio < 3 :
+        for nn in range (nvertex) :
+            bounds_lon[jpj,:,nn] = center_lon[jpj,:] # (North or jpnj = 1), no north fold
+            bounds_lat[jpj,:,nn] = center_lat[jpj,:]
+    
+    # Rotate cells at the north fold
+    if nperio >= 3 :
+        # Working array for location of northfold
+        z_fld = nemo.lbc (np.ones ((jpj, jpi)), nperio=nperio, cd_type=cdgrd, psgn=-1., nemo_4U_bug=nemoUbug)
+        z_fld = np.repeat((z_fld == -1.0)[...,np.newaxis],4,axis=2)
+        # circular shift of 2 indices in bounds third index
+        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
+        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
+        bounds_lon[:,:,:]  = np.where (z_fld, bounds_lon_tmp[:,:,:] , bounds_lon[:,:,:] )
+        bounds_lat[:,:,:]  = np.where (z_fld, bounds_lat_tmp[:,:,:] , bounds_lat[:,:,:] )
+    
+    # Invert cells at the symmetric equator
+    if nperio == 2 :
+        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
+        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
+        bounds_lon [0,:,:] = bounds_lon [0,:,:]
+        bounds_lat [0,:,:] = bounds_lat [0,:,:]
+    
+    #bounds_lon.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
+    #bounds_lat.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
+    bounds_lon.attrs['units']       = 'degrees_east'
+    bounds_lat.attrs['units']       = 'degrees_north'
+    bounds_lon.name = 'bounds_lon_grid_' + cdgrd
+    bounds_lat.name = 'bounds_lat_grid_' + cdgrd
+    # remove _FillValue
+    bounds_lon.encoding['_FillValue'] = None
+    bounds_lat.encoding['_FillValue'] = None
+
+    return bounds_lon, bounds_lat
+#------------------------------------------------------------
+
 bounds_lon_grid_T, bounds_lat_grid_T = set_bounds ('T')
 bounds_lon_grid_U, bounds_lat_grid_U = set_bounds ('U')
@@ -494,7 +422,7 @@
 
 #replace nav_lon nav_lat with variables obtained from NEMO coordinates.nc file
-#by construction nav_lon nav_lat are those of the bathymetry
+#by construction nav_lon nav_lat come from the bathymetry
 for cd_type in ['T', 'U', 'V', 'F', 'W'] :
-    for dir_type in ['lat', 'lon'] :
+    for dir_type in ['lon', 'lat'] :
         coord_name = 'nav_' + dir_type + '_grid_' + cd_type
         ds.coords[coord_name]=locals()[coord_name]
