source: TOOLS/MOSAIX/Build_coordinates_mask.py @ 5940

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

first script version of iPython Notebook to build coordinates, masks, areas and bounds...
File size: 22.4 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
202Bathymetry[0,:] = 0.0 
203
204# Creation of mask_T
205
206if maskbathynemo :
207    # Use same formula as domzgr.
208    mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='int8')
209else :
210    mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='int8')
211
212# Creation of U, V, W, F masks from mask_T
213mask_U = mask_T * mask_T.shift (x_grid_T=-1)
214mask_V = mask_T * mask_T.shift (y_grid_T=-1)
215mask_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)
216mask_W = mask_T
217
218mask_U = nemo.lbc (mask_U, nperio=nperio, cd_type='U', nemo_4U_bug=nemoUbug).astype (dtype='int8')
219mask_V = nemo.lbc (mask_V, nperio=nperio, cd_type='V').astype (dtype='int8')
220mask_F = nemo.lbc (mask_F, nperio=nperio, cd_type='F').astype (dtype='int8')
221mask_W = nemo.lbc (mask_W, nperio=nperio, cd_type='W').astype (dtype='int8')
222
223mask_U = mask_U.rename ( {'y_grid_T'      : 'y_grid_U'      , 'x_grid_T'      : 'x_grid_U', 
224                          'nav_lat_grid_T': 'nav_lat_grid_U', 'nav_lon_grid_T': 'nav_lon_grid_U'} )
225mask_V = mask_V.rename ( {'y_grid_T'      : 'y_grid_V'      , 'x_grid_T'      : 'x_grid_V',
226                          'nav_lat_grid_T': 'nav_lat_grid_V', 'nav_lon_grid_T': 'nav_lon_grid_V'} )
227mask_W = mask_W.rename ( {'y_grid_T'      : 'y_grid_W'      , 'x_grid_T'      : 'x_grid_W',
228                          'nav_lat_grid_T': 'nav_lat_grid_W', 'nav_lon_grid_T': 'nav_lon_grid_W'} )
229mask_F = mask_F.rename ( {'y_grid_T'      : 'y_grid_F'      , 'x_grid_T'      : 'x_grid_F',
230                          'nav_lat_grid_T': 'nav_lat_grid_F', 'nav_lon_grid_T': 'nav_lon_grid_F'} )
231
232mask_T.name = 'mask_T'
233mask_U.name = 'mask_U'
234mask_V.name = 'mask_V'
235mask_F.name = 'mask_F'
236mask_W.name = 'mask_W'
237
238# create masks without duplicate points : maskutil_[TUVWF]
239for cd_type in ['T', 'U', 'V', 'F', 'W'] :
240    MaskName = 'mask_'     + cd_type
241    UtilName = 'maskutil_' + cd_type
242    locals()[UtilName] = nemo.lbc_mask (locals()[MaskName].copy(), nperio=nperio, cd_type=cd_type)
243    locals()[UtilName].name = UtilName
244
245# add masks attributes
246mask_T.attrs    ['cell_measures'] = 'area: area_grid_T'
247mask_U.attrs    ['cell_measures'] = 'area: area_grid_U'
248mask_V.attrs    ['cell_measures'] = 'area: area_grid_V'
249mask_W.attrs    ['cell_measures'] = 'area: area_grid_W'
250mask_F.attrs    ['cell_measures'] = 'area: area_grid_F'
251
252maskutil_T.attrs['cell_measures'] = 'area: area_grid_T'
253maskutil_U.attrs['cell_measures'] = 'area: area_grid_U'
254maskutil_V.attrs['cell_measures'] = 'area: area_grid_V'
255maskutil_W.attrs['cell_measures'] = 'area: area_grid_W'
256maskutil_F.attrs['cell_measures'] = 'area: area_grid_F'
257
258# create grid coordinates from NEMO coordinates.nc file
259if coordperio :
260    nav_lon_grid_T = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='T')
261    nav_lat_grid_T = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='T')
262    nav_lon_grid_U = nemo.lbc (f_coord ['glamu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
263    nav_lat_grid_U = nemo.lbc (f_coord ['gphiu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
264    nav_lon_grid_V = nemo.lbc (f_coord ['glamv'].copy(), nperio=nperio, cd_type='V')
265    nav_lat_grid_V = nemo.lbc (f_coord ['gphiv'].copy(), nperio=nperio, cd_type='V')
266    nav_lon_grid_W = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='W')
267    nav_lat_grid_W = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='W')
268    nav_lon_grid_F = nemo.lbc (f_coord ['glamf'].copy(), nperio=nperio, cd_type='F')
269    nav_lat_grid_F = nemo.lbc (f_coord ['gphif'].copy(), nperio=nperio, cd_type='F')
270else :
271    nav_lon_grid_T = f_coord ['glamt'].copy()
272    nav_lat_grid_T = f_coord ['gphit'].copy()
273    nav_lon_grid_U = f_coord ['glamu'].copy()
274    nav_lat_grid_U = f_coord ['gphiu'].copy()
275    nav_lon_grid_V = f_coord ['glamv'].copy()
276    nav_lat_grid_V = f_coord ['gphiv'].copy()
277    nav_lon_grid_W = f_coord ['glamt'].copy()
278    nav_lat_grid_W = f_coord ['gphit'].copy()
279    nav_lon_grid_F = f_coord ['glamf'].copy()
280    nav_lat_grid_F = f_coord ['gphif'].copy()
281
282nav_lon_grid_T = nav_lon_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )
283nav_lat_grid_T = nav_lat_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )
284nav_lon_grid_U = nav_lon_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )
285nav_lat_grid_U = nav_lat_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )
286nav_lon_grid_V = nav_lon_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )
287nav_lat_grid_V = nav_lat_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )
288nav_lon_grid_W = nav_lon_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )
289nav_lat_grid_W = nav_lat_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )
290nav_lon_grid_F = nav_lon_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )
291nav_lat_grid_F = nav_lat_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )
292
293# remove _FillValue and missing_value
294for cd_type in ['T', 'U', 'V', 'F', 'W'] :
295    for dir_type in ['lat', 'lon'] :
296        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
297        locals()[coord_name].encoding['_FillValue'] = None
298        locals()[coord_name].encoding['missing_value'] = None
299
300nav_lon_grid_T.name = 'nav_lon_grid_T'
301nav_lat_grid_T.name = 'nav_lat_grid_T'
302nav_lon_grid_U.name = 'nav_lon_grid_U'
303nav_lat_grid_U.name = 'nav_lat_grid_U'
304nav_lon_grid_V.name = 'nav_lon_grid_V'
305nav_lat_grid_V.name = 'nav_lat_grid_V'
306nav_lon_grid_F.name = 'nav_lon_grid_F'
307nav_lat_grid_F.name = 'nav_lat_grid_F'
308nav_lon_grid_W.name = 'nav_lon_grid_W'
309nav_lat_grid_W.name = 'nav_lat_grid_W'
310
311# add coordinates attributes
312nav_lon_grid_T.attrs['standard_name'] = 'longitude'
313nav_lon_grid_T.attrs['long_name']     = 'Longitude'
314nav_lon_grid_T.attrs['units']         = 'degrees_east'
315nav_lat_grid_T.attrs['standard_name'] = 'latitude'
316nav_lat_grid_T.attrs['long_name']     = 'Latitude'
317nav_lat_grid_T.attrs['units']         = 'degrees_north'
318
319nav_lon_grid_U.attrs['standard_name'] = 'longitude'
320nav_lon_grid_U.attrs['long_name']     = 'Longitude'
321nav_lon_grid_U.attrs['units']         = 'degrees_east'
322nav_lat_grid_U.attrs['standard_name'] = 'latitude'
323nav_lat_grid_U.attrs['long_name']     = 'Latitude'
324nav_lat_grid_U.attrs['units']         = 'degrees_north'
325
326nav_lon_grid_V.attrs['standard_name'] = 'longitude'
327nav_lon_grid_V.attrs['long_name']     = 'Longitude'
328nav_lon_grid_V.attrs['units']         = 'degrees_east'
329nav_lat_grid_V.attrs['standard_name'] = 'latitude'
330nav_lat_grid_V.attrs['long_name']     = 'Latitude'
331nav_lat_grid_V.attrs['units']         = 'degrees_north'
332
333nav_lon_grid_W.attrs['standard_name'] = 'longitude'
334nav_lon_grid_W.attrs['long_name']     = 'Longitude'
335nav_lon_grid_W.attrs['units']         = 'degrees_east'
336nav_lat_grid_W.attrs['standard_name'] = 'latitude'
337nav_lat_grid_W.attrs['long_name']     = 'Latitude'
338nav_lat_grid_W.attrs['units']         = 'degrees_north'
339
340nav_lon_grid_F.attrs['standard_name'] = 'longitude'
341nav_lon_grid_F.attrs['long_name']     = 'Longitude'
342nav_lon_grid_F.attrs['units']         = 'degrees_east'
343nav_lat_grid_F.attrs['standard_name'] = 'latitude'
344nav_lat_grid_F.attrs['long_name']     = 'Latitude'
345nav_lat_grid_F.attrs['units']         = 'degrees_north'
346
347nav_lon_grid_T.attrs['bounds'] = 'bounds_lon_grid_T'
348nav_lat_grid_T.attrs['bounds'] = 'bounds_lat_grid_T'
349nav_lon_grid_U.attrs['bounds'] = 'bounds_lon_grid_U'
350nav_lat_grid_U.attrs['bounds'] = 'bounds_lat_grid_U'
351nav_lon_grid_V.attrs['bounds'] = 'bounds_lon_grid_V'
352nav_lat_grid_V.attrs['bounds'] = 'bounds_lat_grid_V'
353nav_lon_grid_W.attrs['bounds'] = 'bounds_lon_grid_W'
354nav_lat_grid_W.attrs['bounds'] = 'bounds_lat_grid_W'
355nav_lon_grid_F.attrs['bounds'] = 'bounds_lon_grid_F'
356nav_lat_grid_F.attrs['bounds'] = 'bounds_lat_grid_F'
357
358# create areas variables for grid cells from NEMO coordinates.nc file
359# create new variables e1 e2 to keep f_coord the same
360for cd_type in ['t', 'u', 'v', 'f'] :
361    for axis in ['1', '2'] :
362        coordName = 'e' + axis + cd_type
363        locals()[coordName]=f_coord[coordName].copy()
364        # remove zero values from areas
365        # need to be define for the extended grid south of -80S
366        # some point are undefined but you need to have e1 and e2 .NE. 0
367        locals()[coordName]=xr.where(locals()[coordName] == 0.0, 1.0e2, locals()[coordName])
368
369#Compute cells areas
370
371# Correct areas for straits
372if straits :
373    # ORCA R2 configuration
374    if model == 'ORCA2.3' :
375        # Gibraltar     : e2u reduced to 20 km
376        e2u[101,138:140] = 20.e3
377        # Bab el Mandeb : e2u reduced to 30 km
378        #                e1v reduced to 18 km
379        e1v[87,159] = 18.e3
380        e2u[87,159] = 30.e3
381        # Danish Straits: e2u reduced to 10 km
382        e2u[115,144:146] = 10.e3
383    # ORCA R1 configuration
384    if model == 'eORCA1.2' :
385        # Gibraltar : e2u reduced to 20 km
386        e2u[240,281:283] = 20.e3
387        # Bhosporus : e2u reduced to 10 km
388        e2u[247,313:315] =  10.e3
389        # Lombok : e1v reduced to 13 km
390        e1v[163:165,43] =  13.e3
391        # Sumba : e1v reduced to 8 km
392        e1v[163:165,47] =  8.e3
393        # Ombai : e1v reduced to 13 km
394        e1v[163:165,52] = 13.e3
395        # Timor Passage : e1v reduced to 20 km
396        #e1v[163:165,55] = 20.e3
397        # W Halmahera : e1v reduced to 30 km
398        e1v[180:182,54] = 30.e3
399        # E Halmahera : e1v reduced to 50 km
400        e1v[180:182,57] = 50.e3
401    # ORCA R05 configuration
402    if model == 'ORCA.05' :
403        # Reduced e2u at the Gibraltar Strait
404        e2u[326,562:564] =  20.e3
405        # Reduced e2u at the Bosphore Strait
406        e2u[342,626:628] =  10.e3
407        # Reduced e2u at the Sumba Strait
408        e2u[231,92:94] =  40.e3
409        # Reduced e2u at the Ombai Strait
410        e2u[231,102] =  15.e3
411        # Reduced e2u at the Palk Strait
412        e2u[269,14] =  10.e3
413        # Reduced e1v at the Lombok Strait
414        e1v[231:233,86] =  10.e3
415        # Reduced e1v at the Bab el Mandeb
416        e1v[275,661] =  25.e3
417
418if coordperio :
419    area_grid_T = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='T')
420    area_grid_U = nemo.lbc (e1u*e2u, nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)
421    area_grid_V = nemo.lbc (e1v*e2v, nperio=nperio, cd_type='V')
422    area_grid_W = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='W')
423    area_grid_F = nemo.lbc (e1f*e2f, nperio=nperio, cd_type='F')
424else :
425    area_grid_T = e1t*e2t
426    area_grid_U = e1u*e2u
427    area_grid_V = e1v*e2v
428    area_grid_W = e1t*e2t
429    area_grid_F = e1f*e2f
430
431area_grid_T.name = 'area_grid_T'
432area_grid_U.name = 'area_grid_U'
433area_grid_V.name = 'area_grid_V'
434area_grid_F.name = 'area_grid_F'
435area_grid_W.name = 'area_grid_W'
436
437area_grid_T = area_grid_T.rename ({'y':'y_grid_T', 'x':'x_grid_T'})
438area_grid_U = area_grid_U.rename ({'y':'y_grid_U', 'x':'x_grid_U'})
439area_grid_V = area_grid_V.rename ({'y':'y_grid_V', 'x':'x_grid_V'})
440area_grid_W = area_grid_W.rename ({'y':'y_grid_W', 'x':'x_grid_W'})
441area_grid_F = area_grid_F.rename ({'y':'y_grid_F', 'x':'x_grid_F'})
442
443area_grid_T.attrs['standard_name'] = 'cell_area'
444area_grid_T.attrs['units']         = 'm2'
445area_grid_U.attrs['standard_name'] = 'cell_area'
446area_grid_U.attrs['units']         = 'm2'
447area_grid_V.attrs['standard_name'] = 'cell_area'
448area_grid_V.attrs['units']         = 'm2'
449area_grid_W.attrs['standard_name'] = 'cell_area'
450area_grid_W.attrs['units']         = 'm2'
451area_grid_F.attrs['standard_name'] = 'cell_area'
452area_grid_F.attrs['units']         = 'm2'
453
454for cd_type in ['T', 'U', 'V', 'F', 'W'] :
455    areaName = 'area_grid_' + cd_type
456    locals()[areaName].encoding['_FillValue'] = None
457
458# compute grid cells bounds
459bounds_lon_grid_T, bounds_lat_grid_T = set_bounds ('T')
460bounds_lon_grid_U, bounds_lat_grid_U = set_bounds ('U')
461bounds_lon_grid_V, bounds_lat_grid_V = set_bounds ('V')
462bounds_lon_grid_W, bounds_lat_grid_W = set_bounds ('W')
463bounds_lon_grid_F, bounds_lat_grid_F = set_bounds ('F')
464
465# build xarray dataset to be saved
466ds = xr.Dataset ({
467    'mask_T'      : mask_T     ,
468    'mask_U'      : mask_U     ,
469    'mask_V'      : mask_V     ,
470    'mask_W'      : mask_W     ,
471    'mask_F'      : mask_F     ,
472    'area_grid_T' : area_grid_T,
473    'area_grid_U' : area_grid_U,
474    'area_grid_V' : area_grid_V,
475    'area_grid_W' : area_grid_W,
476    'area_grid_F' : area_grid_F,
477    'maskutil_T'  : maskutil_T ,
478    'maskutil_U'  : maskutil_U ,
479    'maskutil_V'  : maskutil_V ,
480    'maskutil_W'  : maskutil_W ,
481    'maskutil_F'  : maskutil_F ,
482    'bounds_lon_grid_T': bounds_lon_grid_T,
483    'bounds_lat_grid_T': bounds_lat_grid_T,
484    'bounds_lon_grid_U': bounds_lon_grid_U,
485    'bounds_lat_grid_U': bounds_lat_grid_U,
486    'bounds_lon_grid_V': bounds_lon_grid_V,
487    'bounds_lat_grid_V': bounds_lat_grid_V,
488    'bounds_lon_grid_W': bounds_lon_grid_W,
489    'bounds_lat_grid_W': bounds_lat_grid_W,
490    'bounds_lon_grid_F': bounds_lon_grid_F,
491    'bounds_lat_grid_F': bounds_lat_grid_F,
492})
493
494#replace nav_lon nav_lat with variables obtained from NEMO coordinates.nc file
495#by construction nav_lon nav_lat are those of the bathymetry
496for cd_type in ['T', 'U', 'V', 'F', 'W'] :
497    for dir_type in ['lat', 'lon'] :
498        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
499        ds.coords[coord_name]=locals()[coord_name]
500
501ds.attrs['name']         = 'coordinates_mask'
502ds.attrs['description']  = 'coordinates and mask for MOSAIX'
503ds.attrs['title']        = 'coordinates_mask'
504ds.attrs['source']       = 'IPSL Earth system model'
505ds.attrs['group']        = 'ICMC IPSL Climate Modelling Center'
506ds.attrs['Institution']  = 'IPSL https.//www.ipsl.fr'
507ds.attrs['Model']        = model
508ds.attrs['timeStamp']    = '{:%Y-%b-%d %H:%M:%S}'.format (datetime.datetime.now ())
509ds.attrs['history']      = 'Build from ' + n_coord + ' and ' + n_bathy
510ds.attrs['directory']    = os.getcwd     ()
511ds.attrs['user']         = os.getlogin   ()
512ds.attrs['HOSTNAME']     = platform.node ()
513ds.attrs['Python']       = 'Python version: ' +  platform.python_version ()
514ds.attrs['xarray']       = 'xarray version: ' +  xr.__version__
515ds.attrs['OS']           = platform.system  ()
516ds.attrs['release']      = platform.release ()
517ds.attrs['hardware']     = platform.machine ()
518
519ds.attrs['SVN_Author']   = "$Author:  $"
520ds.attrs['SVN_Date']     = "$Date:  $"
521ds.attrs['SVN_Revision'] = "$Revision:  $"
522ds.attrs['SVN_Id']       = "$Id: Build_coordinates_mask.py $"
523ds.attrs['SVN_HeadURL']  = "$HeadURL: http://forge.ipsl.jussieu.fr/igcmg/svn/TOOLS/MOSAIX/Build_coordinates_mask.py $"
524
525# save to output file
526ds.to_netcdf (n_out)
Note: See TracBrowser for help on using the repository browser.