Context Navigation

computegrid.pro @ 173

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

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

Download in other formats:

Original Format