Context Navigation

grad.pro @ 159

Last change on this file since 159 was 159, checked in by smasson, 18 years ago
bugfix in grad.pro
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 11.1 KB

Line
1	;------------------------------------------------------------
2	;------------------------------------------------------------
3	;------------------------------------------------------------
4	;+
5	; NAME:
6	;
7	; PURPOSE:
8	;
9	; CATEGORY:
10	;
11	; CALLING SEQUENCE:
12	;
13	; INPUTS:
14	;
15	; KEYWORD PARAMETERS:
16	;
17	; OUTPUTS:
18	;
19	; COMMON BLOCKS:common.pro
20	;
21	; SIDE EFFECTS:
22	;
23	; RESTRICTIONS:
24	;
25	; EXAMPLE:
26	;
27	; MODIFICATION HISTORY:Sebastien Masson (smasson\@lodyc.jussieu.fr)
28	;
29	; @todo seb: remplir les truc!
30	;
31	;
32	;-
33	;------------------------------------------------------------
34	;------------------------------------------------------------
35	;------------------------------------------------------------
36	FUNCTION grad, field, direc
37	;
38	compile_opt idl2, strictarrsubs
39	;
40	@common
41	;------------------------------------------------------------
42	;
43	IF finite(glamu[0])finite(gphiu[0])finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
44	return, report(['This version of grad is based on Arakawa C-grid.' $
45	, 'U and V grids must therefore be defined'])
46	;
47	res = litchamp(field)
48	taille = size(res)
49	grille, mask, glam, gphi, gdep, nx, ny, nz $
50	, firstx, firsty, firstz, lastx, lasty, lastz
51	if n_elements(valmask) EQ 0 then valmask = 1e20
52	case strupcase(vargrid) of
53	'T':BEGIN
54	case direc of
55	'x':BEGIN
56	divi = e1u[firstx:lastx, firsty:lasty]
57	newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
58	vargrid = 'U'
59	firstxu = firstxt & lastxu = lastxt & nxu = nxt
60	firstyu = firstyt & lastyu = lastyt & nyu = nyt
61	END
62	'y':BEGIN
63	divi = e2v[firstx:lastx, firsty:lasty]
64	newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
65	vargrid = 'V'
66	firstxv = firstxt & lastxv = lastxt & nxv = nxt
67	firstyv = firstyt & lastyv = lastyt & nyv = nyt
68	END
69	'z':BEGIN
70	divi = e3w[firstz:lastz]
71	newmask = mask
72	vargrid = 'W'
73	firstzw = firstzt & lastzw = lastzt & nzw = nzt
74	END
75	ELSE:return, report('Bad definition of direction argument')
76	ENDCASE
77	END
78	'W':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 = e3t[firstz:lastz]
96	newmask = mask
97	vargrid = 'T'
98	firstzt = firstzw & lastzt = lastzw & nzt = nzw
99	END
100	ELSE:return, report('Bad definition of direction argument')
101	endcase
102	END
103	'U':BEGIN
104	case direc of
105	'x':BEGIN
106	divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
107	newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
108	vargrid = 'T'
109	firstxt = firstxu & lastxt = lastxu & nxt = nxu
110	firstyt = firstyu & lastyt = lastyu & nyt = nyu
111	END
112	'y':BEGIN
113	divi = e2f[firstx:lastx, firsty:lasty]
114	newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
115	vargrid = 'F'
116	firstxf = firstxu & lastxf = lastxu & nxf = nxu
117	firstyf = firstyu & lastyf = lastyu & nyf = nyu
118	END
119	'z':BEGIN
120	divi = e3w[firstz:lastz]
121	newmask = mask
122	vargrid = 'W'
123	firstzw = firstzt & lastzw = lastzt & nzw = nzt
124	END
125	ELSE:return, report('Bad definition of direction argument')
126	endcase
127	END
128	'V':BEGIN
129	case direc of
130	'x':BEGIN
131	divi = e1f[firstx:lastx, firsty:lasty]
132	newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
133	vargrid = 'F'
134	firstxf = firstxv & lastxf = lastxv & nxf = nxv
135	firstyf = firstyv & lastyf = lastyv & nyf = nyv
136	END
137	'y':BEGIN
138	divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
139	newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
140	vargrid = 'T'
141	firstxt = firstxv & lastxt = lastxv & nxt = nxv
142	firstyt = firstyv & lastyt = lastyv & nyt = nyv
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	; 'F':BEGIN
154	; case direc of
155	; 'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
156	; 'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
157	; 'z':divi = e3w[firstz:lastz]
158	; ELSE:return, report('Bad definition of direction argument')
159	; endcase
160	; END
161	ELSE:return, report('Bad definition of vargrid')
162	ENDCASE
163	res = fitintobox(res)
164	IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
165	case 1 of
166	;----------------------------------------------------------------------------
167	;xy
168	;----------------------------------------------------------------------------
169	taille[0] EQ 2:BEGIN
170	earth = where(mask[, , firstz] EQ 0)
171	if earth[0] NE -1 then res[earth] = !values.f_nan
172	case direc of
173	'x':BEGIN
174	res = (shift(res, -1, 0)-res)/divi
175	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
176	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)
177	END
178	'y':BEGIN
179	res = (shift(res, 0, -1)-res)/divi
180	res[*, ny-1] = !values.f_nan
181	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1)
182	END
183	ELSE:return, report('Bad definition of direction argument for the type of array')
184	ENDCASE
185	earth = where(newmask[, , firstz] EQ 0)
186	if earth[0] NE -1 then res[earth] = valmask
187	END
188	;----------------------------------------------------------------------------
189	;xyt
190	;----------------------------------------------------------------------------
191	taille[0] EQ 3 AND jpt NE 1:BEGIN
192	earth = where(mask[, , firstz] EQ 0)
193	if earth[0] NE -1 then BEGIN
194	earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nxnylindgen(jpt))
195	res[earth[*]] = !values.f_nan
196	ENDIF
197	divi = divi[*]#replicate(1, jpt)
198	case direc of
199	'x':BEGIN
200	res = (shift(res, -1, 0, 0)-res)/divi
201	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, , ] = !values.f_nan
202	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
203	END
204	'y':BEGIN
205	res = (shift(res, 0, -1, 0)-res)/divi
206	res[, ny-1, ] = !values.f_nan
207	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0)
208	END
209	ELSE:return, report('Bad definition of direction argument for the type of array')
210	ENDCASE
211	earth = where(newmask[, , firstz] EQ 0)
212	if earth[0] NE -1 THEN res[earth] = valmask
213	END
214	;----------------------------------------------------------------------------
215	;xyz
216	;----------------------------------------------------------------------------
217	taille[0] EQ 3 AND jpt EQ 1:BEGIN
218	earth = where(mask EQ 0)
219	if earth[0] NE -1 then res[earth] = !values.f_nan
220	case direc OF
221	'x':BEGIN
222	divi = divi[*]#replicate(1, nz)
223	res = (shift(res, -1, 0, 0)-res)/divi
224	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, , ] = !values.f_nan
225	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
226	END
227	'y':BEGIN
228	divi = divi[*]#replicate(1, nz)
229	res = (shift(res, 0, -1, 0)-res)/divi
230	res[, ny-1, ] = !values.f_nan
231	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0)
232	END
233	'z':BEGIN
234	divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)
235	if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)
236	if vargrid EQ 'W' THEN BEGIN
237	res = (shift(res, 0, 0, 1)-res)/divi
238	res[, , 0] = !values.f_nan
239	ENDIF ELSE BEGIN
240	res = (res-shift(res, 0, 0, -1))/divi
241	res[, , nz-1] = !values.f_nan
242	ENDELSE
243	if earth[0] NE -1 then res[earth] = valmask
244	END
245	ENDCASE
246	END
247	;------------------------------------------------------------
248	;----------------------------------------------------------------------------
249	;xyzt
250	;----------------------------------------------------------------------------
251	taille[0] EQ 4:BEGIN
252	earth = where((mask[*]#replicate(1, jpt)) EQ 0)
253	if earth[0] NE -1 then res[earth] = !values.f_nan
254	case direc OF
255	'x':BEGIN
256	divi = divi[]#replicate(1, nzjpt)
257	res = (shift(res, -1, 0, 0, 0)-res)/divi
258	if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, , , *] = !values.f_nan
259	if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0)
260	END
261	'y':BEGIN
262	divi = divi[]#replicate(1, nzjpt)
263	res = (shift(res, 0, -1, 0, 0)-res)/divi
264	res[, ny-1, , *] = !values.f_nan
265	if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0, 0)
266	END
267	'z':BEGIN
268	divi = replicate(1, nx*ny)#divi
269	divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over)
270	if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt)
271	if vargrid EQ 'W' THEN BEGIN
272	res = (shift(res, 0, 0, 1, 0)-res)/divi
273	res[, , 0, *] = !values.f_nan
274	ENDIF ELSE BEGIN
275	res = (res-shift(res, 0, 0, -1, 0))/divi
276	res[, , nz-1, *] = !values.f_nan
277	ENDELSE
278	END
279	ENDCASE
280	if earth[0] NE -1 then res[earth] = valmask
281	END
282	;------------------------------------------------------------
283	;------------------------------------------------------------
284	endcase
285	varname = 'grad of '+varname
286	varunit = varunit+'/m'
287
288
289	;------------------------------------------------------------
290	return, res
291	end
292
293
294
295
296
297

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

Download in other formats:

Original Format