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

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

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