source: TOOLS/MOSAIX/Build_coordinates_mask.py @ 5942

Last change on this file since 5942 was 5942, checked in by snguyen, 3 years ago

add nperio condition for bathymetry southernmost boundary mask
File size: 22.5 KB
Line 
1### ===========================================================================
2###
3### Build coordnates mask
4###
5### ===========================================================================
6# creates file coordinates_mask.nc used by mosaix from NEMO netcdf files
7# coordinates.nc and bathymetry.nc (on the same grid as coordinates.nc)
8
9import numpy  as np
10import xarray as xr
11import nemo
12import datetime, os, platform
13
14# reproduce periodicity lbclnk bug for U and nperio = 4
15nemoUbug = False
16# change metrics (and bathymetry) for straits
17straits = False
18# periodicity imposed on coordinates, metrics and areas
19coordperio = False
20# function used to compute mask_T from bathymetry
21maskbathynemo = False
22
23# type of grid treated : ORCA 1, 2, paleorca # to be improved
24
25#model   = 'eORCA1.2'
26#model   = 'ORCA2.3'
27#model   = 'ORCA2.4'
28model   = 'paleORCA'
29#model   = 'ORCA.05'
30
31if model == 'eORCA1.2' :
32    n_coord = 'eORCA1.2coordinates.nc' #'eORCA_R1_coordinates_v1.0.nc'
33    n_bathy = 'eORCA1.2bathy_meter.nc' #'eORCA_R1_bathy_meter_v2.2.nc'
34    #n_coord = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_coordinates_v1.0.nc'
35    #n_bathy = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_bathy_meter_v2.2.nc'
36    nperio  = 6
37    n_out   = 'eORCA1.2_coordinates_mask_test.nc'
38
39if model == 'ORCA2.3' :
40    nperio  = 4
41    n_coord = 'coordinates_ORCA2.nc' # 'ORCA2.3_coordinates.nc' # 'coordinates_ORCA2_LIM3_PISCES.nc' #
42    n_bathy = 'bathy_meter_ORCA2.nc' # 'ORCA2.3_bathy_meter.nc' # 'bathy_meter_ORCA2_LIM3_PISCES.nc' #
43    n_out   = 'My_coordinates_xios_ORCA2.nc'
44   
45if model == 'ORCA2.4' :
46    nperio  = 4
47    n_coord = 'ORCA2.4_coordinates.nc'
48    n_bathy = 'ORCA2.4_bathy_meter.nc'
49    n_out   = 'ORCA2.4_coordinates_mask_test.nc'
50
51if model == 'paleORCA' :
52    n_coord = 'coordinates_paleo.nc'
53    n_bathy = 'bathy_meter_paleo.nc'
54    nperio = 6
55    n_out = 'paleo_coordinates_mask_test.nc'
56
57#---------------------------------------------------------------
58# function to generate grid bounds from NEMO coordinates.nc file
59def set_bounds (cdgrd) :
60    '''
61    Constructs lon/lat bounds
62    Bounds are numerated counter clockwise, from bottom left
63    See NEMO file OPA_SRC/IOM/iom.F90, ROUTINE set_grid_bounds, for more details
64    '''
65    # Define offset of coordinate representing bottom-left corner
66    if cdgrd in ['T', 'W'] : 
67        icnr = -1 ; jcnr = -1
68        corner_lon = f_coord ['glamf'].copy() ; corner_lat = f_coord ['gphif'].copy()
69        center_lon = f_coord ['glamt'].copy() ; center_lat = f_coord ['gphit'].copy()
70    if cdgrd == 'U'        : 
71        icnr =  0 ; jcnr = -1
72        corner_lon = f_coord ['glamv'].copy() ; corner_lat = f_coord ['gphiv'].copy()
73        center_lon = f_coord ['glamu'].copy() ; center_lat = f_coord ['gphiu'].copy()
74    if cdgrd == 'V'        : 
75        icnr = -1 ; jcnr =  0
76        corner_lon = f_coord ['glamu'].copy() ; corner_lat = f_coord ['gphiu'].copy()
77        center_lon = f_coord ['glamv'].copy() ; center_lat = f_coord ['gphiv'].copy()
78    if cdgrd == 'F'        : 
79        icnr = -1 ; jcnr = -1
80        corner_lon = f_coord ['glamt'].copy() ; corner_lat = f_coord ['gphit'].copy()
81        center_lon = f_coord ['glamf'].copy() ; center_lat = f_coord ['gphif'].copy()
82     
83    #y_grid, x_grid = corner_lon.shape ;
84    nvertex = 4
85    dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]
86   
87    bounds_lon = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
88    bounds_lat = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
89   
90    idx = [(jcnr,icnr), (jcnr,icnr+1), (jcnr+1,icnr+1), (jcnr+1,icnr)]
91   
92    # Compute cell vertices that can be defined,
93    # and complete with periodicity
94    for nn in range (nvertex) :
95        tmp = np.roll (corner_lon, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
96        bounds_lon[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
97        tmp = np.roll (corner_lat, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
98        bounds_lat[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
99        bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
100        bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
101   
102    # Zero-size cells at closed boundaries if cell points provided,
103    # otherwise they are closed cells with unrealistic bounds
104    if not (nperio == 1 or nperio == 4 or nperio == 6) :
105        for nn in range (nvertex) : 
106            bounds_lon[:,0,nn]   = center_lon[:,0]   # (West or jpni = 1), closed E-W
107            bounds_lat[:,0,nn]   = center_lat[:,0]
108            bounds_lon[:,jpi,nn] = center_lon[:,jpi] # (East or jpni = 1), closed E-W
109            bounds_lat[:,jpi,nn] = center_lat[:,jpi]
110    if nperio != 2 :
111        for nn in range (nvertex) :
112            bounds_lon[0,:,nn] = center_lon[0,:] # (South or jpnj = 1), not symmetric
113            bounds_lat[0,:,nn] = center_lat[0,:]
114    if nperio < 3 :
115        for nn in range (nvertex) :
116            bounds_lon[jpj,:,nn] = center_lon[jpj,:] # (North or jpnj = 1), no north fold
117            bounds_lat[jpj,:,nn] = center_lat[jpj,:]
118   
119    # Rotate cells at the north fold
120    if nperio >= 3 :
121        # Working array for location of northfold
122        z_fld = nemo.lbc (np.ones ((jpj, jpi)), nperio=nperio, cd_type=cdgrd, psgn=-1., nemo_4U_bug=nemoUbug)
123        z_fld = np.repeat((z_fld == -1.0)[...,np.newaxis],4,axis=2)
124        # circular shift of 2 indices in bounds third index
125        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
126        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
127        bounds_lon[:,:,:]  = np.where (z_fld, bounds_lon_tmp[:,:,:] , bounds_lon[:,:,:] )
128        bounds_lat[:,:,:]  = np.where (z_fld, bounds_lat_tmp[:,:,:] , bounds_lat[:,:,:] )
129   
130    # Invert cells at the symmetric equator
131    if nperio == 2 :
132        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
133        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
134        bounds_lon [0,:,:] = bounds_lon [0,:,:]
135        bounds_lat [0,:,:] = bounds_lat [0,:,:]
136   
137    #bounds_lon.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
138    #bounds_lat.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
139    bounds_lon.attrs['units']       = 'degrees_east'
140    bounds_lat.attrs['units']       = 'degrees_north'
141    bounds_lon.name = 'bounds_lon_grid_' + cdgrd
142    bounds_lat.name = 'bounds_lat_grid_' + cdgrd
143    # remove _FillValue
144    bounds_lon.encoding['_FillValue'] = None
145    bounds_lat.encoding['_FillValue'] = None
146
147    return bounds_lon, bounds_lat
148#------------------------------------------------------------
149
150# open input files while removing the time dimension
151f_coord = xr.open_dataset (n_coord, decode_times=False).squeeze()
152f_bathy = xr.open_dataset (n_bathy, decode_times=False).squeeze()
153
154# Suppress time if necessary
155try    :
156    del f_coord['time']
157    print ('time successfully removed')
158except :
159    pass
160    print ('failed to suppress time')
161
162Bathymetry = f_bathy['Bathymetry'].copy()
163
164try :
165    Bathymetry = Bathymetry.rename ({'nav_lat':'nav_lat_grid_T', 'nav_lon':'nav_lon_grid_T'})
166    print ('Bathymetry file: Normal case')
167except : 
168    print ('Bathymetry file: Exception')
169    nav_lon_grid_T = f_bathy ['nav_lon']
170    nav_lat_grid_T = f_bathy ['nav_lat']
171    Bathymetry = xr.DataArray (Bathymetry, coords = { "nav_lat_grid_T": (["y", "x"], nav_lat_grid_T),
172                                                      "nav_lon_grid_T": (["y", "x"], nav_lon_grid_T) } )
173   
174Bathymetry = Bathymetry.rename ({'y':'y_grid_T', 'x':'x_grid_T'})
175
176if straits :
177    # Open straits for usual grids
178    if model == 'ORCA2.3' :
179        # orca_r2: Gibraltar strait open
180        Bathymetry[101,139] = 284.
181        # orca_r2: Bab el Mandeb strait open
182        Bathymetry[87,159] = 137.
183
184# Determine periodicity from grid type
185jpj, jpi = Bathymetry.shape
186try : 
187    if nperio != None :
188        print ('nperio specified : ', nperio)
189except :
190    nperio = None
191    if jpi ==  182 : nperio = 4 # ORCA2. \!/ We choose legacy orca2
192    if jpi ==  362 : nperio = 6 # ORCA1.
193    if jpi == 1442 : nperio = 6 # ORCA025.
194    #
195    if nperio == None :
196        sys.exit ('nperio not found, and cannot by guessed' )
197    else :
198        print ('nperio set as {:d} (deduced from data size {:d} {:d})'.format (nperio, jpj, jpi))
199
200Bathymetry = nemo.lbc (Bathymetry, nperio=nperio, cd_type='T')
201#Suppress ocean points at southernmost position only if nperio in {3, 4, 5, 6}
202if nperio in (3, 4, 5, 6) :
203    Bathymetry[0,:] = 0.0 
204
205# Creation of mask_T
206
207if maskbathynemo :
208    # Use same formula as domzgr.
209    mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='int8')
210else :
211    mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='int8')
212
213# Creation of U, V, W, F masks from mask_T
214mask_U = mask_T * mask_T.shift (x_grid_T=-1)
215mask_V = mask_T * mask_T.shift (y_grid_T=-1)
216mask_F = mask_T * mask_T.shift (y_grid_T=-1) * mask_T.shift (x_grid_T=-1) * mask_T.shift (y_grid_T=-1, x_grid_T=-1)
217mask_W = mask_T
218
219mask_U = nemo.lbc (mask_U, nperio=nperio, cd_type='U', nemo_4U_bug=nemoUbug).astype (dtype='int8')
220mask_V = nemo.lbc (mask_V, nperio=nperio, cd_type='V').astype (dtype='int8')
221mask_F = nemo.lbc (mask_F, nperio=nperio, cd_type='F').astype (dtype='int8')
222mask_W = nemo.lbc (mask_W, nperio=nperio, cd_type='W').astype (dtype='int8')
223
224mask_U = mask_U.rename ( {'y_grid_T'      : 'y_grid_U'      , 'x_grid_T'      : 'x_grid_U', 
225                          'nav_lat_grid_T': 'nav_lat_grid_U', 'nav_lon_grid_T': 'nav_lon_grid_U'} )
226mask_V = mask_V.rename ( {'y_grid_T'      : 'y_grid_V'      , 'x_grid_T'      : 'x_grid_V',
227                          'nav_lat_grid_T': 'nav_lat_grid_V', 'nav_lon_grid_T': 'nav_lon_grid_V'} )
228mask_W = mask_W.rename ( {'y_grid_T'      : 'y_grid_W'      , 'x_grid_T'      : 'x_grid_W',
229                          'nav_lat_grid_T': 'nav_lat_grid_W', 'nav_lon_grid_T': 'nav_lon_grid_W'} )
230mask_F = mask_F.rename ( {'y_grid_T'      : 'y_grid_F'      , 'x_grid_T'      : 'x_grid_F',
231                          'nav_lat_grid_T': 'nav_lat_grid_F', 'nav_lon_grid_T': 'nav_lon_grid_F'} )
232
233mask_T.name = 'mask_T'
234mask_U.name = 'mask_U'
235mask_V.name = 'mask_V'
236mask_F.name = 'mask_F'
237mask_W.name = 'mask_W'
238
239# create masks without duplicate points : maskutil_[TUVWF]
240for cd_type in ['T', 'U', 'V', 'F', 'W'] :
241    MaskName = 'mask_'     + cd_type
242    UtilName = 'maskutil_' + cd_type
243    locals()[UtilName] = nemo.lbc_mask (locals()[MaskName].copy(), nperio=nperio, cd_type=cd_type)
244    locals()[UtilName].name = UtilName
245
246# add masks attributes
247mask_T.attrs    ['cell_measures'] = 'area: area_grid_T'
248mask_U.attrs    ['cell_measures'] = 'area: area_grid_U'
249mask_V.attrs    ['cell_measures'] = 'area: area_grid_V'
250mask_W.attrs    ['cell_measures'] = 'area: area_grid_W'
251mask_F.attrs    ['cell_measures'] = 'area: area_grid_F'
252
253maskutil_T.attrs['cell_measures'] = 'area: area_grid_T'
254maskutil_U.attrs['cell_measures'] = 'area: area_grid_U'
255maskutil_V.attrs['cell_measures'] = 'area: area_grid_V'
256maskutil_W.attrs['cell_measures'] = 'area: area_grid_W'
257maskutil_F.attrs['cell_measures'] = 'area: area_grid_F'
258
259# create grid coordinates from NEMO coordinates.nc file
260if coordperio :
261    nav_lon_grid_T = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='T')
262    nav_lat_grid_T = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='T')
263    nav_lon_grid_U = nemo.lbc (f_coord ['glamu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
264    nav_lat_grid_U = nemo.lbc (f_coord ['gphiu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
265    nav_lon_grid_V = nemo.lbc (f_coord ['glamv'].copy(), nperio=nperio, cd_type='V')
266    nav_lat_grid_V = nemo.lbc (f_coord ['gphiv'].copy(), nperio=nperio, cd_type='V')
267    nav_lon_grid_W = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='W')
268    nav_lat_grid_W = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='W')
269    nav_lon_grid_F = nemo.lbc (f_coord ['glamf'].copy(), nperio=nperio, cd_type='F')
270    nav_lat_grid_F = nemo.lbc (f_coord ['gphif'].copy(), nperio=nperio, cd_type='F')
271else :
272    nav_lon_grid_T = f_coord ['glamt'].copy()
273    nav_lat_grid_T = f_coord ['gphit'].copy()
274    nav_lon_grid_U = f_coord ['glamu'].copy()
275    nav_lat_grid_U = f_coord ['gphiu'].copy()
276    nav_lon_grid_V = f_coord ['glamv'].copy()
277    nav_lat_grid_V = f_coord ['gphiv'].copy()
278    nav_lon_grid_W = f_coord ['glamt'].copy()
279    nav_lat_grid_W = f_coord ['gphit'].copy()
280    nav_lon_grid_F = f_coord ['glamf'].copy()
281    nav_lat_grid_F = f_coord ['gphif'].copy()
282
283nav_lon_grid_T = nav_lon_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )
284nav_lat_grid_T = nav_lat_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )
285nav_lon_grid_U = nav_lon_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )
286nav_lat_grid_U = nav_lat_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )
287nav_lon_grid_V = nav_lon_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )
288nav_lat_grid_V = nav_lat_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )
289nav_lon_grid_W = nav_lon_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )
290nav_lat_grid_W = nav_lat_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )
291nav_lon_grid_F = nav_lon_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )
292nav_lat_grid_F = nav_lat_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )
293
294# remove _FillValue and missing_value
295for cd_type in ['T', 'U', 'V', 'F', 'W'] :
296    for dir_type in ['lat', 'lon'] :
297        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
298        locals()[coord_name].encoding['_FillValue'] = None
299        locals()[coord_name].encoding['missing_value'] = None
300
301nav_lon_grid_T.name = 'nav_lon_grid_T'
302nav_lat_grid_T.name = 'nav_lat_grid_T'
303nav_lon_grid_U.name = 'nav_lon_grid_U'
304nav_lat_grid_U.name = 'nav_lat_grid_U'
305nav_lon_grid_V.name = 'nav_lon_grid_V'
306nav_lat_grid_V.name = 'nav_lat_grid_V'
307nav_lon_grid_F.name = 'nav_lon_grid_F'
308nav_lat_grid_F.name = 'nav_lat_grid_F'
309nav_lon_grid_W.name = 'nav_lon_grid_W'
310nav_lat_grid_W.name = 'nav_lat_grid_W'
311
312# add coordinates attributes
313nav_lon_grid_T.attrs['standard_name'] = 'longitude'
314nav_lon_grid_T.attrs['long_name']     = 'Longitude'
315nav_lon_grid_T.attrs['units']         = 'degrees_east'
316nav_lat_grid_T.attrs['standard_name'] = 'latitude'
317nav_lat_grid_T.attrs['long_name']     = 'Latitude'
318nav_lat_grid_T.attrs['units']         = 'degrees_north'
319
320nav_lon_grid_U.attrs['standard_name'] = 'longitude'
321nav_lon_grid_U.attrs['long_name']     = 'Longitude'
322nav_lon_grid_U.attrs['units']         = 'degrees_east'
323nav_lat_grid_U.attrs['standard_name'] = 'latitude'
324nav_lat_grid_U.attrs['long_name']     = 'Latitude'
325nav_lat_grid_U.attrs['units']         = 'degrees_north'
326
327nav_lon_grid_V.attrs['standard_name'] = 'longitude'
328nav_lon_grid_V.attrs['long_name']     = 'Longitude'
329nav_lon_grid_V.attrs['units']         = 'degrees_east'
330nav_lat_grid_V.attrs['standard_name'] = 'latitude'
331nav_lat_grid_V.attrs['long_name']     = 'Latitude'
332nav_lat_grid_V.attrs['units']         = 'degrees_north'
333
334nav_lon_grid_W.attrs['standard_name'] = 'longitude'
335nav_lon_grid_W.attrs['long_name']     = 'Longitude'
336nav_lon_grid_W.attrs['units']         = 'degrees_east'
337nav_lat_grid_W.attrs['standard_name'] = 'latitude'
338nav_lat_grid_W.attrs['long_name']     = 'Latitude'
339nav_lat_grid_W.attrs['units']         = 'degrees_north'
340
341nav_lon_grid_F.attrs['standard_name'] = 'longitude'
342nav_lon_grid_F.attrs['long_name']     = 'Longitude'
343nav_lon_grid_F.attrs['units']         = 'degrees_east'
344nav_lat_grid_F.attrs['standard_name'] = 'latitude'
345nav_lat_grid_F.attrs['long_name']     = 'Latitude'
346nav_lat_grid_F.attrs['units']         = 'degrees_north'
347
348nav_lon_grid_T.attrs['bounds'] = 'bounds_lon_grid_T'
349nav_lat_grid_T.attrs['bounds'] = 'bounds_lat_grid_T'
350nav_lon_grid_U.attrs['bounds'] = 'bounds_lon_grid_U'
351nav_lat_grid_U.attrs['bounds'] = 'bounds_lat_grid_U'
352nav_lon_grid_V.attrs['bounds'] = 'bounds_lon_grid_V'
353nav_lat_grid_V.attrs['bounds'] = 'bounds_lat_grid_V'
354nav_lon_grid_W.attrs['bounds'] = 'bounds_lon_grid_W'
355nav_lat_grid_W.attrs['bounds'] = 'bounds_lat_grid_W'
356nav_lon_grid_F.attrs['bounds'] = 'bounds_lon_grid_F'
357nav_lat_grid_F.attrs['bounds'] = 'bounds_lat_grid_F'
358
359# create areas variables for grid cells from NEMO coordinates.nc file
360# create new variables e1 e2 to keep f_coord the same
361for cd_type in ['t', 'u', 'v', 'f'] :
362    for axis in ['1', '2'] :
363        coordName = 'e' + axis + cd_type
364        locals()[coordName]=f_coord[coordName].copy()
365        # remove zero values from areas
366        # need to be define for the extended grid south of -80S
367        # some point are undefined but you need to have e1 and e2 .NE. 0
368        locals()[coordName]=xr.where(locals()[coordName] == 0.0, 1.0e2, locals()[coordName])
369
370#Compute cells areas
371
372# Correct areas for straits
373if straits :
374    # ORCA R2 configuration
375    if model == 'ORCA2.3' :
376        # Gibraltar     : e2u reduced to 20 km
377        e2u[101,138:140] = 20.e3
378        # Bab el Mandeb : e2u reduced to 30 km
379        #                e1v reduced to 18 km
380        e1v[87,159] = 18.e3
381        e2u[87,159] = 30.e3
382        # Danish Straits: e2u reduced to 10 km
383        e2u[115,144:146] = 10.e3
384    # ORCA R1 configuration
385    if model == 'eORCA1.2' :
386        # Gibraltar : e2u reduced to 20 km
387        e2u[240,281:283] = 20.e3
388        # Bhosporus : e2u reduced to 10 km
389        e2u[247,313:315] =  10.e3
390        # Lombok : e1v reduced to 13 km
391        e1v[163:165,43] =  13.e3
392        # Sumba : e1v reduced to 8 km
393        e1v[163:165,47] =  8.e3
394        # Ombai : e1v reduced to 13 km
395        e1v[163:165,52] = 13.e3
396        # Timor Passage : e1v reduced to 20 km
397        #e1v[163:165,55] = 20.e3
398        # W Halmahera : e1v reduced to 30 km
399        e1v[180:182,54] = 30.e3
400        # E Halmahera : e1v reduced to 50 km
401        e1v[180:182,57] = 50.e3
402    # ORCA R05 configuration
403    if model == 'ORCA.05' :
404        # Reduced e2u at the Gibraltar Strait
405        e2u[326,562:564] =  20.e3
406        # Reduced e2u at the Bosphore Strait
407        e2u[342,626:628] =  10.e3
408        # Reduced e2u at the Sumba Strait
409        e2u[231,92:94] =  40.e3
410        # Reduced e2u at the Ombai Strait
411        e2u[231,102] =  15.e3
412        # Reduced e2u at the Palk Strait
413        e2u[269,14] =  10.e3
414        # Reduced e1v at the Lombok Strait
415        e1v[231:233,86] =  10.e3
416        # Reduced e1v at the Bab el Mandeb
417        e1v[275,661] =  25.e3
418
419if coordperio :
420    area_grid_T = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='T')
421    area_grid_U = nemo.lbc (e1u*e2u, nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
422    area_grid_V = nemo.lbc (e1v*e2v, nperio=nperio, cd_type='V')
423    area_grid_W = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='W')
424    area_grid_F = nemo.lbc (e1f*e2f, nperio=nperio, cd_type='F')
425else :
426    area_grid_T = e1t*e2t
427    area_grid_U = e1u*e2u
428    area_grid_V = e1v*e2v
429    area_grid_W = e1t*e2t
430    area_grid_F = e1f*e2f
431
432area_grid_T.name = 'area_grid_T'
433area_grid_U.name = 'area_grid_U'
434area_grid_V.name = 'area_grid_V'
435area_grid_F.name = 'area_grid_F'
436area_grid_W.name = 'area_grid_W'
437
438area_grid_T = area_grid_T.rename ({'y':'y_grid_T', 'x':'x_grid_T'})
439area_grid_U = area_grid_U.rename ({'y':'y_grid_U', 'x':'x_grid_U'})
440area_grid_V = area_grid_V.rename ({'y':'y_grid_V', 'x':'x_grid_V'})
441area_grid_W = area_grid_W.rename ({'y':'y_grid_W', 'x':'x_grid_W'})
442area_grid_F = area_grid_F.rename ({'y':'y_grid_F', 'x':'x_grid_F'})
443
444area_grid_T.attrs['standard_name'] = 'cell_area'
445area_grid_T.attrs['units']         = 'm2'
446area_grid_U.attrs['standard_name'] = 'cell_area'
447area_grid_U.attrs['units']         = 'm2'
448area_grid_V.attrs['standard_name'] = 'cell_area'
449area_grid_V.attrs['units']         = 'm2'
450area_grid_W.attrs['standard_name'] = 'cell_area'
451area_grid_W.attrs['units']         = 'm2'
452area_grid_F.attrs['standard_name'] = 'cell_area'
453area_grid_F.attrs['units']         = 'm2'
454
455for cd_type in ['T', 'U', 'V', 'F', 'W'] :
456    areaName = 'area_grid_' + cd_type
457    locals()[areaName].encoding['_FillValue'] = None
458
459# compute grid cells bounds
460bounds_lon_grid_T, bounds_lat_grid_T = set_bounds ('T')
461bounds_lon_grid_U, bounds_lat_grid_U = set_bounds ('U')
462bounds_lon_grid_V, bounds_lat_grid_V = set_bounds ('V')
463bounds_lon_grid_W, bounds_lat_grid_W = set_bounds ('W')
464bounds_lon_grid_F, bounds_lat_grid_F = set_bounds ('F')
465
466# build xarray dataset to be saved
467ds = xr.Dataset ({
468    'mask_T'      : mask_T     ,
469    'mask_U'      : mask_U     ,
470    'mask_V'      : mask_V     ,
471    'mask_W'      : mask_W     ,
472    'mask_F'      : mask_F     ,
473    'area_grid_T' : area_grid_T,
474    'area_grid_U' : area_grid_U,
475    'area_grid_V' : area_grid_V,
476    'area_grid_W' : area_grid_W,
477    'area_grid_F' : area_grid_F,
478    'maskutil_T'  : maskutil_T ,
479    'maskutil_U'  : maskutil_U ,
480    'maskutil_V'  : maskutil_V ,
481    'maskutil_W'  : maskutil_W ,
482    'maskutil_F'  : maskutil_F ,
483    'bounds_lon_grid_T': bounds_lon_grid_T,
484    'bounds_lat_grid_T': bounds_lat_grid_T,
485    'bounds_lon_grid_U': bounds_lon_grid_U,
486    'bounds_lat_grid_U': bounds_lat_grid_U,
487    'bounds_lon_grid_V': bounds_lon_grid_V,
488    'bounds_lat_grid_V': bounds_lat_grid_V,
489    'bounds_lon_grid_W': bounds_lon_grid_W,
490    'bounds_lat_grid_W': bounds_lat_grid_W,
491    'bounds_lon_grid_F': bounds_lon_grid_F,
492    'bounds_lat_grid_F': bounds_lat_grid_F,
493})
494
495#replace nav_lon nav_lat with variables obtained from NEMO coordinates.nc file
496#by construction nav_lon nav_lat are those of the bathymetry
497for cd_type in ['T', 'U', 'V', 'F', 'W'] :
498    for dir_type in ['lat', 'lon'] :
499        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
500        ds.coords[coord_name]=locals()[coord_name]
501
502ds.attrs['name']         = 'coordinates_mask'
503ds.attrs['description']  = 'coordinates and mask for MOSAIX'
504ds.attrs['title']        = 'coordinates_mask'
505ds.attrs['source']       = 'IPSL Earth system model'
506ds.attrs['group']        = 'ICMC IPSL Climate Modelling Center'
507ds.attrs['Institution']  = 'IPSL https.//www.ipsl.fr'
508ds.attrs['Model']        = model
509ds.attrs['timeStamp']    = '{:%Y-%b-%d %H:%M:%S}'.format (datetime.datetime.now ())
510ds.attrs['history']      = 'Build from ' + n_coord + ' and ' + n_bathy
511ds.attrs['directory']    = os.getcwd     ()
512ds.attrs['user']         = os.getlogin   ()
513ds.attrs['HOSTNAME']     = platform.node ()
514ds.attrs['Python']       = 'Python version: ' +  platform.python_version ()
515ds.attrs['xarray']       = 'xarray version: ' +  xr.__version__
516ds.attrs['OS']           = platform.system  ()
517ds.attrs['release']      = platform.release ()
518ds.attrs['hardware']     = platform.machine ()
519
520ds.attrs['SVN_Author']   = "$Author:  $"
521ds.attrs['SVN_Date']     = "$Date:  $"
522ds.attrs['SVN_Revision'] = "$Revision:  $"
523ds.attrs['SVN_Id']       = "$Id: Build_coordinates_mask.py $"
524ds.attrs['SVN_HeadURL']  = "$HeadURL: http://forge.ipsl.jussieu.fr/igcmg/svn/TOOLS/MOSAIX/Build_coordinates_mask.py $"
525
526# save to output file
527ds.to_netcdf (n_out)
Note: See TracBrowser for help on using the repository browser.