Context Navigation

grad.pro @ 168

Last change on this file since 168 was 168, checked in by pinsard, 18 years ago
Main document available on top directory, Source links available in idldoc html output
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 12.3 KB

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

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

Download in other formats:

Original Format