source: trunk/SRC/Computation/grad.pro @ 406

Last change on this file since 406 was 371, checked in by pinsard, 16 years ago

improvements of headers (alignments of IDL prompt in examples)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.2 KB
RevLine 
[2]1;+
[232]2;
[163]3; @file_comments
[314]4; compute the gradient of a variable located on Arakawa C-grid.
[2]5;
[163]6; @categories
[167]7; Calculation
[2]8;
[163]9; @param FIELD
[238]10; The field for which we want to compute the gradient. A 2D (xy),
[242]11; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by
[314]12; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array.
[242]13; Note that the dimension of the array must suit the domain dimension.
[2]14;
[167]15; @param DIREC {type=scalar string}
16; the gradient direction: 'x', 'y', 'z'
[2]17;
[314]18; @keyword MILLION {default=0}{type=scalar: 0 or 1}
19; Activate to multiply returned array by 1.e6
20;
21; @keyword _EXTRA
22; Used to declare that this routine accepts inherited keywords
23;
[238]24; @returns
[314]25; the gradient of the input data with the same size 2D, 3D or 4D
26; array, located on the appropriate grid (see restrictions)
[2]27;
[163]28; @uses
[371]29; <pro>cm_4cal</pro>
30; <pro>cm_4data</pro>
31; <pro>cm_4mesh</pro>
[2]32;
[163]33; @restrictions
[226]34; - Works only for Arakawa C-grid.
[167]35; - When computing the gradient, the result is not on the same grid point
36;   than the input data. In consequence, we update, vargrid and the grid position
37;   parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf],
38;   ny[tuvf], firstz[tw], lastz[tw], nz[tw]).
[226]39; - points that cannot be computed (domain boundaries, coastline) are set to NaN
[167]40; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
[142]41;
[163]42; @examples
43;
[371]44;   IDL> \@tst_initorca2
45;   IDL> plt, grad({arr:gphit,g:'T'}, 'x')
46;   IDL> plt, grad({arr:gphit,g:'T'}, 'y')
47;
[163]48; @history
49; Sebastien Masson (smasson\@lodyc.jussieu.fr)
50;
51; @version
52; $Id$
53;
[2]54;-
[327]55FUNCTION grad, field, direc, MILLION=million, _EXTRA=ex
[114]56;
57  compile_opt idl2, strictarrsubs
58;
[167]59@cm_4cal   ; for jpt
60@cm_4data  ; for varname, vargrid, vardate, varunit, valmask
[226]61@cm_4mesh
[314]62;
63  time1 = systime(1)          ; for key_performance
[2]64;------------------------------------------------------------
[25]65;
[167]66  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
[25]67     return, report(['This version of grad is based on Arakawa C-grid.' $
68                     , 'U and V grids must therefore be defined'])
69;
[167]70  res = litchamp(field)
71  szres = size(res)
72  grille, mask, glam, gphi, gdep, nx, ny, nz $
[226]73          , firstx, firsty, firstz, lastx, lasty, lastz
[167]74;
75  if n_elements(valmask) EQ 0 then valmask = 1.e20
76  varname = 'grad of '+varname
[314]77    IF keyword_set(million) THEN varunit = '1.e6*'+varunit+'/m' ELSE varunit = varunit+'/m'
[167]78  case strupcase(vargrid) of
79    'T':BEGIN
80      case direc of
[226]81        'x':BEGIN
[167]82          divi = e1u[firstx:lastx, firsty:lasty]
83          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
84          vargrid = 'U'
[226]85          firstxu = firstxt & lastxu = lastxt & nxu = nxt
[167]86          firstyu = firstyt & lastyu = lastyt & nyu = nyt
87        END
88        'y':BEGIN
89          divi = e2v[firstx:lastx, firsty:lasty]
90          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
91          vargrid = 'V'
[226]92          firstxv = firstxt & lastxv = lastxt & nxv = nxt
[167]93          firstyv = firstyt & lastyv = lastyt & nyv = nyt
94        END
95        'z':BEGIN
96          divi = e3w[firstz:lastz]
97          newmask = mask
98          vargrid = 'W'
[226]99          firstzw = firstzt & lastzw = lastzt & nzw = nzt
[167]100        END
101        ELSE:return, report('Bad definition of direction argument')
102      ENDCASE
103    END
104    'W':BEGIN
105      case direc of
[226]106        'x':BEGIN
[167]107          divi = e1u[firstx:lastx, firsty:lasty]
108          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
109          vargrid = 'U'
[226]110          firstxu = firstxt & lastxu = lastxt & nxu = nxt
[167]111          firstyu = firstyt & lastyu = lastyt & nyu = nyt
112        END
113        'y':BEGIN
114          divi = e2v[firstx:lastx, firsty:lasty]
115          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
116          vargrid = 'V'
[226]117          firstxv = firstxt & lastxv = lastxt & nxv = nxt
[167]118          firstyv = firstyt & lastyv = lastyt & nyv = nyt
119        END
120        'z':BEGIN
121          divi = e3t[firstz:lastz]
122          newmask = mask
123          vargrid = 'T'
[226]124          firstzt = firstzw & lastzt = lastzw & nzt = nzw
[167]125        END
126        ELSE:return, report('Bad definition of direction argument')
127      endcase
128    END
129    'U':BEGIN
130      case direc of
131        'x':BEGIN
132          divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
133          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
134          vargrid = 'T'
[226]135          firstxt = firstxu & lastxt = lastxu & nxt = nxu
[167]136          firstyt = firstyu & lastyt = lastyu & nyt = nyu
137        END
138        'y':BEGIN
139          divi = e2f[firstx:lastx, firsty:lasty]
140          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
141          vargrid = 'F'
[226]142          firstxf = firstxu & lastxf = lastxu & nxf = nxu
[167]143          firstyf = firstyu & lastyf = lastyu & nyf = nyu
144        END
145        'z':BEGIN
146          divi = e3w[firstz:lastz]
147          newmask = mask
148          vargrid = 'W'
[226]149          firstzw = firstzt & lastzw = lastzt & nzw = nzt
[167]150        END
151        ELSE:return, report('Bad definition of direction argument')
152      endcase
153    END
154    'V':BEGIN
155      case direc of
156        'x':BEGIN
157          divi = e1f[firstx:lastx, firsty:lasty]
158          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
159          vargrid = 'F'
160          firstxf = firstxv & lastxf = lastxv & nxf = nxv
161          firstyf = firstyv & lastyf = lastyv & nyf = nyv
162        END
163        'y':BEGIN
164          divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
165          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
166          vargrid = 'T'
167          firstxt = firstxv & lastxt = lastxv & nxt = nxv
168          firstyt = firstyv & lastyt = lastyv & nyt = nyv
169        END
170        'z':BEGIN
171          divi = e3w[firstz:lastz]
172          newmask = mask
173          vargrid = 'W'
[226]174          firstzw = firstzt & lastzw = lastzt & nzw = nzt
[167]175        END
176        ELSE:return, report('Bad definition of direction argument')
177      endcase
178    END
179    'F':BEGIN
[2]180;          case direc of
[25]181;             'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
182;             'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
183;             'z':divi = e3w[firstz:lastz]
[2]184;             ELSE:return, report('Bad definition of direction argument')
185;          endcase
[167]186      return, report('F grid: case not coded, please contact SAXO team...')
187    END
188    ELSE:return, report('Bad definition of vargrid')
189  ENDCASE
[314]190  IF keyword_set(million) THEN divi = temporary(divi)*1.e-6
[167]191  res = fitintobox(temporary(res))
192  IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
193  case 1 of
[2]194;----------------------------------------------------------------------------
195;xy
196;----------------------------------------------------------------------------
[167]197    szres[0] EQ 2:BEGIN
198      land = where((temporary(mask))[*, *, firstz] EQ 0)
199      if land[0] NE -1 then res[temporary(land)] = !values.f_nan
200      case direc of
[226]201        'x':BEGIN
[167]202          res = (shift(res, -1, 0)-res)/temporary(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(temporary(res), 1, 0)
205        END
[226]206        'y':BEGIN
[167]207          res = (shift(res, 0, -1)-res)/temporary(divi)
208          res[*, ny-1] = !values.f_nan
[226]209          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1)
[167]210        END
211        ELSE:return,  report('Bad definition of direction argument for the type of array')
212      ENDCASE
213      land = where((temporary(newmask))[*, *, firstz] EQ 0)
214      if land[0] NE -1 then res[temporary(land)] = valmask
215    END
[2]216;----------------------------------------------------------------------------
217;xyt
218;----------------------------------------------------------------------------
[167]219    szres[0] EQ 3 AND jpt NE 1:BEGIN
220      land = where((temporary(mask))[*, *, firstz] EQ 0, cnt)
221      if land[0] NE -1 then BEGIN
222        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
223        res[temporary(land)] = !values.f_nan
224      ENDIF
225      divi = (temporary(divi))[*]#replicate(1., jpt)
226      case direc of
[226]227        'x':BEGIN
[167]228          res = (shift(res, -1, 0, 0)-res)/temporary(divi)
229          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
230          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
231        END
[226]232        'y':BEGIN
[167]233          res = (shift(res, 0, -1, 0)-res)/temporary(divi)
234          res[*, ny-1, *] = !values.f_nan
[226]235          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
[167]236        END
237        ELSE:return,  report('Bad definition of direction argument for the type of array')
238      ENDCASE
239      land = where((temporary(newmask))[*, *, firstz] EQ 0, cnt)
240      if land[0] NE -1 then BEGIN
241        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
242        res[temporary(land)] = valmask
243      ENDIF
244    END
[2]245;----------------------------------------------------------------------------
246;xyz
247;----------------------------------------------------------------------------
[167]248    szres[0] EQ 3 AND jpt EQ 1:BEGIN
249      land = where(mask EQ 0)
250      if land[0] NE -1 then res[temporary(land)] = !values.f_nan
251      case direc OF
[226]252        'x':BEGIN
[167]253          divi = (temporary(divi))[*]#replicate(1., nz)
254          res = (shift(res, -1, 0, 0)-res)/temporary(divi)
255          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
256          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
257        END
[226]258        'y':BEGIN
[167]259          divi = (temporary(divi))[*]#replicate(1., nz)
260          res = (shift(res, 0, -1, 0)-res)/temporary(divi)
261          res[*, ny-1, *] = !values.f_nan
[226]262          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
[167]263        END
264        'z':BEGIN
265          divi = replicate(1., nx*ny)#(temporary(divi))[*]
266          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite)
[226]267          if vargrid EQ 'W' THEN BEGIN
[167]268            res = (shift(res, 0, 0, 1)-res)/temporary(divi)
269            res[*, *, 0] = !values.f_nan
270          ENDIF ELSE BEGIN
271            res = (res-shift(res, 0, 0, -1))/temporary(divi)
272            res[*, *, nz-1] = !values.f_nan
273          ENDELSE
274        END
275      ENDCASE
276      land = where(temporary(newmask) EQ 0)
277      if land[0] NE -1 then res[temporary(land)] = valmask
278    END
[2]279;----------------------------------------------------------------------------
280;xyzt
281;----------------------------------------------------------------------------
[167]282    szres[0] EQ 4:BEGIN
283      land = where(temporary(mask) EQ 0, cnt)
284      if land[0] NE -1 then BEGIN
285        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
286        res[temporary(land)] = !values.f_nan
287      ENDIF
288      case direc OF
[226]289        'x':BEGIN
[167]290          divi = (temporary(divi))[*]#replicate(1., nz*jpt)
291          res = (shift(res, -1, 0, 0, 0)-res)/temporary(divi)
292          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan
293          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0, 0)
294        END
[226]295        'y':BEGIN
[167]296          divi = (temporary(divi))[*]#replicate(1., nz*jpt)
297          res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi)
298          res[*, ny-1, *, *] = !values.f_nan
[226]299          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0)
[167]300        END
301        'z':BEGIN
302          divi = replicate(1., nx*ny)#(temporary(divi))[*]
303          divi = (temporary(divi))[*]#replicate(1L, jpt)
304          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite)
[226]305          if vargrid EQ 'W' THEN BEGIN
[167]306            res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi)
307            res[*, *, 0, *] = !values.f_nan
308          ENDIF ELSE BEGIN
309            res = (res-shift(res, 0, 0, -1, 0))/temporary(divi)
310            res[*, *, nz-1, *] = !values.f_nan
311          ENDELSE
312        END
313      ENDCASE
314      land = where(newmask EQ 0, cnt)
315      if land[0] NE -1 then BEGIN
316        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
317        res[temporary(land)] = valmask
318      ENDIF
319    END
320    ELSE:return, report('input array must have 2, 3 or 4 dimensions')
321  ENDCASE
[314]322
323  if keyword_set(key_performance) THEN print, 'time curl', systime(1)-time1
324
[167]325  return, res
326END
Note: See TracBrowser for help on using the repository browser.