Context Navigation

grad.pro @ 231

Last change on this file since 231 was 231, checked in by pinsard, 17 years ago
improvements/corrections of some *.pro headers
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 11.9 KB

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

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

Download in other formats:

Original Format