Context Navigation

computegrid.pro @ 69

Last change on this file since 69 was 69, checked in by smasson, 18 years ago
debug + new xxx
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 38.9 KB

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

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

Download in other formats:

Original Format