Context Navigation

grad.pro @ 325

Last change on this file since 325 was 325, checked in by pinsard, 17 years ago
modification of some headers (+some corrections) to prepare usage of the new idldoc
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 12.2 KB

Line
1	;+
2	;
3	; @file_comments
4	; compute the gradient of a variable located on Arakawa C-grid.
5	;
6	; @categories
7	; Calculation
8	;
9	; @param FIELD
10	; The field for which we want to compute the gradient. A 2D (xy),
11	; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by
12	; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array.
13	; Note that the dimension of the array must suit the domain dimension.
14	;
15	; @param DIREC {type=scalar string}
16	; the gradient direction: 'x', 'y', 'z'
17	;
18	; @keyword MILLION {default=0}{type=scalar: 0 or 1}
19	; Activate to multiply returned array by 1.e6
20	;
21	; @keyword _EXTRA
22	; Used to declare that this routine accepts inherited keywords
23	;
24	; @returns
25	; the gradient of the input data with the same size 2D, 3D or 4D
26	; array, located on the appropriate grid (see restrictions)
27	;
28	; @uses
29	; cm_4cal
30	; cm_4data
31	; cm_4mesh
32	;
33	; @restrictions
34	; - Works only for Arakawa C-grid.
35	; - When computing the gradient, the result is not on the same grid point
36	; than the input data. In consequence, we update, vargrid and the grid position
37	; parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf],
38	; ny[tuvf], firstz[tw], lastz[tw], nz[tw]).
39	; - points that cannot be computed (domain boundaries, coastline) are set to NaN
40	; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
41	;
42	; @examples
43	; IDL> \@tst_initorca2
44	; IDL> plt, grad({arr:gphit,g:'T'}, 'x')
45	; IDL> plt, grad({arr:gphit,g:'T'}, 'y')
46	;
47	; @history
48	; Sebastien Masson (smasson\@lodyc.jussieu.fr)
49	;
50	; @version
51	; $Id$
52	;
53	;-
54	FUNCTION grad, field, direc, MILLION = million, _EXTRA = ex
55	;
56	compile_opt idl2, strictarrsubs
57	;
58	@cm_4cal ; for jpt
59	@cm_4data ; for varname, vargrid, vardate, varunit, valmask
60	@cm_4mesh
61	;
62	time1 = systime(1) ; for key_performance
63	;------------------------------------------------------------
64	;
65	IF finite(glamu[0])finite(gphiu[0])finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
66	return, report(['This version of grad is based on Arakawa C-grid.' $
67	, 'U and V grids must therefore be defined'])
68	;
69	res = litchamp(field)
70	szres = size(res)
71	grille, mask, glam, gphi, gdep, nx, ny, nz $
72	, firstx, firsty, firstz, lastx, lasty, lastz
73	;
74	if n_elements(valmask) EQ 0 then valmask = 1.e20
75	varname = 'grad of '+varname
76	IF keyword_set(million) THEN varunit = '1.e6*'+varunit+'/m' ELSE varunit = varunit+'/m'
77	case strupcase(vargrid) of
78	'T':BEGIN
79	case direc of
80	'x':BEGIN
81	divi = e1u[firstx:lastx, firsty:lasty]
82	newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
83	vargrid = 'U'
84	firstxu = firstxt & lastxu = lastxt & nxu = nxt
85	firstyu = firstyt & lastyu = lastyt & nyu = nyt
86	END
87	'y':BEGIN
88	divi = e2v[firstx:lastx, firsty:lasty]
89	newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
90	vargrid = 'V'
91	firstxv = firstxt & lastxv = lastxt & nxv = nxt
92	firstyv = firstyt & lastyv = lastyt & nyv = nyt
93	END
94	'z':BEGIN
95	divi = e3w[firstz:lastz]
96	newmask = mask
97	vargrid = 'W'
98	firstzw = firstzt & lastzw = lastzt & nzw = nzt
99	END
100	ELSE:return, report('Bad definition of direction argument')
101	ENDCASE
102	END
103	'W':BEGIN
104	case direc of
105	'x':BEGIN
106	divi = e1u[firstx:lastx, firsty:lasty]
107	newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
108	vargrid = 'U'
109	firstxu = firstxt & lastxu = lastxt & nxu = nxt
110	firstyu = firstyt & lastyu = lastyt & nyu = nyt
111	END
112	'y':BEGIN
113	divi = e2v[firstx:lastx, firsty:lasty]
114	newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
115	vargrid = 'V'
116	firstxv = firstxt & lastxv = lastxt & nxv = nxt
117	firstyv = firstyt & lastyv = lastyt & nyv = nyt
118	END
119	'z':BEGIN
120	divi = e3t[firstz:lastz]
121	newmask = mask
122	vargrid = 'T'
123	firstzt = firstzw & lastzt = lastzw & nzt = nzw
124	END
125	ELSE:return, report('Bad definition of direction argument')
126	endcase
127	END
128	'U':BEGIN
129	case direc of
130	'x':BEGIN
131	divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
132	newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
133	vargrid = 'T'
134	firstxt = firstxu & lastxt = lastxu & nxt = nxu
135	firstyt = firstyu & lastyt = lastyu & nyt = nyu
136	END
137	'y':BEGIN
138	divi = e2f[firstx:lastx, firsty:lasty]
139	newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
140	vargrid = 'F'
141	firstxf = firstxu & lastxf = lastxu & nxf = nxu
142	firstyf = firstyu & lastyf = lastyu & nyf = nyu
143	END
144	'z':BEGIN
145	divi = e3w[firstz:lastz]
146	newmask = mask
147	vargrid = 'W'
148	firstzw = firstzt & lastzw = lastzt & nzw = nzt
149	END
150	ELSE:return, report('Bad definition of direction argument')
151	endcase
152	END
153	'V':BEGIN
154	case direc of
155	'x':BEGIN
156	divi = e1f[firstx:lastx, firsty:lasty]
157	newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
158	vargrid = 'F'
159	firstxf = firstxv & lastxf = lastxv & nxf = nxv
160	firstyf = firstyv & lastyf = lastyv & nyf = nyv
161	END
162	'y':BEGIN
163	divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
164	newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
165	vargrid = 'T'
166	firstxt = firstxv & lastxt = lastxv & nxt = nxv
167	firstyt = firstyv & lastyt = lastyv & nyt = nyv
168	END
169	'z':BEGIN
170	divi = e3w[firstz:lastz]
171	newmask = mask
172	vargrid = 'W'
173	firstzw = firstzt & lastzw = lastzt & nzw = nzt
174	END
175	ELSE:return, report('Bad definition of direction argument')
176	endcase
177	END
178	'F':BEGIN
179	; case direc of
180	; 'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
181	; 'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
182	; 'z':divi = e3w[firstz:lastz]
183	; ELSE:return, report('Bad definition of direction argument')
184	; endcase
185	return, report('F grid: case not coded, please contact SAXO team...')
186	END
187	ELSE:return, report('Bad definition of vargrid')
188	ENDCASE
189	IF keyword_set(million) THEN divi = temporary(divi)*1.e-6
190	res = fitintobox(temporary(res))
191	IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
192	case 1 of
193	;----------------------------------------------------------------------------
194	;xy
195	;----------------------------------------------------------------------------
196	szres[0] EQ 2:BEGIN
197	land = where((temporary(mask))[, , firstz] EQ 0)
198	if land[0] NE -1 then res[temporary(land)] = !values.f_nan
199	case direc of
200	'x':BEGIN
201	res = (shift(res, -1, 0)-res)/temporary(divi)
202	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
203	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0)
204	END
205	'y':BEGIN
206	res = (shift(res, 0, -1)-res)/temporary(divi)
207	res[*, ny-1] = !values.f_nan
208	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1)
209	END
210	ELSE:return, report('Bad definition of direction argument for the type of array')
211	ENDCASE
212	land = where((temporary(newmask))[, , firstz] EQ 0)
213	if land[0] NE -1 then res[temporary(land)] = valmask
214	END
215	;----------------------------------------------------------------------------
216	;xyt
217	;----------------------------------------------------------------------------
218	szres[0] EQ 3 AND jpt NE 1:BEGIN
219	land = where((temporary(mask))[, , firstz] EQ 0, cnt)
220	if land[0] NE -1 then BEGIN
221	land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nxnylindgen(jpt))
222	res[temporary(land)] = !values.f_nan
223	ENDIF
224	divi = (temporary(divi))[*]#replicate(1., jpt)
225	case direc of
226	'x':BEGIN
227	res = (shift(res, -1, 0, 0)-res)/temporary(divi)
228	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, , ] = !values.f_nan
229	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
230	END
231	'y':BEGIN
232	res = (shift(res, 0, -1, 0)-res)/temporary(divi)
233	res[, ny-1, ] = !values.f_nan
234	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
235	END
236	ELSE:return, report('Bad definition of direction argument for the type of array')
237	ENDCASE
238	land = where((temporary(newmask))[, , firstz] EQ 0, cnt)
239	if land[0] NE -1 then BEGIN
240	land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nxnylindgen(jpt))
241	res[temporary(land)] = valmask
242	ENDIF
243	END
244	;----------------------------------------------------------------------------
245	;xyz
246	;----------------------------------------------------------------------------
247	szres[0] EQ 3 AND jpt EQ 1:BEGIN
248	land = where(mask EQ 0)
249	if land[0] NE -1 then res[temporary(land)] = !values.f_nan
250	case direc OF
251	'x':BEGIN
252	divi = (temporary(divi))[*]#replicate(1., nz)
253	res = (shift(res, -1, 0, 0)-res)/temporary(divi)
254	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, , ] = !values.f_nan
255	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
256	END
257	'y':BEGIN
258	divi = (temporary(divi))[*]#replicate(1., nz)
259	res = (shift(res, 0, -1, 0)-res)/temporary(divi)
260	res[, ny-1, ] = !values.f_nan
261	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
262	END
263	'z':BEGIN
264	divi = replicate(1., nxny)#(temporary(divi))[]
265	if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite)
266	if vargrid EQ 'W' THEN BEGIN
267	res = (shift(res, 0, 0, 1)-res)/temporary(divi)
268	res[, , 0] = !values.f_nan
269	ENDIF ELSE BEGIN
270	res = (res-shift(res, 0, 0, -1))/temporary(divi)
271	res[, , nz-1] = !values.f_nan
272	ENDELSE
273	END
274	ENDCASE
275	land = where(temporary(newmask) EQ 0)
276	if land[0] NE -1 then res[temporary(land)] = valmask
277	END
278	;----------------------------------------------------------------------------
279	;xyzt
280	;----------------------------------------------------------------------------
281	szres[0] EQ 4:BEGIN
282	land = where(temporary(mask) EQ 0, cnt)
283	if land[0] NE -1 then BEGIN
284	land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nxnynz*lindgen(jpt))
285	res[temporary(land)] = !values.f_nan
286	ENDIF
287	case direc OF
288	'x':BEGIN
289	divi = (temporary(divi))[]#replicate(1., nzjpt)
290	res = (shift(res, -1, 0, 0, 0)-res)/temporary(divi)
291	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, , , *] = !values.f_nan
292	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0, 0)
293	END
294	'y':BEGIN
295	divi = (temporary(divi))[]#replicate(1., nzjpt)
296	res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi)
297	res[, ny-1, , *] = !values.f_nan
298	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0)
299	END
300	'z':BEGIN
301	divi = replicate(1., nxny)#(temporary(divi))[]
302	divi = (temporary(divi))[*]#replicate(1L, jpt)
303	if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite)
304	if vargrid EQ 'W' THEN BEGIN
305	res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi)
306	res[, , 0, *] = !values.f_nan
307	ENDIF ELSE BEGIN
308	res = (res-shift(res, 0, 0, -1, 0))/temporary(divi)
309	res[, , nz-1, *] = !values.f_nan
310	ENDELSE
311	END
312	ENDCASE
313	land = where(newmask EQ 0, cnt)
314	if land[0] NE -1 then BEGIN
315	land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nxnynz*lindgen(jpt))
316	res[temporary(land)] = valmask
317	ENDIF
318	END
319	ELSE:return, report('input array must have 2, 3 or 4 dimensions')
320	ENDCASE
321
322	if keyword_set(key_performance) THEN print, 'time curl', systime(1)-time1
323
324	return, res
325	END

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

Download in other formats:

Original Format