Context Navigation

computegrid.pro @ 207

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

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

Download in other formats:

Original Format