Context Navigation

computegrid.pro @ 488

Last change on this file since 488 was 488, checked in by smasson, 11 years ago
update periodicity check in computegrid
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 50.5 KB

Line
1	;+
2	;
3	; @file_comments
4	; compute the grid parameters (cm_4mesh) common
5	;
6	; domains sizes:
7	; ---------------
8	; jpi, jpj, jpk, jpiglo, jpjglo, jpkglo, jpidta, jpjdta, jpkdta
9	;
10	; domains positions regarding to the original grid:
11	; --------------------------------------------------
12	; ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh
13	; ixmindta, ixmaxdta, iymindta, iymaxdta, izmindta, izmaxdta
14	;
15	; horizontal parameters:
16	; ----------------------
17	; glamt, glamf, gphit, gphit, e1t, e2t
18	;
19	; additional horizontal parameters if FULLCGRID keyword is defined:
20	; -----------------------------------------------------------------
21	; glamu, glamv, gphiu, gphiv, e1u, e1v, e1f, e2u, e2v, e2f
22	;
23	; verticals parameters:
24	; ---------------------
25	; gdept, gdepw, e3t, e3w
26	;
27	; masks:
28	; ------
29	; tmask
30	;
31	; additional masks if FULLCGRID keyword is defined:
32	; -------------------------------------------------
33	; umaskred, vmaskred, fmaskredx, fmaskredy
34	;
35	; triangles_list:
36	; ---------------
37	; triangulation
38	;
39	; key_* variables:
40	; ----------------
41	; key_onearth, key_periodic, key_shift, key_stride, key_partialstep,
42	; key_yreverse, key_zreverse, key_gridtype
43	;
44	; xxx related variables:
45	; ----------------------
46	; ccmeshparameters, ccreadparameters
47	;
48	; @categories
49	; Grid
50	;
51	; @param startx {in}{optional}{type=scalar}
52	; x starting point, optional if [XY]AXIS keyword is used
53	;
54	; @param starty {in}{optional}{type=scalar}
55	; y starting point, optional if [XY]AXIS keyword is used
56	;
57	; @param stepxin {in}{optional}{type=scalar or vector}
58	; x direction step, optional if [XY]AXIS keyword is used, must be > 0
59	; if stepxin is a vector nx is not used
60	;
61	; @param stepyin {in}{optional}{type=scalar or vector}
62	; y direction step, optional if [XY]AXIS keyword is used,
63	; could be > 0 (south to north) or < 0 (north to south)
64	; if stepyin is a vector ny is not used
65	;
66	; @param nxin {in}{optional}{type=scalar}
67	; number of points in x direction,
68	; optional if [XY]AXIS keyword is used or stepxin is a vector
69	;
70	; @param nyin {in}{optional}{type=scalar}
71	; number of points in y direction,
72	; optional if [XY]AXIS keyword is used or stepyin is a vector
73	;
74	; @keyword FULLCGRID {default=0}{type=scalar: 0 or 1}
75	; Activate to specify that you want to compute all the C grid parameters:
76	; definition of glam[uv], gphi[uv], e1[uvf], e2[uvf], [uv]maskred and
77	; fmaskred[xy] will be add to the default computations
78	;
79	; @keyword GLAMBOUNDARY {default=those defined in the file}{type=2 elements vector}
80	; Longitude boundaries that should be used to visualize the data.
81	; lon2 > lon1
82	; lon2 - lon1 le 360
83	; By default, the common (cm_4mesh) variable key_shift will be automatically
84	; defined according to GLAMBOUNDARY.
85	;
86	; @keyword MASK {default=array of 1}{type=2D or 3D array}
87	; Specify the land(0)/sea(1) mask
88	;
89	; @keyword ONEARTH {default=1}{type=scalar: 0 or 1}
90	; Force the manual definition of data localization on the earth or not
91	; 0) if the data are not on the earth
92	; 1) if the data are on earth (in that case we can for example use
93	; the labels 'longitude', 'latitude' in plots).
94	; The resulting value will be stored in the common (cm_4mesh) variable key_onearth
95	; ONEARTH=0 forces PERIODIC=0, SHIFT=0 and is cancelling GLAMBOUNDARY
96	;
97	; @keyword PERIODIC {default=computed by using the first line of glamt}{type=scalar: 0 or 1}
98	; Force the manual definition of the grid zonal periodicity.
99	; The resulting value will be stored in the common (cm_4mesh) variable key_periodic
100	; PERIODIC=0 forces SHIFT=0
101	;
102	; @keyword PLAIN {default=0}{type=scalar: 0 or 1}
103	; Force YREVERSE=0, ZREVERSE=0, PERIODIC=0, SHIFT=0, STRIDE=[1, 1, 1] and
104	; suppress the automatic redefinition of the domain in case of x periodicity overlap,
105	; y periodicity overlap (ORCA type only) and mask border to 0.
106	;
107	; @keyword SHIFT {default=computed according to glamboundary}{type=scalar}
108	; Force the manual definition of the zonal shift that must be apply to the data.
109	; The resulting value will be stored in the common (cm_4mesh) variable key_shift
110	; Note that if key_periodic=0 then in any case key_shift=0.
111	;
112	; @keyword STRCALLING {type=string}
113	; a string containing the calling command used to
114	; call computegrid (this is used by <pro>xxx</pro>)
115	;
116	; @keyword STRIDE {default=[1, 1, 1]}{type=3 elements vector}
117	; Specify the stride in x, y and z direction. The resulting
118	; value will be stored in the common (cm_4mesh) variable key_stride
119	;
120	; @keyword XAXIS {type=1D or 2D array}
121	; Specify longitudes in this case startx, stepx and nx are not used but
122	; could be necessary if the y axis is not defined with yaxis. It must be
123	; possible to sort the first line of xaxis in the increasing order by
124	; shifting its elements.
125	;
126	; @keyword YAXIS {type=1D or 2D array}
127	; Specify latitudes in this case starty, stepy and ny are not used but
128	; starty and stepy could be necessary if the x axis is not defined with
129	; xaxis. It must be sorted in the increasing or decreasing order (along each column if 2d array).
130	;
131	; @keyword XYINDEX activate to specify that the horizontal grid should
132	; be simply defined by using the index of the points
133	; (xaxis = findgen(nx) and yaxis = findgen(ny))
134	; using this keyword forces key_onearth=0
135	;
136	; @keyword XMINMESH {default=0L}{type=scalar}
137	; Define common (cm_4mesh) variables ixminmesh used to define the localization
138	; of the first point of the grid along the x direction in a zoom of the original grid
139	;
140	; @keyword YMINMESH {default=0L}{type=scalar}
141	; Define common (cm_4mesh) variables iyminmesh used to define the localization
142	; of the first point of the grid along the y direction in a zoom of the original grid
143	;
144	; @keyword ZMINMESH {default=0L}{type=scalar}
145	; Define common (cm_4mesh) variables izminmesh used to define the localization
146	; of the first point of the grid along the z direction in a zoom of the original grid
147	;
148	; @keyword XMAXMESH {default=jpiglo-1}{type=scalar}
149	; Define common (cm_4mesh) variables ixmaxmesh used to define the localization
150	; of the last point of the grid along the x direction in a zoom of the original grid
151	; Note that if XMAXMESH < 0 then ixmaxmesh is defined as ixmaxmesh = jpiglo -1 + xmaxmesh
152	;
153	; @keyword YMAXMESH {default=jpjglo-1}{type=scalar}
154	; Define common (cm_4mesh) variables iymaxmesh used to define the localization
155	; of the last point of the grid along the y direction in a zoom of the original grid
156	; Note that if YMAXMESH < 0 then iymaxmesh is defined as iymaxmesh = jpjglo -1 + ymaxmesh
157	;
158	; @keyword ZMAXMESH {default=jpkglo-1}{type=scalar}
159	; Define common (cm_4mesh) variables izmaxmesh used to define the localization
160	; of the last point of the grid along the z direction in a zoom of the original grid
161	; Note that if ZMAXMESH < 0 then izmaxmesh is defined as izmaxmesh = jpkglo -1 + maxmesh
162	;
163	; @keyword FBASE2TBASE
164	; Activate when the model is a C grid based on a F point
165	; (with a F point at the bottom-left corner and a T point at the
166	; upper-right corner). In this case, we ignore
167	; - the first line of F and V points
168	; - the last line of T and U points
169	; - if the grid is not x-periodic, the first column of F and U points
170	; - if the grid is not x-periodic, the last column of T and V points.
171	; => we are back to a C grid based on T point as for OPA model.
172	; Note that in that case, key_gridtype = 'c_f' and not 'c' (-> used in read_ncdf)
173	; Note that activate FBASE2TBASE forces FULLCGRID=1
174	;
175	; @keyword UBASE2TBASE
176	; Activate when the model is a C grid based on a U point
177	; (with a U point at the bottom-left corner and a T point at the
178	; upper-right corner). In this case, we ignore
179	; - if the grid is not x-periodic, the first column of F and U points
180	; - if the grid is not x-periodic, the last column of T and V points.
181	; => we are back to a C grid based on T point as for OPA model.
182	; Note that in that case, key_gridtype = 'c_u' and not 'c' (-> used in read_ncdf)
183	; Note that activate UBASE2TBASE forces FULLCGRID=1
184	;
185	; @keyword VBASE2TBASE
186	; Activate when the model is a C grid based on a V point
187	; (with a V point at the bottom-left corner and a T point at the
188	; upper-right corner). In this case, we ignore
189	; - the first line of F and V points
190	; - the last line of T and U points
191	; => we are back to a C grid based on T point as for OPA model.
192	; Note that in that case, key_gridtype = 'c_v' and not 'c' (-> used in read_ncdf)
193	; Note that activate VBASE2TBASE forces FULLCGRID=1
194	;
195	; @keyword ROMSH {type=2D array}
196	; This array is the final bathymetry at RHO-points. It is stored in the common
197	; variable (cm_4mesh) romszinfos.h
198	; Used when the model is a ROMS C-grid with one more point
199	; in longitude for T and V grid and one more point in latitude
200	; for T and U grid. In this case, we ignore
201	; - the last line of T and U points
202	; - the last column of T and V points.
203	; => we are back to a C grid based on T point as for OPA model.
204	; Note that activate ROMSH forces FULLCGRID=1
205	;
206	; @keyword STRCALLING {type=scalar string}
207	; Used by xxx...
208	;
209	; @keyword YREVERSE {default=computed according to gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]}{type=scalar}
210	; Force the manual definition of the y reverse that must be apply to the data.
211	; The resulting value will be stored in the common (cm_4mesh) variable key_yreverse
212	;
213	; @keyword ZAXIS {type=1D}
214	; Specify the vertical axis. Must be sorted in the increasing or decreasing order
215	;
216	; @keyword ZREVERSE {default=computed according to gdept[0] GT gdept[1]}{type=scalar}
217	; Force the manual definition of the z reverse that must be apply to the data.
218	; The resulting value will be stored in the common (cm_4mesh) variable key_zreverse
219	;
220	; @keyword _EXTRA
221	; not used in the present case ...
222	;
223	; @uses
224	; <pro>cm_4mesh</pro>
225	; <pro>cm_4data</pro>
226	; <pro>cm_4cal</pro>
227	;
228
229	; @restrictions
230	; if the grid has x/y periodicity overlap and/or if
231	; the mask has 0 everywhere at the border (like a closed sea) and
232	; if (we did not activate /plain and xminmesh, xmaxmesh, yminmesh,
233	; ymaxmesh keywords are defined to their default values), we redefine
234	; xminmesh, xmaxmesh, yminmesh, ymaxmesh in order to remove the
235	; overlapping part and/or to open the domain (avoid it be forced
236	; to use cell_fill = 1).
237	;
238	; FUV points definition is not exact if the grid is irregular
239	;
240	; @history
241	; Sebastien Masson (smasson\@lodyc.jussieu.fr)
242	; 2000-04-20
243	; Sept 2004, several bug fix to suit C grid type...
244	; Aug 2005, rewritte almost everything...
245	;
246	; @version
247	; $Id$
248	;
249	;-
250	PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $
251	, XAXIS=xaxis, YAXIS=yaxis, ZAXIS=zaxis $
252	, MASK=mask, GLAMBOUNDARY=glamboundary $
253	, XMINMESH=xminmesh, XMAXMESH=xmaxmesh $
254	, YMINMESH=yminmesh, YMAXMESH=ymaxmesh $
255	, ZMINMESH=zminmesh, ZMAXMESH=zmaxmesh $
256	, ONEARTH=onearth, PERIODIC=periodic $
257	, PLAIN=plain, SHIFT=shift, STRIDE=stride $
258	, YREVERSE=yreverse, ZREVERSE=zreverse $
259	, FULLCGRID=fullcgrid, XYINDEX=xyindex $
260	, UBASE2TBASE=ubase2tbase, VBASE2TBASE=vbase2tbase $
261	, FBASE2TBASE=fbase2tbase $
262	, STRCALLING=strcalling, ROMSH=romsh, _EXTRA=ex
263	;
264	compile_opt idl2, strictarrsubs
265	;
266	@cm_4mesh
267	@cm_4data
268	@cm_4cal
269	IF NOT keyword_set(key_forgetold) THEN BEGIN
270	@updatenew
271	@updatekwd
272	ENDIF
273	;---------------------------------------------------------
274	;------------------------------------------------------------
275	time1 = systime(1) ; for key_performance
276	;------------------------------------------------------------
277	;
278	;====================================================
279	; Check input parameters
280	;====================================================
281	;
282	; xaxis related parameters
283	;
284	if n_elements(xaxis) NE 0 then BEGIN
285	CASE (size(xaxis))[0] OF
286	0:nx = 1L
287	1:nx = (size(xaxis))[1]
288	2:nx = (size(xaxis))[1]
289	ENDCASE
290	ENDIF ELSE BEGIN
291	IF n_elements(startx) EQ 0 THEN BEGIN
292	dummy = report('If xaxis is not given, startx must be defined')
293	return
294	ENDIF
295	CASE n_elements(stepxin) OF
296	0:BEGIN
297	dummy = report('If xaxis is not given, stepxin must be defined')
298	return
299	END
300	1:BEGIN
301	IF n_elements(nxin) EQ 0 THEN BEGIN
302	dummy = report('If xaxis is not given and stepxin has only one element, nx must be defined')
303	return
304	ENDIF ELSE nx = nxin
305	END
306	ELSE:nx = n_elements(stepxin)
307	ENDCASE
308	ENDELSE
309	;
310	; yaxis related parameters
311	;
312	if n_elements(yaxis) NE 0 then BEGIN
313	CASE (size(yaxis))[0] OF
314	0:ny = 1L
315	1:ny = (size(yaxis))[1]
316	2:ny = (size(yaxis))[2]
317	ENDCASE
318	ENDIF ELSE BEGIN
319	IF n_elements(starty) EQ 0 THEN BEGIN
320	dummy = report('If yaxis is not given, starty must be defined')
321	return
322	ENDIF
323	CASE n_elements(stepyin) OF
324	0:BEGIN
325	dummy = report('If yaxis is not given, stepyin must be defined')
326	return
327	END
328	1:BEGIN
329	IF n_elements(nyin) EQ 0 THEN BEGIN
330	dummy = report('If yaxis is not given and stepyin has only one element, ny must be defined')
331	return
332	ENDIF ELSE ny = nyin
333	END
334	ELSE:ny = n_elements(stepyin)
335	ENDCASE
336	ENDELSE
337	;
338	; zaxis related parameters
339	;
340	if n_elements(zaxis) NE 0 then BEGIN
341	CASE (size(zaxis))[0] OF
342	0:nz = 1L
343	1:nz = (size(zaxis))[1]
344	ELSE:BEGIN
345	ras = report( 'not coded')
346	stop
347	END
348	ENDCASE
349	ENDIF ELSE nz = 1L
350	;
351	;====================================================
352	; Others automatic definitions...
353	;====================================================
354	;
355	IF keyword_set(romsh) THEN fullcgrid = 1
356	CASE 1 OF
357	keyword_set(fbase2tbase):key_gridtype = 'c_f'
358	keyword_set(ubase2tbase):key_gridtype = 'c_u'
359	keyword_set(vbase2tbase):key_gridtype = 'c_v'
360	else:key_gridtype = 'c'
361	ENDCASE
362	IF strlen(key_gridtype) EQ 3 THEN fullcgrid = 1
363	;
364	IF n_elements(xminmesh) NE 0 OR n_elements(xmaxmesh) NE 0 THEN BEGIN
365	IF n_elements(xminmesh) EQ 0 THEN xminmesh = ixminmesh
366	IF n_elements(xmaxmesh) EQ 0 THEN xmaxmesh = ixmaxmesh
367	IF nx EQ jpi AND xminmesh EQ ixminmesh AND xmaxmesh EQ ixmaxmesh THEN xalreadycut = 1
368	ENDIF
369	IF keyword_set(xalreadycut) THEN BEGIN
370	xmin = 0
371	xmax = jpi - 1
372	nxx = jpi
373	ENDIF ELSE BEGIN
374	jpiglo = long(nx)
375	IF n_elements(xminmesh) NE 0 THEN ixminmesh = long(xminmesh[0]) ELSE ixminmesh = 0l
376	IF n_elements(xmaxmesh) NE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh = jpiglo-1
377	IF ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh
378	ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1)
379	ixminmesh = 0 > ixminmesh < ixmaxmesh
380	jpi = ixmaxmesh-ixminmesh+1
381	xmin = ixminmesh
382	xmax = ixmaxmesh
383	nxx = jpiglo
384	ENDELSE
385
386	IF n_elements(yminmesh) NE 0 OR n_elements(ymaxmesh) NE 0 THEN BEGIN
387	IF n_elements(yminmesh) EQ 0 THEN yminmesh = iyminmesh
388	IF n_elements(ymaxmesh) EQ 0 THEN ymaxmesh = iymaxmesh
389	IF ny EQ jpj AND yminmesh EQ iyminmesh AND ymaxmesh EQ iymaxmesh THEN yalreadycut = 1
390	ENDIF
391	IF keyword_set(yalreadycut) THEN BEGIN
392	ymin = 0
393	ymax = jpj - 1
394	nyy = jpj
395	ENDIF ELSE BEGIN
396	jpjglo = long(ny)
397	IF n_elements(yminmesh) NE 0 THEN iyminmesh = long(yminmesh[0]) ELSE iyminmesh = 0l
398	IF n_elements(ymaxmesh) NE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh = jpjglo-1
399	IF key_gridtype EQ 'c_v' OR key_gridtype EQ 'c_f' THEN iymaxmesh = iymaxmesh-1
400	IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh
401	iymaxmesh = 0 > iymaxmesh < (jpjglo-1)
402	iyminmesh = 0 > iyminmesh < iymaxmesh
403	jpj = iymaxmesh-iyminmesh+1
404	ymin = iyminmesh
405	ymax = iymaxmesh
406	nyy = jpjglo
407	ENDELSE
408
409	IF n_elements(zminmesh) NE 0 AND n_elements(zmaxmesh) NE 0 THEN BEGIN
410	IF nz EQ jpk AND zminmesh EQ izminmesh AND zmaxmesh EQ izmaxmesh THEN zalreadycut = 1
411	ENDIF
412	IF keyword_set(zalreadycut) THEN BEGIN
413	zmin = 0
414	zmax = jpk - 1
415	nzz = jpk
416	ENDIF ELSE BEGIN
417	jpkglo = long(nz)
418	IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh = 0l
419	IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh = jpkglo-1
420	IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh
421	izmaxmesh = 0 > izmaxmesh < (jpkglo-1)
422	izminmesh = 0 > izminmesh < izmaxmesh
423	jpk = izmaxmesh-izminmesh+1
424	zmin = izminmesh
425	zmax = izmaxmesh
426	nzz = jpkglo
427	ENDELSE
428	;
429	; impact of plain keyword:
430	;
431	IF keyword_set(plain) THEN BEGIN
432	yreverse = 0
433	zreverse = 0
434	periodic = 0
435	shift = 0
436	stride = [1, 1, 1]
437	ENDIF
438	;
439	; avoid basics errors...
440	;
441	jpidta = jpiglo
442	jpjdta = jpjglo
443	jpkdta = jpkglo
444	ixmindta = 0
445	ixmaxdta = jpidta-1
446	iymindta = 0
447	iymaxdta = jpjdta-1
448	izmindta = 0
449	izmaxdta = jpkdta-1
450	;
451	key_partialstep = 0
452	if n_elements(stride) eq 3 then key_stride = stride $
453	ELSE key_stride = [1, 1, 1]
454	;
455	; check xyindex and its consequences
456	;
457	if keyword_set(xyindex) then onearth = 0
458	;
459	; check onearth and its consequences
460	;
461	IF n_elements(onearth) EQ 0 THEN key_onearth = 1b $
462	ELSE key_onearth = keyword_set(onearth)
463	IF NOT key_onearth THEN BEGIN
464	periodic = 0
465	shift = 0
466	ENDIF
467
468	r = 6371000.
469	;
470	;====================================================
471	; X direction : glamt
472	;====================================================
473	;
474	; def of glamt
475	;
476	if n_elements(xaxis) NE 0 then BEGIN
477	if keyword_set(xyindex) THEN glamt = findgen(jpiglo) ELSE glamt = xaxis
478	ENDIF ELSE BEGIN
479	if keyword_set(xyindex) THEN stepx = 1. ELSE stepx = stepxin
480	CASE 1 OF
481	n_elements(stepx):glamt = startx + findgen(jpiglo)*stepx
482	size(stepx, /n_dimensions):glamt = startx + total(stepx, /cumulative)
483	ELSE:BEGIN
484	dummy = report('Wrong definition of stepx...')
485	return
486	END
487	ENDCASE
488	ENDELSE
489	;
490	; apply glamboundary
491	;
492	IF keyword_set(glamboundary) AND key_onearth THEN BEGIN
493	IF glamboundary[0] GE glamboundary[1] THEN stop
494	IF glamboundary[1]-glamboundary[0] GT 360 THEN stop
495	glamt = glamt MOD 360
496	smaller = where(glamt LT glamboundary[0])
497	if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360.
498	bigger = where(glamt GE glamboundary[1])
499	if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360.
500	ENDIF
501	;
502	; force glamt to have 2 dimensions
503	;
504	IF n_elements(glamt) EQ nxx*nyy THEN glamt = reform(glamt, nxx, nyy, /over) $
505	ELSE glamt = reform(glamt, nxx, /over)
506	CASE size(glamt, /n_dimensions) OF
507	1:BEGIN
508	IF n_elements(glamt) EQ 1 THEN glamt = replicate(glamt[0], jpi, jpj) $
509	ELSE glamt = glamt[xmin:xmax]#replicate(1, jpj)
510	END
511	2:glamt = glamt[xmin:xmax, ymin:ymax]
512	ENDCASE
513	; keep 2d array even with degenerated dimension
514	IF jpj EQ 1 THEN glamt = reform(glamt, jpi, jpj, /over)
515	;
516	;====================================================
517	; Y direction : gphit
518	;====================================================
519	;
520	; def of gphit
521	;
522	if n_elements(yaxis) NE 0 THEN BEGIN
523	if keyword_set(xyindex) THEN gphit = findgen(jpjglo) ELSE gphit = yaxis
524	ENDIF ELSE BEGIN
525	if keyword_set(xyindex) THEN stepy = 1. ELSE stepy = stepyin
526	CASE 1 OF
527	n_elements(stepy):gphit = starty + findgen(jpjglo)*stepy
528	size(stepy, /n_dimensions):gphit = starty + total(stepy, /cumulative)
529	ELSE:BEGIN
530	dummy = report('Wrong definition of stepy...')
531	return
532	END
533	ENDCASE
534	ENDELSE
535	;
536	; force gphit to have 2 dimensions
537	;
538	IF n_elements(gphit) EQ nxx*nyy THEN gphit = reform(gphit, nxx, nyy, /over) $
539	ELSE gphit = reform(gphit, nyy, /over)
540	CASE size(gphit, /n_dimensions) OF
541	1:BEGIN
542	IF n_elements(gphit) EQ 1 THEN gphit = replicate(gphit[0], jpi, jpj) $
543	ELSE gphit = replicate(1, jpi)#gphit[ymin:ymax]
544	END
545	2:gphit = gphit[xmin:xmax, ymin:ymax]
546	ENDCASE
547	; keep 2d array even with degenerated dimension
548	IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over)
549	;
550	;====================================================
551	; check y periodicity... Only according to ORCA grid
552	;====================================================
553	; check the periodicity if iyminmesh and iymaxmesh have the default definitions...
554	IF NOT keyword_set(plain) AND key_onearth EQ 1 $
555	AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 THEN BEGIN
556
557	CASE 1 OF
558	ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
559	AND array_equal(gphit[1:, jpj-1], reverse(gphit[1:, jpj-3])) EQ 1:BEGIN
560	; T pivot
561	ymaxmesh = -1
562	recall = 1
563	END
564	ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $
565	AND array_equal(gphit[, jpj-1], reverse(shift(gphit[, jpj-3], -1))) EQ 1:BEGIN
566	; T pivot
567	ymaxmesh = -1
568	recall = 1
569	END
570	ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
571	AND array_equal(gphit[, jpj-1], reverse(gphit[, jpj-2])) EQ 1:BEGIN
572	; F pivot
573	ymaxmesh = -1
574	recall = 1
575	END
576	ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $
577	AND array_equal(gphit[, jpj-1], reverse(gphit[, jpj-2])) EQ 1:BEGIN
578	; F pivot
579	ymaxmesh = -1
580	recall = 1
581	END
582	ELSE:
583	ENDCASE
584	ENDIF
585	;
586	;====================================================
587	; check x periodicity...
588	;====================================================
589	IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic)
590	; check the periodicity if ixminmesh and ixmaxmesh have the default definitions...
591	IF NOT keyword_set(plain) AND NOT keyword_set(forcenoperio) AND key_onearth EQ 1 $
592	AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 AND jpi GE 3 THEN BEGIN
593	CASE 0 OF
594	total((glamt[0, ] - glamt[jpi-2, ]) MOD 360) $
595	+ total((glamt[1, ] - glamt[jpi-1, ]) MOD 360):BEGIN
596	xminmesh = 1
597	xmaxmesh = -1
598	recall = 1
599	END
600	total((glamt[0, ] - glamt[jpi-2, ]) MOD 360):BEGIN
601	xminmesh = 1
602	recall = 1
603	END
604	total((glamt[1, ] - glamt[jpi-1, ]) MOD 360):BEGIN
605	xmaxmesh = -1
606	recall = 1
607	END
608	total((glamt[0, ] - glamt[jpi-1, ]) MOD 360):BEGIN
609	xmaxmesh = -1
610	recall = 1
611	END
612	ELSE:
613	ENDCASE
614	ENDIF
615	;====================================================
616	; recall computegrid if needed...
617	;====================================================
618	IF keyword_set(recall) THEN BEGIN
619	computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $
620	, MASK = mask, GLAMBOUNDARY = glamboundary $
621	, XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
622	, YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
623	, ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
624	, PERIODIC = periodic, SHIFT = shift, STRIDE = stride $
625	, FULLCGRID = fullcgrid, XYINDEX = xyindex $
626	, STRCALLING = strcalling $
627	, ROMSH = romsh, _extra = ex
628	return
629	ENDIF
630	;
631	;====================================================
632	; def key_yreverse
633	;====================================================
634	;
635	IF n_elements(yreverse) EQ 0 THEN BEGIN
636	IF jpj GT 1 THEN BEGIN
637	IF total(gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]) GT jpj/2 THEN key_yreverse = 1 ELSE key_yreverse = 0
638	ENDIF ELSE key_yreverse = 0
639	ENDIF ELSE key_yreverse = yreverse
640	IF keyword_set(key_yreverse) THEN BEGIN
641	gphit = reverse(gphit, 2)
642	glamt = reverse(glamt, 2)
643	ENDIF
644	;
645	;====================================================
646	; def of key_shift
647	;====================================================
648	;
649	; definition of key_shift by shifting the array to have the min
650	; values of glamt[*, 0] in glamt[0, 0]
651	;
652	IF n_elements(shift) EQ 0 THEN BEGIN
653	IF jpi GT 1 then BEGIN
654	xtest = glamt[*, 0]
655	key_shift = (where(xtest EQ min(xtest)))[0]
656	IF key_shift NE 0 THEN key_shift = jpi - key_shift
657	ENDIF ELSE key_shift = 0
658	ENDIF ELSE key_shift = shift
659	;
660	;====================================================
661	; def of key_periodic
662	;====================================================
663	;
664	IF n_elements(periodic) EQ 0 THEN BEGIN
665	IF jpi GT 1 THEN BEGIN
666	xtest = shift(glamt[*, 0], key_shift)
667	; check that xtest is now sorted in the increasing order
668	IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN
669	ras = report(['WARNING: we cannot sort the xaxis with a simple shift...', $
670	'we force key_periodic = 0 and key_shift = 0', $
671	'only horizontal plot may be ok...'])
672	key_periodic = 0
673	xnotsorted = 1
674	ENDIF ELSE BEGIN
675	key_periodic = (xtest[jpi-1]+2*(xtest[jpi-1]-xtest[jpi-2])) $
676	GE (xtest[0]+360)
677	ENDELSE
678	ENDIF ELSE key_periodic = 0
679	ENDIF ELSE key_periodic = keyword_set(periodic)
680	;
681	; update key_shift
682	;
683	key_shift = key_shift * (key_periodic EQ 1)
684	;
685	IF (key_gridtype EQ 'c_u' OR key_gridtype EQ 'c_f') AND NOT keyword_set(key_periodic) THEN BEGIN
686	ixmaxmesh = ixmaxmesh-1
687	jpi = jpi-1
688	ENDIF
689	;
690	;====================================================
691	; apply key_shift
692	;====================================================
693	;
694	if keyword_set(key_shift) then BEGIN
695	glamt = shift(glamt, key_shift, 0)
696	gphit = shift(gphit, key_shift, 0)
697	IF jpj EQ 1 THEN BEGIN
698	glamt = reform(glamt, jpi, jpj, /over)
699	gphit = reform(gphit, jpi, jpj, /over)
700	ENDIF
701	ENDIF
702	;
703	;====================================================
704	; Are we using a "regular" grid (that can be described
705	; with x vector and y vector)?
706	;====================================================
707	;
708	; to get faster, we first test the most basic cases before
709	; testing the full array.
710	;
711	CASE 1 OF
712	keyword_set(xyindex):key_irregular = 0b
713	jpi EQ 1 OR jpj EQ 1:key_irregular = 0b
714	n_elements(xaxis) EQ 0 AND n_elements(yaxis) EQ 0:key_irregular = 0b
715	size(reform(xaxis), /n_dimensions) EQ 1 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b
716	n_elements(xaxis) EQ 0 AND size(reform(yaxis), /n_dimensions) EQ 1:key_irregular = 0b
717	n_elements(yaxis) EQ 0 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b
718	array_equal(glamt[, 0], glamt[, jpj-1]) EQ 0:key_irregular = 1b
719	array_equal(gphit[0, ], gphit[jpi-1, ]) EQ 0:key_irregular = 1b
720	array_equal(glamt, glamt[*, 0]#replicate(1, jpj)) EQ 0:key_irregular = 1b
721	array_equal(gphit, replicate(1, jpi)#(gphit[0, ])[]) EQ 0:key_irregular = 1b
722	ELSE:key_irregular = 0b
723	ENDCASE
724	;
725	;====================================================
726	; def of glamf: defined as the middle of T(i,j) T(i+1,j+1)
727	;====================================================
728	;
729	IF jpi GT 1 THEN BEGIN
730	; we must compute stepxf: x distance between T(i,j) T(i+1,j+1)
731	CASE 1 OF
732	n_elements(stepx):stepxf = stepx
733	size(stepx, /n_dimensions):stepxf = stepx#replicate(1, jpj)
734	ELSE:BEGIN
735	if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
736	OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
737	stepxf = (glamt + 720) MOD 360
738	IF jpj EQ 1 THEN stepxf = reform(stepxf, jpi, jpj, /over)
739	stepxf = shift(stepxf, -1, -1) - stepxf
740	stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ]
741	stepxf = min(abs(stepxf), dimension = 3)
742	IF NOT keyword_set(key_periodic) THEN $
743	stepxf[jpi-1, ] = stepxf[jpi-2, ]
744	ENDIF ELSE BEGIN
745	stepxf = shift(glamt, -1, -1) - glamt
746	IF keyword_set(key_periodic) THEN $
747	stepxf[jpi-1, ] = 360 + stepxf[jpi-1, ] $
748	ELSE stepxf[jpi-1, ] = stepxf[jpi-2, ]
749	ENDELSE
750	IF jpj GT 1 THEN BEGIN
751	stepxf[, jpj-1] = stepxf[, jpj-2]
752	stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2]
753	ENDIF
754	END
755	ENDCASE
756	glamf = glamt + 0.5 * stepxf
757	ENDIF ELSE glamf = glamt + 0.5
758	;
759	IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN
760	IF NOT keyword_set(glamboundary) THEN BEGIN
761	bigger = where(glamf GE min(glamt)+360)
762	glamf[bigger] = glamf[bigger]-360.
763	ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
764	ENDIF
765	;
766	IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over)
767	;
768	;====================================================
769	; def of gphif: defined as the middle of T(i,j) T(i+1,j+1)
770	;====================================================
771	;
772	IF jpj GT 1 THEN BEGIN
773	; we must compute stepyf: y distance between T(i,j) T(i+1,j+1)
774	CASE 1 OF
775	n_elements(stepy):stepyf = stepy
776	size(stepy, /n_dimensions):stepyf = replicate(1, jpi)#stepy
777	ELSE:BEGIN
778	stepyf = shift(gphit, -1, -1) - gphit
779	stepyf[, jpj-1] = stepyf[, jpj-2]
780	IF jpi GT 1 THEN BEGIN
781	if NOT keyword_set(key_periodic) THEN $
782	stepyf[jpi-1, ] = stepyf[jpi-2, ]
783	stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
784	ENDIF
785	END
786	ENDCASE
787	gphif = gphit + 0.5 * stepyf
788	ENDIF ELSE gphif = gphit + 0.5
789	IF key_onearth THEN gphif = -90. > gphif < 90.
790	;
791	IF jpj EQ 1 THEN gphif = reform(gphif, jpi, jpj, /over)
792	;
793	;====================================================
794	; e1t: x distance between U(i-1,j) and U(i,j)
795	;====================================================
796	;
797	; -\|----\|---*---\|
798	;
799	IF jpi GT 1 THEN BEGIN
800	IF n_elements(stepx) NE 1 THEN BEGIN
801	IF keyword_set(irregular) THEN BEGIN
802	; we must compute stepxu: x distance between T(i,j) T(i+1,j)
803	IF keyword_set(key_periodic) THEN BEGIN
804	stepxu = (glamt + 720) MOD 360
805	stepxu = shift(stepxu, -1, 0) - stepxu
806	stepxu = [ [[stepxu]], [[stepxu + 360]], [[stepxu - 360]] ]
807	stepxu = min(abs(stepxu), dimension = 3)
808	ENDIF ELSE BEGIN
809	stepxu = shift(glamt, -1, 0) - glamt
810	stepxu[jpi-1, ] = stepxf[jpi-2, ]
811	ENDELSE
812	ENDIF ELSE stepxu = stepxf
813	IF jpj EQ 1 THEN stepxu = reform(stepxu, jpi, jpj, /over)
814	e1t = 0.5*(stepxu+shift(stepxu, 1, 0))
815	IF NOT keyword_set(key_periodic) THEN $
816	e1t[0, ] = e1t[1, ]
817	ENDIF ELSE e1t = replicate(stepx, jpi, jpj)
818	ENDIF ELSE e1t = replicate(1b, jpi, jpj)
819	;
820	IF jpj EQ 1 THEN e1t = reform(e1t, jpi, jpj, /over)
821	;
822	;====================================================
823	; e2t: y distance between V(i,j-1) and V(i,j)
824	;====================================================
825	;
826	IF jpj GT 1 THEN BEGIN
827	; we must compute stepyv: y distance between T(i,j) T(i,j+1)
828	IF n_elements(stepy) NE 1 THEN BEGIN
829	IF keyword_set(key_irregular) THEN BEGIN
830	stepyv = shift(gphit, 0, -1) - gphit
831	stepyv[, jpj-1] = stepyv[, jpj-2]
832	ENDIF ELSE stepyv = stepyf
833	e2t = 0.5*(stepyv+shift(stepyv, 0, 1))
834	e2t[, 0] = e2t[, 1]
835	ENDIF ELSE e2t = replicate(stepy, jpi, jpj)
836	ENDIF ELSE e2t = replicate(1b, jpi, jpj)
837	;
838	IF key_onearth THEN e2t = r * !pi/180. * temporary(e2t)
839	;
840	IF jpj EQ 1 THEN e2t = reform(e2t, jpi, jpj, /over)
841	;
842	;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
843	IF keyword_set(fullcgrid) THEN BEGIN
844	;
845	;====================================================
846	; def of glamu: defined as the middle of T(i,j) T(i+1,j)
847	;====================================================
848	;
849	IF keyword_set(irregular) THEN BEGIN
850	glamu = glamt + 0.5 * stepxu
851	IF keyword_set(glamboundary) AND key_onearth THEN $
852	glamu = glamboundary[0] > temporary(glamu) < glamboundary[1]
853	ENDIF ELSE glamu = glamf
854	;
855	IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over)
856	;
857	;====================================================
858	; def of gphiu: defined as the middle of T(i,j) T(i+1,j)
859	;====================================================
860	;
861	IF jpi GT 1 THEN BEGIN
862	; we must compute stepyu: y distance between T(i+1,j) T(i,j)
863	IF keyword_set(key_irregular) THEN BEGIN
864	stepyu = shift(gphit, -1, 0) - gphit
865	IF NOT keyword_set(key_periodic) THEN $
866	stepyu[jpi-1, ] = stepyu[jpi-2, ]
867	gphiu = gphit + 0.5 * stepyu
868	ENDIF ELSE gphiu = gphit
869	ENDIF ELSE gphiu = gphit
870	IF key_onearth THEN gphiu = -90. > gphiu < 90.
871	;
872	IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over)
873	;
874	;====================================================
875	; def of glamv: defined as the middle of T(i,j) T(i,j+1)
876	;====================================================
877	;
878	IF jpj GT 1 THEN BEGIN
879	; we must compute stepxv: x distance between T(i,j) T(i,j+1)
880	IF keyword_set(irregular) THEN BEGIN
881	IF keyword_set(key_periodic) THEN BEGIN
882	stepxv = (glamt + 720) MOD 360
883	stepxv = shift(stepxv, 0, -1) - stepxv
884	stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ]
885	stepxv = min(abs(stepxv), dimension = 3)
886	ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt
887	stepxv[, jpj-1] = stepxv[, jpj-2]
888	glamv = glamt + 0.5 * stepxv
889	IF keyword_set(glamboundary) AND key_onearth THEN $
890	glamv = glamboundary[0] > temporary(glamv) < glamboundary[1]
891	ENDIF ELSE glamv = glamt
892	ENDIF ELSE glamv = glamt
893	;
894	;====================================================
895	; def of gphiv: defined as the middle of T(i,j) T(i,j+1)
896	;====================================================
897	;
898	IF keyword_set(key_irregular) THEN $
899	gphiv = gphit + 0.5 * stepyv $
900	ELSE gphiv = gphif
901	IF key_onearth THEN gphiv = -90. > gphiv < 90.
902	;
903	IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over)
904	;
905	;====================================================
906	; e1u: x distance between T(i,j) and T(i+1,j)
907	;====================================================
908	;
909	IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $
910	e1u = stepxu ELSE e1u = e1t
911	;
912	;====================================================
913	; e2u: y distance between F(i,j-1) and F(i,j)
914	;====================================================
915	;
916	IF keyword_set(key_irregular) THEN BEGIN
917	e2u = gphif - shift(gphif, 0, 1)
918	e2u[, 0] = e2u[, 1]
919	IF key_onearth THEN e2u = r * !pi/180. * temporary(e2u)
920	ENDIF ELSE e2u = e2t
921	;
922	IF jpj EQ 1 THEN e2u = reform(e2u, jpi, jpj, /over)
923	;
924	;====================================================
925	; e1v: x distance between F(i-1,j) and F(i,j)
926	;====================================================
927	;
928	IF keyword_set(irregular) THEN BEGIN
929	IF keyword_set(key_periodic) THEN BEGIN
930	e1v = (glamf + 720) MOD 360
931	e1v = e1v - shift(e1v, 1, 0)
932	e1v = [ [[e1v]], [[e1v + 360]], [[e1v - 360]] ]
933	e1v = min(abs(e1v), dimension = 3)
934	ENDIF ELSE BEGIN
935	e1v = glamf - shift(glamf, 1, 0)
936	e1v[0, ] = stepxf[1, ]
937	ENDELSE
938	ENDIF ELSE e1v = e1t
939	;
940	IF jpj EQ 1 THEN e1v = reform(e1v, jpi, jpj, /over)
941	;
942	;====================================================
943	; e2v: y distance between T(i,j) and T(i+1,j)
944	;====================================================
945	;
946	IF jpj GT 1 and n_elements(stepy) NE 1 THEN BEGIN
947	e2v = stepyv
948	IF key_onearth THEN e2v = r * !pi/180. * temporary(e2v)
949	ENDIF ELSE e2v = e2t
950	;
951	;====================================================
952	; e1f: x distance between V(i,j) and V(i+1,j)
953	;====================================================
954	;
955	IF keyword_set(irregular) THEN BEGIN
956	IF keyword_set(key_periodic) THEN BEGIN
957	e1f = (glamv + 720) MOD 360
958	e1f = shift(e1f, -1, 0) - e1f
959	e1f = [ [[e1f]], [[e1f + 360]], [[e1f - 360]] ]
960	e1f = min(abs(e1f), dimension = 3)
961	ENDIF ELSE BEGIN
962	e1f = shift(glamv, -1, 0) - glamt
963	e1f[jpi-1, ] = stepxf[jpi-2, ]
964	ENDELSE
965	ENDIF ELSE e1f = e1u
966	;
967	IF jpj EQ 1 THEN e1f = reform(e1f, jpi, jpj, /over)
968	;
969	;====================================================
970	; e2f: y distance between U(i,j) and U(i,j+1)
971	;====================================================
972	;
973	IF keyword_set(key_irregular) THEN BEGIN
974	e2f = shift(gphiu, 0, -1) - gphiu
975	e2f[, jpj-1] = e2f[, jpj-2]
976	IF key_onearth THEN e2f = r * !pi/180. * temporary(e2f)
977	ENDIF ELSE e2f = e2v
978	;
979	IF jpj EQ 1 THEN e2f = reform(e2f, jpi, jpj, /over)
980	;
981	ENDIF
982	;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
983	;
984	;
985	;====================================================
986	; e1[tuvf] from degree to meters
987	;====================================================
988	;
989	IF keyword_set(key_onearth) THEN BEGIN
990	e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit)
991	IF keyword_set(fullcgrid) THEN BEGIN
992	e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu)
993	e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv)
994	e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif)
995	ENDIF
996	ENDIF
997	;
998	IF jpj EQ 1 THEN BEGIN
999	e1t = reform(e1t, jpi, jpj, /over)
1000	IF keyword_set(fullcgrid) THEN BEGIN
1001	e1u = reform(e1u, jpi, jpj, /over)
1002	e1v = reform(e1v, jpi, jpj, /over)
1003	e1f = reform(e1f, jpi, jpj, /over)
1004	ENDIF
1005	ENDIF
1006	;
1007	;====================================================
1008	; if not fullcgrid: make sure we don't use glam[uv], gphi[uv], e[12][uvf]
1009	;====================================================
1010	;
1011	IF NOT keyword_set(fullcgrid) THEN BEGIN
1012	glamu = !values.f_nan & glamv = !values.f_nan
1013	gphiu = !values.f_nan & gphiv = !values.f_nan
1014	e1u = !values.f_nan & e1v = !values.f_nan & e1f = !values.f_nan
1015	e2u = !values.f_nan & e2v = !values.f_nan & e2f = !values.f_nan
1016	firstxu = !values.f_nan & lastxu = !values.f_nan & nxu = !values.f_nan
1017	firstyu = !values.f_nan & lastyu = !values.f_nan & nyu = !values.f_nan
1018	firstxv = !values.f_nan & lastxv = !values.f_nan & nxv = !values.f_nan
1019	firstyv = !values.f_nan & lastyv = !values.f_nan & nyv = !values.f_nan
1020	ENDIF
1021	;
1022	;====================================================
1023	; Z direction
1024	;====================================================
1025	;
1026	; z axis
1027	;
1028	CASE n_elements(zaxis) OF
1029	0:BEGIN
1030	gdept = 0.
1031	key_zreverse = 0
1032	END
1033	1:BEGIN
1034	gdept = zaxis
1035	key_zreverse = 0
1036	END
1037	ELSE:BEGIN
1038	gdept = zaxis[zmin:zmax]
1039	IF n_elements(zreverse) EQ 0 THEN BEGIN
1040	IF jpk GT 1 THEN BEGIN
1041	if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0
1042	ENDIF ELSE key_zreverse = 0
1043	ENDIF ELSE key_zreverse = zreverse
1044	IF keyword_set(key_zreverse) THEN gdept = reverse(gdept)
1045	END
1046	ENDCASE
1047	;
1048	if n_elements(gdept) GT 1 then BEGIN
1049	stepz = shift(gdept, -1)-gdept
1050	stepz[jpk-1] = stepz[jpk-2]
1051	gdepw = 0. > (gdept-stepz/2.)
1052	ENDIF ELSE BEGIN
1053	stepz = 1.
1054	gdepw = gdept
1055	ENDELSE
1056	IF keyword_set(romsh) THEN gdepw = gdept
1057	;
1058	;====================================================
1059	; e3[tw]:
1060	;====================================================
1061	;
1062	e3t = stepz
1063	IF n_elements(stepz) GT 1 THEN BEGIN
1064	e3w = 0.5*(stepz+shift(stepz, 1))
1065	e3w[0] = 0.5*e3t[0]
1066	ENDIF ELSE e3w = e3t
1067	;
1068	;====================================================
1069	; Mask
1070	;====================================================
1071	;
1072	; default mask eq 1
1073	if NOT keyword_set(mask) then tmask = -1 ELSE tmask = mask
1074	;
1075	if tmask[0] NE -1 then BEGIN
1076	tmask = byte(temporary(tmask))
1077	IF n_elements(mask) EQ nxx*nyy AND nzz GT 1 THEN BEGIN
1078	tmask = tmask[*]#replicate(1b, nzz)
1079	tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
1080	ENDIF
1081	IF nxx EQ 1 OR nyy EQ 1 OR nzz EQ 1 THEN tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
1082	tmask = tmask[xmin:xmax, ymin:ymax, zmin:zmax]
1083	IF jpi EQ 1 OR jpj EQ 1 OR jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1084	if key_shift NE 0 then tmask = shift(tmask, key_shift, 0, 0)
1085	; because tmask = reverse(tmask, 2) is not working if the 3rd
1086	; dimension of tmask = 1, we call reform.
1087	IF jpk EQ 1 THEN tmask = reform(tmask, /over)
1088	IF key_yreverse EQ 1 THEN tmask = reverse(tmask, 2)
1089	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1090	IF key_zreverse EQ 1 THEN tmask = reverse(tmask, 3)
1091	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1092	IF keyword_set(fullcgrid) THEN BEGIN
1093	IF keyword_set(key_periodic) THEN BEGIN
1094	msk = tmask*shift(tmask, -1, 0, 0)
1095	umaskred = msk[jpi-1, , ]
1096	ENDIF ELSE umaskred = tmask[jpi-1, , ]
1097	vmaskred = tmask[, jpj-1, ]
1098	ENDIF
1099	fmaskredy = tmask[jpi-1, , ]
1100	fmaskredx = tmask[, jpj-1, ]
1101	ENDIF ELSE BEGIN
1102	tmask = replicate(1b, jpi, jpj, jpk)
1103	IF keyword_set(fullcgrid) THEN BEGIN
1104	umaskred = replicate(1b, jpj, jpk)
1105	vmaskred = replicate(1b, jpi, jpk)
1106	ENDIF
1107	fmaskredy = replicate(1b, jpj, jpk)
1108	fmaskredx = replicate(1b, jpi, jpk)
1109	ENDELSE
1110	;
1111	IF jpi GT 2 AND jpj GT 2 AND NOT keyword_set(plain) $
1112	;; AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
1113	;; AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 $
1114	AND total(tmask[, 0, ]) EQ 0 AND total(tmask[, jpj-1, ]) EQ 0 $
1115	AND total(tmask[0, , ]) EQ 0 AND total(tmask[jpi-1, , ]) EQ 0 THEN BEGIN
1116	IF NOT keyword_set(key_periodic) THEN BEGIN
1117	xminmesh = ixminmesh + 1
1118	xmaxmesh = ixmaxmesh - 1
1119	ENDIF
1120	yminmesh = iyminmesh + 1
1121	ymaxmesh = iymaxmesh - 1
1122	; come back to the original grid before calling computegrid with the
1123	; new parameters...
1124	IF keyword_set(yreverse) THEN BEGIN
1125	gphit = reverse(gphit, 2)
1126	glamt = reverse(glamt, 2)
1127	tmask = reverse(tmask, 2)
1128	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1129	ENDIF
1130	IF keyword_set(zreverse) THEN BEGIN
1131	zaxis = reverse(zaxis)
1132	tmask = reverse(tmask, 3)
1133	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1134	ENDIF
1135	IF key_shift NE 0 THEN BEGIN
1136	glamt = shift(glamt, -key_shift, 0)
1137	gphit = shift(gphit, -key_shift, 0)
1138	tmask = shift(tmask, -key_shift, 0, 0)
1139	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1140	ENDIF
1141	IF ixminmesh NE 0 THEN BEGIN
1142	glamt = [fltarr(ixminmesh, jpj), glamt]
1143	gphit = [fltarr(ixminmesh, jpj), gphit]
1144	tmask = [fltarr(ixminmesh, jpj), tmask]
1145	jpi = jpi+ixminmesh
1146	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1147	ENDIF
1148	IF ixmaxmesh NE jpiglo-1 THEN BEGIN
1149	glamt = [glamt, fltarr(jpiglo-1-ixmaxmesh, jpj)]
1150	gphit = [gphit, fltarr(jpiglo-1-ixmaxmesh, jpj)]
1151	tmask = [tmask, fltarr(jpiglo-1-ixmaxmesh, jpj)]
1152	jpi = jpi+jpiglo-1-ixmaxmesh
1153	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1154	ENDIF
1155	IF iyminmesh NE 0 THEN BEGIN
1156	glamt = [[fltarr(jpi, iyminmesh)], [glamt]]
1157	gphit = [[fltarr(jpi, iyminmesh)], [gphit]]
1158	tmask = [[fltarr(jpi, iyminmesh)], [tmask]]
1159	jpj = jpj+iyminmesh
1160	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1161	ENDIF
1162	IF iymaxmesh NE jpjglo-1 THEN BEGIN
1163	glamt = [[glamt], [fltarr(jpi, jpjglo-1-iymaxmesh)]]
1164	gphit = [[gphit], [fltarr(jpi, jpjglo-1-iymaxmesh)]]
1165	tmask = [[tmask], [fltarr(jpi, jpjglo-1-iymaxmesh)]]
1166	jpj = jpj+jpjglo-1-iymaxmesh
1167	IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
1168	ENDIF
1169	computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $
1170	, MASK = tmask, GLAMBOUNDARY = glamboundary $
1171	, XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
1172	, YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
1173	, ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
1174	, ONEARTH = onearth, PERIODIC = key_periodic $
1175	, PLAIN = plain, SHIFT = key_shift, STRIDE = key_stride $
1176	, FULLCGRID = fullcgrid, XYINDEX = xyindex $
1177	, UBASE2TBASE = ubase2tbase, VBASE2TBASE = vbase2tbase $
1178	, FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
1179	, ROMSH = romsh, _extra = ex
1180	return
1181	ENDIF
1182	;
1183	IF NOT keyword_set(fullcgrid) THEN BEGIN
1184	umaskred = !values.f_nan
1185	vmaskred = !values.f_nan
1186	ENDIF
1187	;
1188	;====================================================
1189	; stride...
1190	;====================================================
1191	;
1192	IF total(key_stride) GT 3 THEN BEGIN
1193	IF key_shift NE 0 THEN BEGIN
1194	; for explanation, see header of read_ncdf_varget.pro
1195	jpiright = key_shift
1196	jpileft = jpi - key_shift - ( (key_stride[0]-1)-((key_shift-1) MOD key_stride[0]) )
1197	jpi = ((jpiright-1)/key_stride[0]+1) + ((jpileft-1)/key_stride[0]+1)
1198	ENDIF ELSE jpi = (jpi-1)/key_stride[0]+1
1199	jpj = (jpj-1)/key_stride[1]+1
1200	jpk = (jpk-1)/key_stride[2]+1
1201	;
1202	glamt = (temporary(glamt))[0::stride[0], 0::stride[1]]
1203	gphit = (temporary(gphit))[0::stride[0], 0::stride[1]]
1204	e1t = (temporary(e1t))[0::stride[0], 0::stride[1]]
1205	e2t = (temporary(e2t))[0::stride[0], 0::stride[1]]
1206	tmask = (temporary(tmask))[0::stride[0], 0::stride[1], 0:*:stride[2]]
1207	gdept = gdept[0:*:stride[2]]
1208	gdepw = gdepw[0:*:stride[2]]
1209	e3t = e3t[0:*:stride[2]]
1210	e3w = e3w[0:*:stride[2]]
1211	; we must recompute glamf and gphif...
1212	IF jpi GT 1 THEN BEGIN
1213	if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
1214	OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
1215	stepxf = (glamt + 720) MOD 360
1216	stepxf = shift(stepxf, -1, -1) - stepxf
1217	stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ]
1218	stepxf = min(abs(stepxf), dimension = 3)
1219	IF NOT keyword_set(key_periodic) THEN $
1220	stepxf[jpi-1, ] = stepxf[jpi-2, ]
1221	ENDIF ELSE BEGIN
1222	stepxf = shift(glamt, -1, -1) - glamt
1223	IF keyword_set(key_periodic) THEN $
1224	stepxf[jpi-1, ] = 360 + stepxf[jpi-1, ] $
1225	ELSE stepxf[jpi-1, ] = stepxf[jpi-2, ]
1226	ENDELSE
1227	IF jpj GT 1 THEN BEGIN
1228	stepxf[, jpj-1] = stepxf[, jpj-2]
1229	stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2]
1230	ENDIF
1231	glamf = glamt + 0.5 * stepxf
1232	ENDIF ELSE glamf = glamt + 0.5
1233	IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN
1234	IF NOT keyword_set(glamboundary) THEN BEGIN
1235	bigger = where(glamf GE min(glamt)+360)
1236	glamf[bigger] = glamf[bigger]-360.
1237	ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
1238	ENDIF
1239	IF jpj GT 1 THEN BEGIN
1240	; we must compute stepyf: y distance between T(i,j) T(i+1,j+1)
1241	stepyf = shift(gphit, -1, -1) - gphit
1242	stepyf[, jpj-1] = stepyf[, jpj-2]
1243	IF jpi GT 1 THEN BEGIN
1244	if NOT keyword_set(key_periodic) THEN $
1245	stepyf[jpi-1, ] = stepyf[jpi-2, ]
1246	stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
1247	ENDIF
1248	gphif = gphit + 0.5 * stepyf
1249	ENDIF ELSE gphif = gphit + 0.5
1250	;
1251	IF jpj EQ 1 THEN BEGIN
1252	glamt = reform(glamt, jpi, jpj, /over)
1253	gphit = reform(gphit, jpi, jpj, /over)
1254	glamf = reform(glamf, jpi, jpj, /over)
1255	gphif = reform(gphif, jpi, jpj, /over)
1256	e1t = reform(e1t, jpi, jpj, /over)
1257	e2t = reform(e2t, jpi, jpj, /over)
1258	ENDIF
1259	;
1260	IF keyword_set(fullcgrid) THEN BEGIN
1261	glamu = (temporary(glamu))[0::stride[0], 0::stride[1]]
1262	gphiu = (temporary(gphiu))[0::stride[0], 0::stride[1]]
1263	e1u = (temporary(e1u))[0::stride[0], 0::stride[1]]
1264	e2u = (temporary(e2u))[0::stride[0], 0::stride[1]]
1265	glamv = (temporary(glamv))[0::stride[0], 0::stride[1]]
1266	gphiv = (temporary(gphiv))[0::stride[0], 0::stride[1]]
1267	e1v = (temporary(e1v))[0::stride[0], 0::stride[1]]
1268	e2v = (temporary(e2v))[0::stride[0], 0::stride[1]]
1269	e1f = (temporary(e1f))[0::stride[0], 0::stride[1]]
1270	e2f = (temporary(e2f))[0::stride[0], 0::stride[1]]
1271	umaskred = (temporary(umaskred))[0, 0::stride[1], 0::stride[2]]
1272	vmaskred = (temporary(vmaskred))[0::stride[0], 0, 0::stride[2]]
1273	fmaskredy = (temporary(fmaskredy))[0, 0::stride[1], 0::stride[2]]
1274	fmaskredx = (temporary(fmaskredx))[0::stride[0], 0, 0::stride[2]]
1275	IF jpj EQ 1 THEN BEGIN
1276	glamu = reform(glamu, jpi, jpj, /over)
1277	gphiu = reform(gphiu, jpi, jpj, /over)
1278	e1u = reform(e1u, jpi, jpj, /over)
1279	e2u = reform(e2u, jpi, jpj, /over)
1280	glamv = reform(glamv, jpi, jpj, /over)
1281	gphiv = reform(gphiv, jpi, jpj, /over)
1282	e1v = reform(e1v, jpi, jpj, /over)
1283	e2v = reform(e2v, jpi, jpj, /over)
1284	e1f = reform(e1f, jpi, jpj, /over)
1285	e2f = reform(e2f, jpi, jpj, /over)
1286	ENDIF
1287	ENDIF
1288	ENDIF
1289	;
1290	;====================================================
1291	; apply all the grid parameters
1292	;====================================================
1293	;
1294	@updateold
1295	domdef
1296	;
1297	;====================================================
1298	; Triangulation
1299	;====================================================
1300	;
1301	IF total(tmask) EQ jpijpjjpk $
1302	AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
1303	ELSE BEGIN
1304	; are we using ORCA2 ?
1305	IF jpiglo EQ 182 AND jpi EQ 180 AND jpjglo EQ 149 AND jpj EQ 148 THEN $
1306	triangles_list = triangule() ELSE triangles_list = triangule(/keep_cont)
1307	ENDELSE
1308	;
1309	;====================================================
1310	; time axis (default definition)
1311	;====================================================
1312	;
1313	IF n_elements(time) EQ 0 OR n_elements(jpt) EQ 0 THEN BEGIN
1314	jpt = 1
1315	time = 0
1316	ENDIF
1317	;
1318	IF NOT keyword_set(key_forgetold) THEN BEGIN
1319	@updateold
1320	ENDIF
1321	;====================================================
1322	; grid parameters used by xxx
1323	;====================================================
1324	;
1325	IF NOT keyword_set(strcalling) THEN BEGIN
1326	IF n_elements(ccmeshparameters) EQ 0 THEN strcalling = 'computegrid' $
1327	ELSE strcalling = ccmeshparameters.filename
1328	ENDIF
1329	IF n_elements(glamt) GE 2 THEN BEGIN
1330	glaminfo = moment(glamt)
1331	IF finite(glaminfo[2]) EQ 0 THEN glaminfo = glaminfo[0:1]
1332	gphiinfo = moment(gphit)
1333	IF finite(gphiinfo[2]) EQ 0 THEN gphiinfo = gphiinfo[0:1]
1334	ENDIF ELSE BEGIN
1335	glaminfo = glamt
1336	gphiinfo = gphit
1337	ENDELSE
1338	IF keyword_set(romsh) THEN $
1339	romszinfos = {h:romsh[xmin:xmax, ymin:ymax], zeta:replicate(0., jpi, jpj), theta_s:-1, theta_b:-1, hc:-1} $
1340	ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}
1341
1342	ccmeshparameters = {filename:strcalling $
1343	, glaminfo:float(string(glaminfo, format = '(E11.4)')) $
1344	, gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $
1345	, jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $
1346	, jpi:jpi, jpj:jpj, jpk:jpk $
1347	, ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $
1348	, iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
1349	, izminmesh:izminmesh, izmaxmesh:izmaxmesh $
1350	, key_shift:key_shift, key_periodic:key_periodic $
1351	, key_stride:key_stride, key_gridtype:key_gridtype $
1352	, key_yreverse:key_yreverse, key_zreverse:key_zreverse $
1353	, key_partialstep:key_partialstep, key_onearth:key_onearth}
1354
1355	ccreadparameters = {funclec_name:'read_ncdf' $
1356	, jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $
1357	, ixmindta:ixmindta, ixmaxdta:ixmaxdta $
1358	, iymindta:iymindta, iymaxdta:iymaxdta $
1359	, izmindta:izmindta, izmaxdta:izmaxdta}
1360	;------------------------------------------------------------
1361	IF keyword_set(key_performance) EQ 1 THEN $
1362	print, 'time computegrid', systime(1)-time1
1363	;------------------------------------------------------------
1364	return
1365	end

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

Download in other formats:

Original Format