source: trunk/grad.pro @ 2

Last change on this file since 2 was 2, checked in by pinsard, 18 years ago

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