Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

Build_coordinates_mask.ipynb

Last change on this file was 5941, checked in by snguyen, 3 years ago
post final version ;-) which applies zero to southernmost boundary only if nperio = (3, 4, 5, 6)
Property svn:keywords set to `Date Revision HeadURL Author Id`
File size: 32.5 KB

Line
1	{
2	"cells": [
3	{
4	"cell_type": "markdown",
5	"id": "10fcba33",
6	"metadata": {},
7	"source": [
8	"# Constructs ```coordinates_mask``` file\n",
9	"The ```coordinates_mask``` file is needed to run MOSAIX, to build interpolation weights for the coupled model"
10	]
11	},
12	{
13	"cell_type": "markdown",
14	"id": "b1cbb871",
15	"metadata": {},
16	"source": [
17	"Olivier Marti - olivier.marti@lsce.ipsl - 10/06/2021"
18	]
19	},
20	{
21	"cell_type": "code",
22	"execution_count": 1,
23	"id": "731d770c",
24	"metadata": {},
25	"outputs": [],
26	"source": [
27	"import xarray as xr\n",
28	"import numpy as np\n",
29	"import nemo\n",
30	"import datetime, os, platform"
31	]
32	},
33	{
34	"cell_type": "markdown",
35	"id": "fb78f0e5",
36	"metadata": {},
37	"source": [
38	"# Various Switches To Deal With Bugs and Implementation Choices"
39	]
40	},
41	{
42	"cell_type": "code",
43	"execution_count": 2,
44	"id": "2bd430d4",
45	"metadata": {},
46	"outputs": [],
47	"source": [
48	"# reproduce periodicity bug for U and nperio = 4\n",
49	"nemoUbug = False\n",
50	"# change metrics (and bathymetry) for straits\n",
51	"straits = False\n",
52	"# periodicity imposed on coordinates, metrics and areas\n",
53	"coordperio = False\n",
54	"# function used to compute mask_T from bathymetry\n",
55	"maskbathynemo = False"
56	]
57	},
58	{
59	"cell_type": "markdown",
60	"id": "26422e3f",
61	"metadata": {},
62	"source": [
63	"# Input files, type of ORCA grid, and output file"
64	]
65	},
66	{
67	"cell_type": "code",
68	"execution_count": 3,
69	"id": "deb66d3d",
70	"metadata": {},
71	"outputs": [],
72	"source": [
73	"#model = 'eORCA1.2'\n",
74	"#model = 'ORCA2.3'\n",
75	"model = 'paleORCA'\n",
76	"#model = 'ORCA2.4'\n",
77	"#model = 'ORCA.05'http://localhost:8888/notebooks/Build_coordinates_mask_new_2409.ipynb#"
78	]
79	},
80	{
81	"cell_type": "code",
82	"execution_count": 4,
83	"id": "2475c8b7",
84	"metadata": {},
85	"outputs": [],
86	"source": [
87	"if model == 'eORCA1.2' :\n",
88	" n_coord = 'eORCA1.2coordinates.nc' # 'eORCA_R1_coordinates_v1.0.nc' \n",
89	" n_bathy = 'eORCA1.2bathy_meter.nc' # 'eORCA_R1_bathy_meter_v2.2.nc'\n",
90	" #n_coord = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_coordinates_v1.0.nc'\n",
91	" #n_bathy = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_bathy_meter_v2.2.nc'\n",
92	" nperio = 6\n",
93	" n_out = 'My_eORCA_R1_coordinates_mask_test.nc'\n",
94	"\n",
95	"if model == 'ORCA2.3' :\n",
96	" nperio = 4\n",
97	" n_coord = 'coordinates.nc' # 'ORCA2.3_coordinates.nc' # 'coordinates_ORCA2_LIM3_PISCES.nc' #\n",
98	" n_bathy = 'bathy_meter_closea_topo.nc' # 'ORCA2.3_bathy_meter.nc' # 'bathy_meter_ORCA2_LIM3_PISCES.nc' #\n",
99	" n_out = 'My_test.nc'\n",
100	" \n",
101	"if model == 'ORCA2.4' :\n",
102	" nperio = 4\n",
103	" n_coord = 'ORCA2.4_coordinates.nc'\n",
104	" n_bathy = 'ORCA2.4_bathy_meter.nc'\n",
105	" n_out = 'ORCA2.4_coordinates_mask_test.nc'\n",
106	" \n",
107	"if model == 'paleORCA' :\n",
108	" n_coord = 'coordinates_paleo.nc'\n",
109	" n_bathy = 'bathy_meter_paleo.nc'\n",
110	" nperio = 6\n",
111	" n_out = 'My_coordinates_mask_paleo.nc'\n",
112	"\n"
113	]
114	},
115	{
116	"cell_type": "markdown",
117	"id": "8cd952a6",
118	"metadata": {},
119	"source": [
120	"## Open input files\n",
121	"Do not decode time which might be buggy in some files\n",
122	"Suppress singleton dimensions if necessary"
123	]
124	},
125	{
126	"cell_type": "code",
127	"execution_count": 5,
128	"id": "0b7c7279",
129	"metadata": {
130	"scrolled": false
131	},
132	"outputs": [],
133	"source": [
134	"f_coord = xr.open_dataset (n_coord, decode_times=False).squeeze()\n",
135	"f_bathy = xr.open_dataset (n_bathy, decode_times=False).squeeze()"
136	]
137	},
138	{
139	"cell_type": "code",
140	"execution_count": 6,
141	"id": "9fb92919",
142	"metadata": {},
143	"outputs": [
144	{
145	"name": "stdout",
146	"output_type": "stream",
147	"text": [
148	"failed to suppress time\n"
149	]
150	}
151	],
152	"source": [
153	"# Suppress time if necessary\n",
154	"try :\n",
155	" del f_coord['time']\n",
156	" print ('success')\n",
157	"except :\n",
158	" pass\n",
159	" print ('failed to suppress time')"
160	]
161	},
162	{
163	"cell_type": "markdown",
164	"id": "d61fb6b2",
165	"metadata": {},
166	"source": [
167	"## Creates masks"
168	]
169	},
170	{
171	"cell_type": "markdown",
172	"id": "299c0a66",
173	"metadata": {},
174	"source": [
175	"### Read bathymetry "
176	]
177	},
178	{
179	"cell_type": "code",
180	"execution_count": 7,
181	"id": "9dfd90e6",
182	"metadata": {
183	"scrolled": true
184	},
185	"outputs": [
186	{
187	"name": "stdout",
188	"output_type": "stream",
189	"text": [
190	"Normal case\n"
191	]
192	}
193	],
194	"source": [
195	"Bathymetry = f_bathy['Bathymetry'].copy()\n",
196	"\n",
197	"try :\n",
198	" Bathymetry = Bathymetry.rename ({'nav_lat':'nav_lat_grid_T', 'nav_lon':'nav_lon_grid_T'})\n",
199	" print ('Normal case')\n",
200	"except : \n",
201	" print ('Exception')\n",
202	" nav_lon_grid_T = f_bathy ['nav_lon']\n",
203	" nav_lat_grid_T = f_bathy ['nav_lat']\n",
204	" Bathymetry = xr.DataArray (Bathymetry, coords = { \"nav_lat_grid_T\": ([\"y\", \"x\"], nav_lat_grid_T),\n",
205	" \"nav_lon_grid_T\": ([\"y\", \"x\"], nav_lon_grid_T) } )\n",
206	" \n",
207	"Bathymetry = Bathymetry.rename ({'y':'y_grid_T', 'x':'x_grid_T'})"
208	]
209	},
210	{
211	"cell_type": "code",
212	"execution_count": 8,
213	"id": "2606d14b",
214	"metadata": {},
215	"outputs": [],
216	"source": [
217	"if straits :\n",
218	" # Open straits for ORCA2 grid\n",
219	" if model == 'ORCA2.3' :\n",
220	" # orca_r2: Gibraltar strait open\n",
221	" Bathymetry[101,139] = 284.\n",
222	" # orca_r2: Bab el Mandeb strait open\n",
223	" Bathymetry[87,159] = 137."
224	]
225	},
226	{
227	"cell_type": "markdown",
228	"id": "71c9c821",
229	"metadata": {},
230	"source": [
231	"### Determine which periodicity is needed"
232	]
233	},
234	{
235	"cell_type": "code",
236	"execution_count": 9,
237	"id": "7e7a49d5",
238	"metadata": {},
239	"outputs": [
240	{
241	"name": "stdout",
242	"output_type": "stream",
243	"text": [
244	"nperio specified : 6\n"
245	]
246	}
247	],
248	"source": [
249	"jpj, jpi = Bathymetry.shape\n",
250	"try : \n",
251	" if nperio != None :\n",
252	" print ('nperio specified : ', nperio)\n",
253	"except :\n",
254	" nperio = None\n",
255	" if jpi == 182 : nperio = 4 # ORCA2. \\!/ We choose legacy orca2\n",
256	" if jpi == 362 : nperio = 6 # ORCA1.\n",
257	" if jpi == 1442 : nperio = 6 # ORCA025.\n",
258	" #\n",
259	" if nperio == None :\n",
260	" sys.exit ('nperio not found, and cannot by guessed' )\n",
261	" else :\n",
262	" print ('nperio set as {:d} (deduced from data size {:d} {:d})'.format (nperio, jpj, jpi))"
263	]
264	},
265	{
266	"cell_type": "code",
267	"execution_count": 10,
268	"id": "34cdcfe4",
269	"metadata": {},
270	"outputs": [],
271	"source": [
272	"Bathymetry = nemo.lbc (Bathymetry, nperio=nperio, cd_type='T')\n",
273	"#Suppress ocean points at southernmost position for periodicity (3, 4, 5, 6)\n",
274	"if nperio in (3 , 4 , 5 , 6) :\n",
275	" Bathymetry[0,:] = 0.0 "
276	]
277	},
278	{
279	"cell_type": "markdown",
280	"id": "0da91538",
281	"metadata": {},
282	"source": [
283	"### Creates ```mask_T```"
284	]
285	},
286	{
287	"cell_type": "code",
288	"execution_count": 11,
289	"id": "49a654ef",
290	"metadata": {},
291	"outputs": [],
292	"source": [
293	"if maskbathynemo :\n",
294	" # Use same formula as domzgr. Should be equivalent to Bathimetry > 0.0. \n",
295	" mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='f4')\n",
296	"else :\n",
297	" mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='f4')"
298	]
299	},
300	{
301	"cell_type": "markdown",
302	"id": "327d48ea",
303	"metadata": {},
304	"source": [
305	"### Creates other masks"
306	]
307	},
308	{
309	"cell_type": "code",
310	"execution_count": 12,
311	"id": "f0b656a1",
312	"metadata": {},
313	"outputs": [],
314	"source": [
315	"mask_U = mask_T * mask_T.shift (x_grid_T=-1)\n",
316	"mask_V = mask_T * mask_T.shift (y_grid_T=-1)\n",
317	"mask_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)\n",
318	"mask_W = mask_T\n",
319	"\n",
320	"mask_U = nemo.lbc (mask_U, nperio=nperio, cd_type='U', nemo_4U_bug=nemoUbug).astype (dtype='f4')\n",
321	"mask_V = nemo.lbc (mask_V, nperio=nperio, cd_type='V').astype (dtype='f4')\n",
322	"mask_F = nemo.lbc (mask_F, nperio=nperio, cd_type='F').astype (dtype='f4')\n",
323	"mask_W = nemo.lbc (mask_W, nperio=nperio, cd_type='W').astype (dtype='f4')\n",
324	"\n",
325	"mask_U = mask_U.rename ( {'y_grid_T' : 'y_grid_U' , 'x_grid_T' : 'x_grid_U', \n",
326	" 'nav_lat_grid_T': 'nav_lat_grid_U', 'nav_lon_grid_T': 'nav_lon_grid_U'} )\n",
327	"mask_V = mask_V.rename ( {'y_grid_T' : 'y_grid_V' , 'x_grid_T' : 'x_grid_V',\n",
328	" 'nav_lat_grid_T': 'nav_lat_grid_V', 'nav_lon_grid_T': 'nav_lon_grid_V'} )\n",
329	"mask_W = mask_W.rename ( {'y_grid_T' : 'y_grid_W' , 'x_grid_T' : 'x_grid_W',\n",
330	" 'nav_lat_grid_T': 'nav_lat_grid_W', 'nav_lon_grid_T': 'nav_lon_grid_W'} )\n",
331	"mask_F = mask_F.rename ( {'y_grid_T' : 'y_grid_F' , 'x_grid_T' : 'x_grid_F',\n",
332	" 'nav_lat_grid_T': 'nav_lat_grid_F', 'nav_lon_grid_T': 'nav_lon_grid_F'} )\n",
333	"\n",
334	"mask_T.name = 'mask_T'\n",
335	"mask_U.name = 'mask_U'\n",
336	"mask_V.name = 'mask_V'\n",
337	"mask_F.name = 'mask_F'\n",
338	"mask_W.name = 'mask_W'"
339	]
340	},
341	{
342	"cell_type": "markdown",
343	"id": "5cb36da6",
344	"metadata": {},
345	"source": [
346	"### Masks with duplicate points removed"
347	]
348	},
349	{
350	"cell_type": "code",
351	"execution_count": 13,
352	"id": "274e4321",
353	"metadata": {},
354	"outputs": [],
355	"source": [
356	"for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
357	" MaskName = 'mask_' + cd_type\n",
358	" UtilName = 'maskutil_' + cd_type\n",
359	" locals()[UtilName] = nemo.lbc_mask (locals()[MaskName].copy(), nperio=nperio, cd_type=cd_type).astype (dtype='f4')\n",
360	" locals()[UtilName].name = UtilName\n",
361	" locals()[UtilName].encoding['_FillValue'] = None"
362	]
363	},
364	{
365	"cell_type": "markdown",
366	"id": "69ec3daa",
367	"metadata": {},
368	"source": [
369	"### Add attributes"
370	]
371	},
372	{
373	"cell_type": "code",
374	"execution_count": 14,
375	"id": "f46e1c93",
376	"metadata": {},
377	"outputs": [],
378	"source": [
379	"mask_T.attrs ['cell_measures'] = 'area: area_grid_T'\n",
380	"mask_U.attrs ['cell_measures'] = 'area: area_grid_U'\n",
381	"mask_V.attrs ['cell_measures'] = 'area: area_grid_V'\n",
382	"mask_W.attrs ['cell_measures'] = 'area: area_grid_W'\n",
383	"mask_F.attrs ['cell_measures'] = 'area: area_grid_F'\n",
384	"\n",
385	"maskutil_T.attrs['cell_measures'] = 'area: area_grid_T'\n",
386	"maskutil_U.attrs['cell_measures'] = 'area: area_grid_U'\n",
387	"maskutil_V.attrs['cell_measures'] = 'area: area_grid_V'\n",
388	"maskutil_W.attrs['cell_measures'] = 'area: area_grid_W'\n",
389	"maskutil_F.attrs['cell_measures'] = 'area: area_grid_F'"
390	]
391	},
392	{
393	"cell_type": "markdown",
394	"id": "580e2776",
395	"metadata": {},
396	"source": [
397	"## Creates grid coordinates and surfaces"
398	]
399	},
400	{
401	"cell_type": "code",
402	"execution_count": 15,
403	"id": "a7ce1490",
404	"metadata": {},
405	"outputs": [],
406	"source": [
407	"if coordperio :\n",
408	" nav_lon_grid_T = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='T')\n",
409	" nav_lat_grid_T = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='T')\n",
410	" nav_lon_grid_U = nemo.lbc (f_coord ['glamu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)\n",
411	" nav_lat_grid_U = nemo.lbc (f_coord ['gphiu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)\n",
412	" nav_lon_grid_V = nemo.lbc (f_coord ['glamv'].copy(), nperio=nperio, cd_type='V')\n",
413	" nav_lat_grid_V = nemo.lbc (f_coord ['gphiv'].copy(), nperio=nperio, cd_type='V')\n",
414	" nav_lon_grid_W = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='W')\n",
415	" nav_lat_grid_W = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='W')\n",
416	" nav_lon_grid_F = nemo.lbc (f_coord ['glamf'].copy(), nperio=nperio, cd_type='F')\n",
417	" nav_lat_grid_F = nemo.lbc (f_coord ['gphif'].copy(), nperio=nperio, cd_type='F')\n",
418	"else :\n",
419	" nav_lon_grid_T = f_coord ['glamt'].copy()\n",
420	" nav_lat_grid_T = f_coord ['gphit'].copy()\n",
421	" nav_lon_grid_U = f_coord ['glamu'].copy()\n",
422	" nav_lat_grid_U = f_coord ['gphiu'].copy()\n",
423	" nav_lon_grid_V = f_coord ['glamv'].copy()\n",
424	" nav_lat_grid_V = f_coord ['gphiv'].copy()\n",
425	" nav_lon_grid_W = f_coord ['glamt'].copy()\n",
426	" nav_lat_grid_W = f_coord ['gphit'].copy()\n",
427	" nav_lon_grid_F = f_coord ['glamf'].copy()\n",
428	" nav_lat_grid_F = f_coord ['gphif'].copy()\n",
429	"\n",
430	"nav_lon_grid_T = nav_lon_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )\n",
431	"nav_lat_grid_T = nav_lat_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )\n",
432	"nav_lon_grid_U = nav_lon_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )\n",
433	"nav_lat_grid_U = nav_lat_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )\n",
434	"nav_lon_grid_V = nav_lon_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )\n",
435	"nav_lat_grid_V = nav_lat_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )\n",
436	"nav_lon_grid_W = nav_lon_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )\n",
437	"nav_lat_grid_W = nav_lat_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )\n",
438	"nav_lon_grid_F = nav_lon_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )\n",
439	"nav_lat_grid_F = nav_lat_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )\n",
440	"\n",
441	"# remove _FillValue and missing_value\n",
442	"for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
443	" for dir_type in ['lat', 'lon'] :\n",
444	" coord_name = 'nav_' + dir_type + '_grid_' + cd_type\n",
445	" locals()[coord_name].encoding['_FillValue'] = None\n",
446	" locals()[coord_name].encoding['missing_value'] = None"
447	]
448	},
449	{
450	"cell_type": "code",
451	"execution_count": 16,
452	"id": "516b4c5c",
453	"metadata": {},
454	"outputs": [],
455	"source": [
456	"nav_lon_grid_T.name = 'nav_lon_grid_T'\n",
457	"nav_lat_grid_T.name = 'nav_lat_grid_T'\n",
458	"nav_lon_grid_U.name = 'nav_lon_grid_U'\n",
459	"nav_lat_grid_U.name = 'nav_lat_grid_U'\n",
460	"nav_lon_grid_V.name = 'nav_lon_grid_V'\n",
461	"nav_lat_grid_V.name = 'nav_lat_grid_V'\n",
462	"nav_lon_grid_F.name = 'nav_lon_grid_F'\n",
463	"nav_lat_grid_F.name = 'nav_lat_grid_F'\n",
464	"nav_lon_grid_W.name = 'nav_lon_grid_W'\n",
465	"nav_lat_grid_W.name = 'nav_lat_grid_W'"
466	]
467	},
468	{
469	"cell_type": "markdown",
470	"id": "cffad9b4",
471	"metadata": {},
472	"source": [
473	"### Add attributes"
474	]
475	},
476	{
477	"cell_type": "code",
478	"execution_count": 17,
479	"id": "de0e6514",
480	"metadata": {},
481	"outputs": [],
482	"source": [
483	"nav_lon_grid_T.attrs['standard_name'] = 'longitude'\n",
484	"nav_lon_grid_T.attrs['long_name'] = 'Longitude'\n",
485	"nav_lon_grid_T.attrs['units'] = 'degrees_east'\n",
486	"nav_lat_grid_T.attrs['standard_name'] = 'latitude'\n",
487	"nav_lat_grid_T.attrs['long_name'] = 'Latitude'\n",
488	"nav_lat_grid_T.attrs['units'] = 'degrees_north'\n",
489	"\n",
490	"nav_lon_grid_U.attrs['standard_name'] = 'longitude'\n",
491	"nav_lon_grid_U.attrs['long_name'] = 'Longitude'\n",
492	"nav_lon_grid_U.attrs['units'] = 'degrees_east'\n",
493	"nav_lat_grid_U.attrs['standard_name'] = 'latitude'\n",
494	"nav_lat_grid_U.attrs['long_name'] = 'Latitude'\n",
495	"nav_lat_grid_U.attrs['units'] = 'degrees_north'\n",
496	"\n",
497	"nav_lon_grid_V.attrs['standard_name'] = 'longitude'\n",
498	"nav_lon_grid_V.attrs['long_name'] = 'Longitude'\n",
499	"nav_lon_grid_V.attrs['units'] = 'degrees_east'\n",
500	"nav_lat_grid_V.attrs['standard_name'] = 'latitude'\n",
501	"nav_lat_grid_V.attrs['long_name'] = 'Latitude'\n",
502	"nav_lat_grid_V.attrs['units'] = 'degrees_north'\n",
503	"\n",
504	"nav_lon_grid_W.attrs['standard_name'] = 'longitude'\n",
505	"nav_lon_grid_W.attrs['long_name'] = 'Longitude'\n",
506	"nav_lon_grid_W.attrs['units'] = 'degrees_east'\n",
507	"nav_lat_grid_W.attrs['standard_name'] = 'latitude'\n",
508	"nav_lat_grid_W.attrs['long_name'] = 'Latitude'\n",
509	"nav_lat_grid_W.attrs['units'] = 'degrees_north'\n",
510	"\n",
511	"nav_lon_grid_F.attrs['standard_name'] = 'longitude'\n",
512	"nav_lon_grid_F.attrs['long_name'] = 'Longitude'\n",
513	"nav_lon_grid_F.attrs['units'] = 'degrees_east'\n",
514	"nav_lat_grid_F.attrs['standard_name'] = 'latitude'\n",
515	"nav_lat_grid_F.attrs['long_name'] = 'Latitude'\n",
516	"nav_lat_grid_F.attrs['units'] = 'degrees_north'\n",
517	"\n",
518	"nav_lon_grid_T.attrs['bounds'] = 'bounds_lon_grid_T'\n",
519	"nav_lat_grid_T.attrs['bounds'] = 'bounds_lat_grid_T'\n",
520	"nav_lon_grid_U.attrs['bounds'] = 'bounds_lon_grid_U'\n",
521	"nav_lat_grid_U.attrs['bounds'] = 'bounds_lat_grid_U'\n",
522	"nav_lon_grid_V.attrs['bounds'] = 'bounds_lon_grid_V'\n",
523	"nav_lat_grid_V.attrs['bounds'] = 'bounds_lat_grid_V'\n",
524	"nav_lon_grid_W.attrs['bounds'] = 'bounds_lon_grid_W'\n",
525	"nav_lat_grid_W.attrs['bounds'] = 'bounds_lat_grid_W'\n",
526	"nav_lon_grid_F.attrs['bounds'] = 'bounds_lon_grid_F'\n",
527	"nav_lat_grid_F.attrs['bounds'] = 'bounds_lat_grid_F'"
528	]
529	},
530	{
531	"cell_type": "code",
532	"execution_count": 18,
533	"id": "8e4d997c",
534	"metadata": {},
535	"outputs": [],
536	"source": [
537	"# create new variables e1 e2 to keep f_coord the same\n",
538	"for cd_type in ['t', 'u', 'v', 'f'] :\n",
539	" for axis in ['1', '2'] :\n",
540	" coordName = 'e' + axis + cd_type\n",
541	" locals()[coordName]=f_coord[coordName].copy()\n",
542	"# remove zero values from areas\n",
543	"# need to be define for the extended grid south of -80S\n",
544	"# some point are undefined but you need to have e1 and e2 .NE. 0\n",
545	" locals()[coordName]=xr.where(locals()[coordName] == 0.0, 1.0e2, locals()[coordName])"
546	]
547	},
548	{
549	"cell_type": "code",
550	"execution_count": 19,
551	"id": "175bae75",
552	"metadata": {},
553	"outputs": [],
554	"source": [
555	"# Correct areas for straits\n",
556	"if straits :\n",
557	" # ORCA R2 configuration\n",
558	" if model == 'ORCA2.3' :\n",
559	" # Gibraltar : e2u reduced to 20 km\n",
560	" e2u[101,138:140] = 20.e3\n",
561	" # Bab el Mandeb : e2u reduced to 30 km\n",
562	" # e1v reduced to 18 km\n",
563	" e1v[87,159] = 18.e3\n",
564	" e2u[87,159] = 30.e3\n",
565	" # Danish Straits: e2u reduced to 10 km\n",
566	" e2u[115,144:146] = 10.e3\n",
567	" # ORCA R1 configuration\n",
568	" if model == 'eORCA1.2' :\n",
569	" # Gibraltar : e2u reduced to 20 km\n",
570	" e2u[240,281:283] = 20.e3\n",
571	" # Bhosporus : e2u reduced to 10 km\n",
572	" e2u[247,313:315] = 10.e3\n",
573	" # Lombok : e1v reduced to 13 km\n",
574	" e1v[163:165,43] = 13.e3\n",
575	" # Sumba : e1v reduced to 8 km\n",
576	" e1v[163:165,47] = 8.e3\n",
577	" # Ombai : e1v reduced to 13 km\n",
578	" e1v[163:165,52] = 13.e3\n",
579	" # Timor Passage : e1v reduced to 20 km\n",
580	" #e1v[163:165,55] = 20.e3\n",
581	" # W Halmahera : e1v reduced to 30 km\n",
582	" e1v[180:182,54] = 30.e3\n",
583	" # E Halmahera : e1v reduced to 50 km\n",
584	" e1v[180:182,57] = 50.e3\n",
585	" # ORCA R05 configuration\n",
586	" if model == 'ORCA.05' :\n",
587	" # Reduced e2u at the Gibraltar Strait\n",
588	" e2u[326,562:564] = 20.e3\n",
589	" # Reduced e2u at the Bosphore Strait\n",
590	" e2u[342,626:628] = 10.e3\n",
591	" # Reduced e2u at the Sumba Strait\n",
592	" e2u[231,92:94] = 40.e3\n",
593	" # Reduced e2u at the Ombai Strait\n",
594	" e2u[231,102] = 15.e3\n",
595	" # Reduced e2u at the Palk Strait\n",
596	" e2u[269,14] = 10.e3\n",
597	" # Reduced e1v at the Lombok Strait\n",
598	" e1v[231:233,86] = 10.e3\n",
599	" # Reduced e1v at the Bab el Mandeb\n",
600	" e1v[275,661] = 25.e3"
601	]
602	},
603	{
604	"cell_type": "code",
605	"execution_count": 20,
606	"id": "4f6cfbca",
607	"metadata": {},
608	"outputs": [],
609	"source": [
610	"if coordperio :\n",
611	" area_grid_T = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='T')\n",
612	" area_grid_U = nemo.lbc (e1u*e2u, nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)\n",
613	" area_grid_V = nemo.lbc (e1v*e2v, nperio=nperio, cd_type='V')\n",
614	" area_grid_W = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='W')\n",
615	" area_grid_F = nemo.lbc (e1f*e2f, nperio=nperio, cd_type='F')\n",
616	"else :\n",
617	" area_grid_T = e1t*e2t\n",
618	" area_grid_U = e1u*e2u\n",
619	" area_grid_V = e1v*e2v\n",
620	" area_grid_W = e1t*e2t\n",
621	" area_grid_F = e1f*e2f\n",
622	"\n",
623	"\n",
624	"area_grid_T.name = 'area_grid_T'\n",
625	"area_grid_U.name = 'area_grid_U'\n",
626	"area_grid_V.name = 'area_grid_V'\n",
627	"area_grid_F.name = 'area_grid_F'\n",
628	"area_grid_W.name = 'area_grid_W'\n",
629	"\n",
630	"area_grid_T = area_grid_T.rename ({'y':'y_grid_T', 'x':'x_grid_T'})\n",
631	"area_grid_U = area_grid_U.rename ({'y':'y_grid_U', 'x':'x_grid_U'})\n",
632	"area_grid_V = area_grid_V.rename ({'y':'y_grid_V', 'x':'x_grid_V'})\n",
633	"area_grid_W = area_grid_W.rename ({'y':'y_grid_W', 'x':'x_grid_W'})\n",
634	"area_grid_F = area_grid_F.rename ({'y':'y_grid_F', 'x':'x_grid_F'})\n",
635	"\n",
636	"area_grid_T.attrs['standard_name'] = 'cell_area'\n",
637	"area_grid_T.attrs['units'] = 'm2'\n",
638	"#area_grid_T.attrs['coordinates'] = 'nav_lat_grid_T nav_lon_grid_T'\n",
639	"area_grid_U.attrs['standard_name'] = 'cell_area'\n",
640	"area_grid_U.attrs['units'] = 'm2'\n",
641	"#area_grid_U.attrs['coordinates'] = 'nav_lat_grid_U nav_lon_grid_U'\n",
642	"area_grid_V.attrs['standard_name'] = 'cell_area'\n",
643	"area_grid_V.attrs['units'] = 'm2'\n",
644	"#area_grid_V.attrs['coordinates'] = 'nav_lat_grid_V nav_lon_grid_V'\n",
645	"area_grid_W.attrs['standard_name'] = 'cell_area'\n",
646	"area_grid_W.attrs['units'] = 'm2'\n",
647	"#area_grid_W.attrs['coordinates'] = 'nav_lat_grid_W nav_lon_grid_W'\n",
648	"area_grid_F.attrs['standard_name'] = 'cell_area'\n",
649	"area_grid_F.attrs['units'] = 'm2'\n",
650	"#area_grid_F.attrs['coordinates'] = 'nav_lat_grid_F nav_lon_grid_F'\n",
651	"\n",
652	"for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
653	" areaName = 'area_grid_' + cd_type\n",
654	" locals()[areaName].encoding['_FillValue'] = None"
655	]
656	},
657	{
658	"cell_type": "markdown",
659	"id": "3c6f7b08",
660	"metadata": {},
661	"source": [
662	"## Construct bounds"
663	]
664	},
665	{
666	"cell_type": "code",
667	"execution_count": 21,
668	"id": "69c404d9",
669	"metadata": {},
670	"outputs": [],
671	"source": [
672	"def set_bounds (cdgrd) :\n",
673	" '''\n",
674	" Constructs lon/lat bounds\n",
675	" Bounds are numerated counter clockwise, from bottom left\n",
676	" See NEMO file OPA_SRC/IOM/iom.F90, ROUTINE set_grid_bounds, for more details\n",
677	" '''\n",
678	" # Define offset of coordinate representing bottom-left corner\n",
679	" if cdgrd in ['T', 'W'] : \n",
680	" icnr = -1 ; jcnr = -1\n",
681	" corner_lon = f_coord ['glamf'].copy() ; corner_lat = f_coord ['gphif'].copy()\n",
682	" center_lon = f_coord ['glamt'].copy() ; center_lat = f_coord ['gphit'].copy()\n",
683	" if cdgrd == 'U' : \n",
684	" icnr = 0 ; jcnr = -1\n",
685	" corner_lon = f_coord ['glamv'].copy() ; corner_lat = f_coord ['gphiv'].copy()\n",
686	" center_lon = f_coord ['glamu'].copy() ; center_lat = f_coord ['gphiu'].copy()\n",
687	" if cdgrd == 'V' : \n",
688	" icnr = -1 ; jcnr = 0\n",
689	" corner_lon = f_coord ['glamu'].copy() ; corner_lat = f_coord ['gphiu'].copy()\n",
690	" center_lon = f_coord ['glamv'].copy() ; center_lat = f_coord ['gphiv'].copy()\n",
691	" if cdgrd == 'F' : \n",
692	" icnr = -1 ; jcnr = -1\n",
693	" corner_lon = f_coord ['glamt'].copy() ; corner_lat = f_coord ['gphit'].copy()\n",
694	" center_lon = f_coord ['glamf'].copy() ; center_lat = f_coord ['gphif'].copy()\n",
695	" \n",
696	" #y_grid, x_grid = corner_lon.shape ;\n",
697	" nvertex = 4\n",
698	" dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]\n",
699	" \n",
700	" bounds_lon = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)\n",
701	" bounds_lat = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)\n",
702	" \n",
703	" idx = [(jcnr,icnr), (jcnr,icnr+1), (jcnr+1,icnr+1), (jcnr+1,icnr)]\n",
704	" \n",
705	" # Compute cell vertices that can be defined, \n",
706	" # and complete with periodicity\n",
707	" for nn in range (nvertex) :\n",
708	" tmp = np.roll (corner_lon, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))\n",
709	" bounds_lon[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]\n",
710	" tmp = np.roll (corner_lat, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))\n",
711	" bounds_lat[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]\n",
712	" bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)\n",
713	" bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)\n",
714	" \n",
715	" # Zero-size cells at closed boundaries if cell points provided,\n",
716	" # otherwise they are closed cells with unrealistic bounds\n",
717	" if not (nperio == 1 or nperio == 4 or nperio == 6) :\n",
718	" for nn in range (nvertex) : \n",
719	" bounds_lon[:,0,nn] = center_lon[:,0] # (West or jpni = 1), closed E-W\n",
720	" bounds_lat[:,0,nn] = center_lat[:,0]\n",
721	" bounds_lon[:,jpi,nn] = center_lon[:,jpi] # (East or jpni = 1), closed E-W\n",
722	" bounds_lat[:,jpi,nn] = center_lat[:,jpi]\n",
723	" if nperio != 2 :\n",
724	" for nn in range (nvertex) :\n",
725	" bounds_lon[0,:,nn] = center_lon[0,:] # (South or jpnj = 1), not symmetric\n",
726	" bounds_lat[0,:,nn] = center_lat[0,:]\n",
727	" if nperio < 3 :\n",
728	" for nn in range (nvertex) :\n",
729	" bounds_lon[jpj,:,nn] = center_lon[jpj,:] # (North or jpnj = 1), no north fold\n",
730	" bounds_lat[jpj,:,nn] = center_lat[jpj,:]\n",
731	" \n",
732	" # Rotate cells at the north fold\n",
733	" if nperio >= 3 :\n",
734	" # Working array for location of northfold\n",
735	" z_fld = nemo.lbc (np.ones ((jpj, jpi)), nperio=nperio, cd_type=cdgrd, psgn=-1., nemo_4U_bug=nemoUbug)\n",
736	" z_fld = np.repeat((z_fld == -1.0)[...,np.newaxis],4,axis=2)\n",
737	" # circular shift of 2 indices in bounds third index\n",
738	" bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)\n",
739	" bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)\n",
740	" bounds_lon[:,:,:] = np.where (z_fld, bounds_lon_tmp[:,:,:] , bounds_lon[:,:,:] )\n",
741	" bounds_lat[:,:,:] = np.where (z_fld, bounds_lat_tmp[:,:,:] , bounds_lat[:,:,:] )\n",
742	" \n",
743	" # Invert cells at the symmetric equator\n",
744	" if nperio == 2 :\n",
745	" bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)\n",
746	" bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)\n",
747	" bounds_lon [0,:,:] = bounds_lon [0,:,:]\n",
748	" bounds_lat [0,:,:] = bounds_lat [0,:,:]\n",
749	" \n",
750	" #bounds_lon.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd\n",
751	" #bounds_lat.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd\n",
752	" bounds_lon.attrs['units'] = 'degrees_east'\n",
753	" bounds_lat.attrs['units'] = 'degrees_north'\n",
754	" bounds_lon.name = 'bounds_lon_grid_' + cdgrd\n",
755	" bounds_lat.name = 'bounds_lat_grid_' + cdgrd\n",
756	" # remove _FillValue\n",
757	" bounds_lon.encoding['_FillValue'] = None\n",
758	" bounds_lat.encoding['_FillValue'] = None\n",
759	"\n",
760	" return bounds_lon, bounds_lat"
761	]
762	},
763	{
764	"cell_type": "code",
765	"execution_count": 22,
766	"id": "87a8f592",
767	"metadata": {},
768	"outputs": [],
769	"source": [
770	"bounds_lon_grid_T, bounds_lat_grid_T = set_bounds ('T')\n",
771	"bounds_lon_grid_U, bounds_lat_grid_U = set_bounds ('U')\n",
772	"bounds_lon_grid_V, bounds_lat_grid_V = set_bounds ('V')\n",
773	"bounds_lon_grid_W, bounds_lat_grid_W = set_bounds ('W')\n",
774	"bounds_lon_grid_F, bounds_lat_grid_F = set_bounds ('F')"
775	]
776	},
777	{
778	"cell_type": "markdown",
779	"id": "5dbf22b3",
780	"metadata": {},
781	"source": [
782	"## Save Data in a NetCDF file"
783	]
784	},
785	{
786	"cell_type": "code",
787	"execution_count": 23,
788	"id": "ca87fa1b",
789	"metadata": {},
790	"outputs": [],
791	"source": [
792	"ds = xr.Dataset ({\n",
793	" 'mask_T' : mask_T ,\n",
794	" 'mask_U' : mask_U ,\n",
795	" 'mask_V' : mask_V ,\n",
796	" 'mask_W' : mask_W ,\n",
797	" 'mask_F' : mask_F ,\n",
798	" 'area_grid_T' : area_grid_T,\n",
799	" 'area_grid_U' : area_grid_U,\n",
800	" 'area_grid_V' : area_grid_V,\n",
801	" 'area_grid_W' : area_grid_W,\n",
802	" 'area_grid_F' : area_grid_F,\n",
803	" 'maskutil_T' : maskutil_T ,\n",
804	" 'maskutil_U' : maskutil_U ,\n",
805	" 'maskutil_V' : maskutil_V ,\n",
806	" 'maskutil_W' : maskutil_W ,\n",
807	" 'maskutil_F' : maskutil_F ,\n",
808	" 'bounds_lon_grid_T': bounds_lon_grid_T,\n",
809	" 'bounds_lat_grid_T': bounds_lat_grid_T,\n",
810	" 'bounds_lon_grid_U': bounds_lon_grid_U,\n",
811	" 'bounds_lat_grid_U': bounds_lat_grid_U,\n",
812	" 'bounds_lon_grid_V': bounds_lon_grid_V,\n",
813	" 'bounds_lat_grid_V': bounds_lat_grid_V,\n",
814	" 'bounds_lon_grid_W': bounds_lon_grid_W,\n",
815	" 'bounds_lat_grid_W': bounds_lat_grid_W,\n",
816	" 'bounds_lon_grid_F': bounds_lon_grid_F,\n",
817	" 'bounds_lat_grid_F': bounds_lat_grid_F,\n",
818	"})\n",
819	"#replace nav_lon nav_lat with variables obtained from NEMO coordinates.nc file\n",
820	"#by construction nav_lon nav_lat are those of the bathymetry \n",
821	"for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
822	" for dir_type in ['lat', 'lon'] :\n",
823	" coord_name = 'nav_' + dir_type + '_grid_' + cd_type\n",
824	" ds.coords[coord_name]=locals()[coord_name]\n",
825	"\n",
826	"ds.attrs['name'] = 'coordinates_mask'\n",
827	"ds.attrs['description'] = 'coordinates and mask for MOSAIX'\n",
828	"ds.attrs['title'] = 'coordinates_mask'\n",
829	"ds.attrs['source'] = 'IPSL Earth system model'\n",
830	"ds.attrs['group'] = 'ICMC IPSL Climate Modelling Center'\n",
831	"ds.attrs['Institution'] = 'IPSL https.//www.ipsl.fr'\n",
832	"ds.attrs['Model'] = model\n",
833	"ds.attrs['timeStamp'] = '{:%Y-%b-%d %H:%M:%S}'.format (datetime.datetime.now ())\n",
834	"ds.attrs['history'] = 'Build from ' + n_coord + ' and ' + n_bathy\n",
835	"ds.attrs['directory'] = os.getcwd ()\n",
836	"ds.attrs['user'] = os.getlogin ()\n",
837	"ds.attrs['HOSTNAME'] = platform.node ()\n",
838	"ds.attrs['Python'] = 'Python version: ' + platform.python_version ()\n",
839	"ds.attrs['xarray'] = 'xarray version: ' + xr.__version__\n",
840	"ds.attrs['OS'] = platform.system ()\n",
841	"ds.attrs['release'] = platform.release ()\n",
842	"ds.attrs['hardware'] = platform.machine ()\n",
843	"\n",
844	"ds.attrs['SVN_Author'] = \"$Author$\"\n",
845	"ds.attrs['SVN_Date'] = \"$Date$\"\n",
846	"ds.attrs['SVN_Revision'] = \"$Revision$\"\n",
847	"ds.attrs['SVN_Id'] = \"$Id$\"\n",
848	"ds.attrs['SVN_HeadURL'] = \"$HeadURL$\""
849	]
850	},
851	{
852	"cell_type": "code",
853	"execution_count": 24,
854	"id": "41391c22",
855	"metadata": {},
856	"outputs": [],
857	"source": [
858	"ds.to_netcdf (n_out)"
859	]
860	},
861	{
862	"cell_type": "markdown",
863	"id": "75143732",
864	"metadata": {},
865	"source": [
866	"# That's all folk's !"
867	]
868	},
869	{
870	"cell_type": "code",
871	"execution_count": null,
872	"id": "88457886",
873	"metadata": {},
874	"outputs": [],
875	"source": []
876	}
877	],
878	"metadata": {
879	"kernelspec": {
880	"display_name": "Python 3",
881	"language": "python",
882	"name": "python3"
883	},
884	"language_info": {
885	"codemirror_mode": {
886	"name": "ipython",
887	"version": 3
888	},
889	"file_extension": ".py",
890	"mimetype": "text/x-python",
891	"name": "python",
892	"nbconvert_exporter": "python",
893	"pygments_lexer": "ipython3",
894	"version": "3.8.10"
895	}
896	},
897	"nbformat": 4,
898	"nbformat_minor": 5
899	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format