Context Navigation

computegrid.pro @ 210

Last change on this file since 210 was 209, checked in by smasson, 17 years ago
bugfix + introduce C grid based on F, U and V points
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 46.4 KB

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

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

Download in other formats:

Original Format