source: trunk/grad.pro

Last change on this file was 48, checked in by pinsard, 10 years ago

fix thanks to coding rules
File size: 9.2 KB
Line 
1;+
2;
3; @history
4; Sebastien Masson (smasson\@lodyc.jussieu.fr)
5;
6;-
7FUNCTION grad, field, direc
8;
9@common
10;
11   res = litchamp(field)
12   taille=size(res)
13   grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz
14   if n_elements(valmask) EQ 0 then valmask = 1e20
15   case strupcase(vargrid) of
16      'T':BEGIN
17         case direc of
18            'x':BEGIN
19               divi = e1u[premierx:dernierx, premiery:derniery]
20               newmask = (umask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
21               vargrid = 'U'
22               domdef, glamt[premierx, 0], glamu[dernierx, 0] $
23                , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U']
24            END
25            'y':BEGIN
26               divi = e2v[premierx:dernierx, premiery:derniery]
27               newmask = (vmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
28               vargrid = 'V'
29               domdef, glamt[premierx, 0], glamv[dernierx, 0] $
30                , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V']
31            END
32            'z':BEGIN
33               divi = e3w[premierz:dernierz]
34               newmask = mask
35               vargrid = 'W'
36            END
37            ELSE : return, report('Bad definition of direction argument')
38         ENDCASE
39      END
40      'W':BEGIN
41         case direc of
42            'x':divi = e1u[premierx:dernierx, premiery:derniery]
43            'y':divi = e2v[premierx:dernierx, premiery:derniery]
44            'z':BEGIN
45               divi = e3t[premierz:dernierz]
46               newmask = mask
47               vargrid = 'T'
48            END
49            ELSE : return, report('Bad definition of direction argument')
50         endcase
51      END
52      'U':BEGIN
53         case direc of
54            'x':BEGIN
55               divi = (shift(e1t, -1, 0))[premierx:dernierx, premiery:derniery]
56               newmask = tmask[premierx:dernierx, premiery:derniery, premierz:dernierz]
57               vargrid = 'T'
58               domdef, glamt[premierx, 0], glamu[dernierx] $
59                , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U']
60            END
61            'y':BEGIN
62               divi = e2f[premierx:dernierx, premiery:derniery]
63               newmask = (fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
64               vargrid = 'F'
65               domdef, glamu[premierx, 0], glamf[dernierx, 0] $
66                , gphiu[0, premiery], gphif[0, derniery], grille = ['U','F']
67            END
68            'z':BEGIN
69               divi = e3w[premierz:dernierz]
70               newmask = mask
71               vargrid = 'W'
72            END
73            ELSE : return, report('Bad definition of direction argument')
74         endcase
75      END
76      'V':BEGIN
77         case direc of
78            'x':BEGIN
79               divi = e1f[premierx:dernierx, premiery:derniery]
80               newmask = (fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
81               vargrid = 'F'
82               domdef, glamv[premierx, 0], glamf[dernierx, 0] $
83                , gphiv[0, premiery], gphif[0, derniery], grille = ['V','F']
84            END
85            'y':BEGIN
86               divi = (shift(e2t, 0, -1))[premierx:dernierx, premiery:derniery]
87               newmask = tmask[premierx:dernierx, premiery:derniery, premierz:dernierz]
88               vargrid = 'T'
89               domdef, glamt[premierx, 0], glamv[dernierx, 0] $
90                , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V']
91            END
92            'z':BEGIN
93               divi = e3w[premierz:dernierz]
94               newmask = mask
95               vargrid = 'W'
96            END
97            ELSE : return, report('Bad definition of direction argument')
98         endcase
99      END
100;       'F':BEGIN
101;          case direc of
102;             'x':divi = (shift(e1v, -1, 0))[premierx:dernierx, premiery:derniery]
103;             'y':divi = (shift(e2u, 0, -1))[premierx:dernierx, premiery:derniery]
104;             'z':divi = e3w[premierz:dernierz]
105;             ELSE : return, report('Bad definition of direction argument')
106;          endcase
107;       END
108      ELSE : return, report('Bad definition of vargrid')
109   ENDCASE
110   res = fitintobox(res)
111   case 1 of
112;----------------------------------------------------------------------------
113;xy
114;----------------------------------------------------------------------------
115      taille[0] EQ 2:BEGIN
116         earth = where(mask[*, *, premierz] EQ 0)
117         if earth[0] NE -1 then res[earth] = !values.f_nan
118         case direc of
119            'x':BEGIN
120               res = (shift(res, -1, 0)-res)/divi
121               if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
122               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)
123            END
124            'y':BEGIN
125               res = (shift(res, 0, -1)-res)/divi
126               res[*, ny-1] = !values.f_nan
127               if vargrid EQ 'T' OR vargrid EQ 'U' then res =   shift(res, 0, 1)
128            END
129            ELSE : return,  report('Bad definition of direction argument for the type of array')
130         ENDCASE
131         earth = where(newmask[*, *, premierz] EQ 0)
132         if earth[0] NE -1 then res[earth] = valmask
133      END
134;----------------------------------------------------------------------------
135;xyt
136;----------------------------------------------------------------------------
137      taille[0] EQ 3 AND jpt NE 1:BEGIN
138         earth = where(mask[*, *, premierz] EQ 0)
139         if earth[0] NE -1 then BEGIN
140            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
141            res[earth[*]] = !values.f_nan
142         ENDIF
143         divi = divi[*]#replicate(1, jpt)
144         case direc of
145            'x':BEGIN
146               res = (shift(res, -1, 0, 0)-res)/divi
147               if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
148               if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0)
149            END
150            'y':BEGIN
151               res = (shift(res, 0, -1, 0)-res)/divi
152               res[*, ny-1, *] = !values.f_nan
153               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)
154            END
155            ELSE : return,  report('Bad definition of direction argument for the type of array')
156         ENDCASE
157         earth = where(newmask[*, *, premierz] EQ 0)
158         if earth[0] NE -1 then BEGIN
159            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
160            res[earth] = valmask
161         ENDIF
162      END
163;----------------------------------------------------------------------------
164;xyz
165;----------------------------------------------------------------------------
166      taille[0] EQ 3 AND jpt EQ 1:BEGIN
167         earth = where(mask EQ 0)
168         if earth[0] NE -1 then res[earth] = !values.f_nan
169         case direc OF
170            'x':BEGIN
171               divi = divi[*]#replicate(1, nz)
172               res = (shift(res, -1, 0, 0)-res)/divi
173               if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
174               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
175            END
176            'y':BEGIN
177               divi = divi[*]#replicate(1, nz)
178               res = (shift(res, 0, -1, 0)-res)/divi
179               res[*, ny-1, *] = !values.f_nan
180               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)
181            END
182            'z':BEGIN
183               divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)
184               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)
185               if vargrid EQ 'W' THEN BEGIN
186                  res = (shift(res, 0, 0, 1)-res)/divi
187                  res[*, *, 0] = !values.f_nan
188               ENDIF ELSE BEGIN
189                  res = (res-shift(res, 0, 0, -1))/divi
190                  res[*, *, nz-1] = !values.f_nan
191               ENDELSE
192               if earth[0] NE -1 then res[earth] = valmask
193            END
194         ENDCASE
195      END
196;------------------------------------------------------------
197;----------------------------------------------------------------------------
198;xyzt
199;----------------------------------------------------------------------------
200      taille[0] EQ 4:BEGIN
201         earth = where((mask[*]#replicate(1, jpt)) EQ 0)
202         if earth[0] NE -1 then res[earth] = !values.f_nan
203         case direc OF
204            'z':BEGIN
205               divi = replicate(1, nx*ny)#divi
206               divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over)
207               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt)
208               if vargrid EQ 'W' THEN BEGIN
209                  res = (shift(res, 0, 0, 1, 0)-res)/divi
210                  res[*, *, 0, *] = !values.f_nan
211               ENDIF ELSE BEGIN
212                  res = (res-shift(res, 0, 0, -1, 0))/divi
213                  res[*, *, nz-1, *] = !values.f_nan
214               ENDELSE
215               if earth[0] NE -1 then res[earth] = valmask
216            END
217         ENDCASE
218      END
219;------------------------------------------------------------
220;------------------------------------------------------------
221   endcase
222   varname = 'grad of '+varname
223   varunit = varunit+'/m'
224
225
226;------------------------------------------------------------
227   return, res
228end
Note: See TracBrowser for help on using the repository browser.