source: trunk/SRC/ToBeReviewed/CALCULS/grad.pro @ 134

Last change on this file since 134 was 134, checked in by navarro, 18 years ago

change *.pro file properties (del eof-style, del executable, set keywords Id
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.3 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;-
30;------------------------------------------------------------
31;------------------------------------------------------------
32;------------------------------------------------------------
33FUNCTION grad, field, direc
34;
35  compile_opt idl2, strictarrsubs
36;
37@common
38;------------------------------------------------------------
39;
40   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
41     return, report(['This version of grad is based on Arakawa C-grid.' $
42                     , 'U and V grids must therefore be defined'])
43;
44   res = litchamp(field)
45   taille=size(res)
46   grille, mask, glam, gphi, gdep, nx, ny,nz,firstx,firsty,firstz,lastx, lasty, lastz     
47   if n_elements(valmask) EQ 0 then valmask = 1e20
48   case strupcase(vargrid) of
49      'T':BEGIN
50         case direc of
51            'x':BEGIN
52               divi = e1u[firstx:lastx, firsty:lasty]
53               newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
54               vargrid = 'U'
55               domdef, glamt[firstx, 0], glamu[lastx, 0] $
56                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U']
57            END
58            'y':BEGIN
59               divi = e2v[firstx:lastx, firsty:lasty]
60               newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
61               vargrid = 'V'
62               domdef, glamt[firstx, 0], glamv[lastx, 0] $
63                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V']
64            END
65            'z':BEGIN
66               divi = e3w[firstz:lastz]
67               newmask = mask
68               vargrid = 'W'
69            END
70            ELSE:return, report('Bad definition of direction argument')
71         ENDCASE
72      END
73      'W':BEGIN
74         case direc of
75            'x':divi = e1u[firstx:lastx, firsty:lasty]
76            'y':divi = e2v[firstx:lastx, firsty:lasty]
77            'z':BEGIN
78               divi = e3t[firstz:lastz]
79               newmask = mask
80               vargrid = 'T'
81            END
82            ELSE:return, report('Bad definition of direction argument')
83         endcase
84      END
85      'U':BEGIN
86         case direc of
87            'x':BEGIN
88               divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
89               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
90               vargrid = 'T'
91               domdef, glamt[firstx, 0], glamu[lastx] $
92                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U']
93            END
94            'y':BEGIN
95               divi = e2f[firstx:lastx, firsty:lasty]
96               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
97               vargrid = 'F'
98               domdef, glamu[firstx, 0], glamf[lastx, 0] $
99                , gphiu[0, firsty], gphif[0, lasty], gridtype = ['U','F']
100            END
101            'z':BEGIN
102               divi = e3w[firstz:lastz]
103               newmask = mask
104               vargrid = 'W'
105            END
106            ELSE:return, report('Bad definition of direction argument')
107         endcase
108      END
109      'V':BEGIN
110         case direc of
111            'x':BEGIN
112               divi = e1f[firstx:lastx, firsty:lasty]
113               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
114               vargrid = 'F'
115               domdef, glamv[firstx, 0], glamf[lastx, 0] $
116                , gphiv[0, firsty], gphif[0, lasty], gridtype = ['V','F']
117            END
118            'y':BEGIN
119               divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
120               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
121               vargrid = 'T'
122               domdef, glamt[firstx, 0], glamv[lastx, 0] $
123                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V']
124            END
125            'z':BEGIN
126               divi = e3w[firstz:lastz]
127               newmask = mask
128               vargrid = 'W'
129            END
130            ELSE:return, report('Bad definition of direction argument')
131         endcase
132      END
133;       'F':BEGIN
134;          case direc of
135;             'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
136;             'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
137;             'z':divi = e3w[firstz:lastz]
138;             ELSE:return, report('Bad definition of direction argument')
139;          endcase
140;       END
141      ELSE:return, report('Bad definition of vargrid')
142   ENDCASE
143   res = fitintobox(res)
144   case 1 of
145;----------------------------------------------------------------------------
146;xy
147;----------------------------------------------------------------------------
148      taille[0] EQ 2:BEGIN
149         earth = where(mask[*, *, firstz] EQ 0)
150         if earth[0] NE -1 then res[earth] = !values.f_nan
151         case direc of
152            'x':BEGIN
153               res = (shift(res, -1, 0)-res)/divi
154               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
155               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)
156            END
157            'y':BEGIN
158               res = (shift(res, 0, -1)-res)/divi
159               res[*, ny-1] = !values.f_nan
160               if vargrid EQ 'T' OR vargrid EQ 'U' then res =   shift(res, 0, 1)           
161            END
162            ELSE:return,  report('Bad definition of direction argument for the type of array')
163         ENDCASE
164         earth = where(newmask[*, *, firstz] EQ 0)
165         if earth[0] NE -1 then res[earth] = valmask
166      END
167;----------------------------------------------------------------------------
168;xyt
169;----------------------------------------------------------------------------
170      taille[0] EQ 3 AND jpt NE 1:BEGIN
171         earth = where(mask[*, *, firstz] EQ 0)
172         if earth[0] NE -1 then BEGIN
173            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
174            res[earth[*]] = !values.f_nan
175         ENDIF
176         divi = divi[*]#replicate(1, jpt)
177         case direc of
178            'x':BEGIN
179               res = (shift(res, -1, 0, 0)-res)/divi
180               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
181               if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0)
182            END
183            'y':BEGIN
184               res = (shift(res, 0, -1, 0)-res)/divi
185               res[*, ny-1, *] = !values.f_nan
186               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)       
187            END
188            ELSE:return,  report('Bad definition of direction argument for the type of array')
189         ENDCASE
190         earth = where(newmask[*, *, firstz] EQ 0)
191         if earth[0] NE -1 THEN res[earth] = valmask
192      END
193;----------------------------------------------------------------------------
194;xyz
195;----------------------------------------------------------------------------
196      taille[0] EQ 3 AND jpt EQ 1:BEGIN
197         earth = where(mask EQ 0)
198         if earth[0] NE -1 then res[earth] = !values.f_nan
199         case direc OF
200            'x':BEGIN
201               divi = divi[*]#replicate(1, nz)
202               res = (shift(res, -1, 0, 0)-res)/divi
203               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
204               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
205            END
206            'y':BEGIN
207               divi = divi[*]#replicate(1, nz)
208               res = (shift(res, 0, -1, 0)-res)/divi
209               res[*, ny-1, *] = !values.f_nan
210               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)       
211            END
212            'z':BEGIN
213               divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)
214               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)
215               if vargrid EQ 'W' THEN BEGIN
216                  res = (shift(res, 0, 0, 1)-res)/divi
217                  res[*, *, 0] = !values.f_nan
218               ENDIF ELSE BEGIN
219                  res = (res-shift(res, 0, 0, -1))/divi
220                  res[*, *, nz-1] = !values.f_nan
221               ENDELSE
222               if earth[0] NE -1 then res[earth] = valmask
223            END
224         ENDCASE
225      END
226;------------------------------------------------------------
227;----------------------------------------------------------------------------
228;xyzt
229;----------------------------------------------------------------------------
230      taille[0] EQ 4:BEGIN
231         earth = where((mask[*]#replicate(1, jpt)) EQ 0)
232         if earth[0] NE -1 then res[earth] = !values.f_nan
233         case direc OF
234            'x':BEGIN
235               divi = divi[*]#replicate(1, nz*jpt)
236               res = (shift(res, -1, 0, 0, 0)-res)/divi
237               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan
238               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0)
239            END
240            'y':BEGIN
241               divi = divi[*]#replicate(1, nz*jpt)
242               res = (shift(res, 0, -1, 0, 0)-res)/divi
243               res[*, ny-1, *, *] = !values.f_nan
244               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0, 0)       
245            END
246            'z':BEGIN
247               divi = replicate(1, nx*ny)#divi
248               divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over)
249               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt)
250               if vargrid EQ 'W' THEN BEGIN
251                  res = (shift(res, 0, 0, 1, 0)-res)/divi
252                  res[*, *, 0, *] = !values.f_nan
253               ENDIF ELSE BEGIN
254                  res = (res-shift(res, 0, 0, -1, 0))/divi
255                  res[*, *, nz-1, *] = !values.f_nan
256               ENDELSE
257            END
258         ENDCASE
259         if earth[0] NE -1 then res[earth] = valmask
260      END
261;------------------------------------------------------------
262;------------------------------------------------------------
263   endcase
264   varname = 'grad of '+varname
265   varunit = varunit+'/m'
266
267
268;------------------------------------------------------------
269   return, res
270end
271
272
273
274
275
276
Note: See TracBrowser for help on using the repository browser.