source: trunk/SRC/ToBeReviewed/CALCULS/moyenne.pro @ 163

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

header improvements : type of parameters and keywords, default values, spell checking + idldoc assistant (IDL online_help)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 23.6 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5;
6; @file_comments 
7; averages a 2- or 3-d field over a selected
8;           geographical area and along one ore more
9;           selected axes (x, y or z)
10;
11; @categories
12;
13; @param TAB {in}{required}
14; 2 or 3d field
15;
16; @param DIREC {in}{required}
17; 'x' 'y' 'z' 'xy' 'xz' 'yz' or 'xyz'
18;
19; @keyword BOXZOOM 
20; [xmin,xmax,ymin,ymax (,(zmin,)zmax)] to more details, see domdef
21; boxzoom can have 5 forms:
22;  [vert2],
23;  [vert1, vert2],
24;  [lon1, lon2, lat1, lat2], 
25;  [lon1, lon2, lat1, lat2, vert2],
26;  [lon1, lon2, lat1, lat2, vert1,vert2]
27;
28; @keyword NAN
29; not a number, we activate it if we want to average without considerate some masked values of TAB.
30; If masked values of TAB are values consecrated by IDL(!values.f_nan), we just have to put NAN.
31; If masked values of TAB are valued a (a must be different of 1, corresponding to nan =
32; !values.f_nan and of 0, which desactivate nan). We have to put NAN=a.
33; Comment: In output, result points which are NAN will be valued a or !values.f_nan.
34;           
35; @keyword NODOMDEF
36; We activate it if we do not want to pass in domdef even if the keyword boxzoom
37; is present (like when grossemoyenne is called via checkfield)
38;
39; @keyword INTEGRATION
40; To make an integral rather than an average
41;
42; @keyword WDEPTH
43; to specify that the field is at W depth instead of T
44; depth (automatically activated if vargrid eq 'W')
45;           
46; @returns
47; An array
48;
49; @uses
50; common
51; domdef
52;
53; @restrictions Put values corresponding to land at 1.e20
54;
55; @history
56;  Jerome Vialard (jv\@lodyc.jussieu.fr)
57;                       2/7/98
58;                       Sebastien Masson (smasson\@lodyc.jussieu.fr)
59;                       14/8/98
60;                       15/1/98
61;                       11/3/99 adaptation for NAN
62;                       28/7/99 Averages 1d arrays
63;
64; @version
65; $Id$
66;
67;-
68
69function moyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $
70                  , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $
71                  , _extra = ex
72;---------------------------------------------------------
73;
74  compile_opt idl2, strictarrsubs
75;
76@cm_4mesh
77@cm_4data
78@cm_4cal
79  IF NOT keyword_set(key_forgetold) THEN BEGIN
80@updatenew
81@updatekwd
82  ENDIF
83;---------------------------------------------------------
84  tempsun = systime(1)          ; To key_performance
85;------------------------------------------------------------
86;I) preliminaries
87;------------------------------------------------------------
88  dirt = 0
89  dirx = 0
90  diry = 0
91  dirz = 0
92;------------------------------------------------------------
93; I.1) Directions(s) we follow to integrate
94;------------------------------------------------------------
95  if ( strpos(direc, 't') ge 0 ) then dirt = 1
96  if ( strpos(direc, 'x') ge 0 ) then dirx = 1
97  if ( strpos(direc, 'y') ge 0 ) then diry = 1
98  if ( strpos(direc, 'z') ge 0 ) then dirz = 1
99  if (dirx eq 0 and diry eq 0 and dirz eq 0) then return, tab
100;------------------------------------------------------------
101; I.2) verification of the input array's size
102;------------------------------------------------------------
103  taille = size(tab)
104  case 1 of
105    taille[0] eq 1 :dim = '1d'
106    taille[0] eq 2 :BEGIN
107      dim = '2d'
108      if dirx eq 0 and diry eq 0 then return, tab
109    END
110    taille[0] eq 3 :BEGIN
111      dim = '3d'
112      if dirx eq 0 and diry eq 0 and dirz eq 0 then return, tab
113    END
114    else : return, report('Le tableau d''entree doit etre a 2 ou 3 dimensions s''il contient une dim temporelle utiliser grossemoyenne')
115  endcase
116;------------------------------------------------------------
117;   I.3) Obtainment of scale's factors and of the mask on the subdomain concernedby the average.
118; Redefinition of the domain ajusted at boxzoom (at 6 elements)
119; This will allowed us to calculate only in the domain concerned by the average.
120; Domdef, followed by grid give us all arrays of the grid on the subdomain
121;------------------------------------------------------------
122  if keyword_set(boxzoom) then BEGIN
123    Case 1 Of
124      N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
125      N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
126      N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
127      N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
128      N_Elements(Boxzoom) Eq 6:bte = Boxzoom
129      Else: return, report('Mauvaise Definition de Boxzoom')
130    endcase
131    if NOT keyword_set(nodomdef) then BEGIN
132      savedbox = 1b
133      saveboxparam, 'boxparam4moyenne.dat'
134      domdef, bte, GRIDTYPE = vargrid, _extra = ex
135    ENDIF
136  ENDIF
137;---------------------------------------------------------------
138; attribution of the mask and of longitude and latitude arrays...
139;---------------------------------------------------------------
140  IF vargrid EQ 'W' THEN wdepth = 1
141  grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth
142;------------------------------------------------------------
143;------------------------------------------------------------
144; II) Case of the 1d array
145;------------------------------------------------------------
146;------------------------------------------------------------
147  if dim EQ '1d' then BEGIN
148    if n_elements(tab) NE nx*ny AND n_elements(tab) NE nx*ny*nz then BEGIN
149      if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
150      return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
151    ENDIF
152    case 1 of
153      nx EQ 1 AND ny EQ 1:BEGIN ;vector following z
154        case n_elements(tab) of
155          jpk:res = tab[firstz:lastz]
156          nz:res = tab
157          ELSE:BEGIN
158            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
159            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
160          END
161        ENDCASE
162        if dirz EQ 1 then BEGIN
163          dim = '3d'
164          taille = size(reform(res, nx, ny, nz))
165        ENDIF ELSE BEGIN
166          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
167          return, res
168        ENDELSE
169      END
170      ny EQ 1:BEGIN             ;vector following x
171        case n_elements(tab) of
172          jpi:res = tab[firstx:lastx]
173          nx:res = tab
174          ELSE:BEGIN
175            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
176            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
177          END
178        ENDCASE
179        if dirx EQ 1 then BEGIN
180          dim = '2d'
181          taille = size(reform(res, nx, ny))
182        ENDIF ELSE BEGIN
183          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
184          return, res
185        ENDELSE
186      END
187      nx EQ 1:BEGIN             ;vector following y
188        case n_elements(tab) of
189          jpj:res = tab[firsty:lasty]
190          ny:res = tab
191          ELSE:BEGIN
192            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
193            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
194          END
195        ENDCASE
196        if diry EQ 1 then BEGIN
197          dim = '2d'
198          taille = size(reform(res, nx, ny))
199        ENDIF ELSE BEGIN
200          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
201          return, res
202        ENDELSE
203      END
204    endcase
205  endif
206;------------------------------------------------------------
207;------------------------------------------------------------
208; II) Case of the 2d array
209;------------------------------------------------------------
210;------------------------------------------------------------
211  if (dim eq '2d') then begin
212;---------------------------------------------------------------
213;   II.1) verification of the coherence of the array's size to average
214; Verification of the coherence between the array's size and the domain defined by domdef
215; The input array must have either the total domain's size (jpi,jpj) or this
216; one of the reduced domain (nx,ny)
217;---------------------------------------------------------------
218    case 1 of
219      taille[1] eq jpi and taille[2] eq jpj: $
220        res = tab[firstx:lastx, firsty:lasty]
221      taille[1] eq  nx and taille[2] eq  ny:res = tab
222      else:BEGIN
223        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
224        return, report('Probleme d''adequation entre les tailles du domaine nx*ny '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1))
225      END
226    ENDCASE
227    if keyword_set(nan) NE 0 then BEGIN
228      if nan NE 1 then BEGIN    ; If nan is not !values.f_nan
229; we put it at !values.f_nan
230        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $
231        ELSE notanumber = where(abs(res) GT abs(nan)/10.)
232        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan
233      ENDIF
234    ENDIF
235;---------------------------------------------------------------
236; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO AVERAGE = 1,
237; AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE reform(...,nx,ny,...) WHICH CAN
238; LOOK USELESS AT THE BEGINNING
239;---------------------------------------------------------------
240    if nx EQ 1 OR ny EQ 1 then BEGIN
241      res = reform(res, nx, ny, /over)
242      e1 =  reform(e1, nx, ny, /over)
243      e2 =  reform(e2, nx, ny, /over)
244    endif
245    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $
246      mask =  reform(mask, nx, ny, nz, /over)
247;---------------------------------------------------------------
248; II.3) Different types of average
249;---------------------------------------------------------------
250    mask = mask[*, *, 0]
251    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1
252    case 1 of
253      (dirx eq 1) and (diry eq 0) : begin
254        e = e1*mask
255        if keyword_set(integration) then divi = 1 $
256        else begin
257          divi = e
258          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
259          if ny EQ 1 then divi = reform(divi, nx, ny, /over)
260          divi = total(divi, 1)
261        endelse
262        res = res*e
263        if ny EQ 1 then res = reform(res, nx, ny, /over)
264        res = total(res, 1, nan = nan)/(divi > 1.)
265        if msknan[0] NE -1 then begin
266          testnan  = msknan*mask
267          if ny EQ 1 then testnan = reform(testnan, nx, ny, /over)
268          testnan = total(testnan, 1)+(total(mask, 1) EQ 0)
269        endif
270      end
271      (dirx eq 0) and (diry eq 1) : begin
272        e = e2*mask
273        if keyword_set(integration) then divi = 1 $
274        else begin
275          divi = e
276          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
277          if ny EQ 1 then divi = reform(divi, nx, ny, /over)
278          divi = total(divi, 2)
279        endelse
280        res = res*e
281        if ny EQ 1 then res = reform(res, nx, ny, /over)
282        res = total(res, 2, nan = nan)/(divi > 1.)
283        if msknan[0] NE -1 then begin
284          testnan  = msknan*mask
285          if ny EQ 1 then testnan  = reform(testnan, nx, ny, /over)
286          testnan = total(testnan, 2)+(total(mask, 2) EQ 0)
287        endif
288      end
289      (dirx eq 1) and (diry eq 1) : begin
290        if keyword_set(integration) then divi = 1 else BEGIN
291          IF msknan[0] NE -1 THEN divi = total(e1*e2*mask*msknan) $
292          ELSE divi = total(e1*e2*mask)
293        ENDELSE
294        res = total(res*e1*e2*mask, nan = nan)/(divi > 1.)
295        if msknan[0] NE -1 then begin
296          testnan  = msknan*mask
297          testnan = total(testnan)+(total(mask) EQ 0)
298        endif
299      end
300    endcase
301  endif
302;------------------------------------------------------------
303;------------------------------------------------------------
304; III) Case 3d arrays series (tab4d)
305;------------------------------------------------------------
306;------------------------------------------------------------
307  if (dim eq '3d') then begin
308;---------------------------------------------------------------
309; III.1) Verification of the coherence of the array to average size
310; Verification of the coherence between the array's size and the domain
311; defind by domdef
312; The input array must have either the total domain size (jpi,jpj,jpk)
313; or this one of the reduced domain (nx,ny,ny)
314;---------------------------------------------------------------
315    case 1 of
316      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk: $
317        res = tab[firstx:lastx, firsty:lasty, firstz:lastz]
318      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz: $
319        res = tab[firstx:lastx, firsty:lasty, *]
320      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz :res = tab
321      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk : $
322        res = tab[*, *, firstz:lastz]
323      else:BEGIN
324        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
325        return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1))
326      END
327    endcase
328    if keyword_set(nan) NE 0 then BEGIN
329      if nan NE 1 then BEGIN    ; if nan is not !values.f_nan
330; we put it at !values.f_nan
331        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $
332        ELSE notanumber = where(abs(res) GT abs(nan)/10.)
333        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan
334      ENDIF
335    ENDIF
336;---------------------------------------------------------------
337; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO AVERAGE = 1,
338; AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE reform(...,nx,ny,...) WHICH CAN
339; LOOK USELESS AT THE BEGINNING
340;---------------------------------------------------------------
341    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN
342      res = reform(res, nx, ny, nz, /over)
343      e1 =  reform(e1, nx, ny, /over)
344      e2 =  reform(e2, nx, ny, /over)
345    endif
346    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $
347      mask =  reform(mask, nx, ny, nz, /over)
348    IF keyword_set(key_partialstep) THEN BEGIN
349; the top of the ocean floor is
350      IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $
351      ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
352; we suppress columns with only ocean or land
353      good = where(bottom NE 0 AND bottom NE nz)
354; the bottom of the ocean in 3D index is:
355      bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny
356      IF good[0] NE -1 THEN bottom = bottom[good] $
357      ELSE bottom = -1
358    ENDIF ELSE bottom = -1
359;---------------------------------------------------------------
360; III.2) different average types
361;---------------------------------------------------------------
362    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1
363    case 1 of
364      (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin
365        e13 = e1[*]#replicate(1., nz)
366        e13 = reform(e13, nx, ny, nz, /over)
367        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
368          AND nx NE 1 THEN BEGIN
369          IF msknan[0] EQ -1 THEN BEGIN
370            msknan = replicate(1b, nx, ny, nz)
371            nan = 1
372          endif
373          msknan[bottom] = 0
374          res[bottom] = !values.f_nan
375        ENDIF
376        if keyword_set(integration) then divi = 1 else begin
377          divi = e13*mask
378          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
379          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
380          divi = total(divi, 1)
381        ENDELSE
382        res = res*e13*mask
383        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
384        res = total(res, 1, nan = nan)/(divi > 1.)
385        e13 = 1
386        if msknan[0] NE -1 then begin
387          testnan  = msknan*mask
388          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
389          testnan = total(testnan, 1)+(total(mask, 1) EQ 0)
390        endif
391      end
392      (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin
393        e23 = e2[*]#replicate(1., nz)
394        e23 = reform(e23, nx, ny, nz, /over)
395        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
396          AND ny NE 1 THEN BEGIN
397          IF msknan[0] EQ -1 THEN BEGIN
398            msknan = replicate(1b, nx, ny, nz)
399            nan = 1
400          endif
401          msknan[bottom] = 0
402          res[bottom] = !values.f_nan
403        ENDIF
404        if keyword_set(integration) then divi = 1 else begin
405          divi = e23*mask
406          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
407          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
408          divi = total(divi, 2)
409        ENDELSE
410        res = res*e23*mask
411        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
412        res = total(res, 2, nan = nan)/(divi > 1.)
413        e23 = 1
414        if msknan[0] NE -1 then begin
415          testnan  = msknan*mask
416          if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over)
417          testnan = total(testnan, 2)+(total(mask, 2) EQ 0)
418        endif
419      end
420      (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin
421        e33 = replicate(1, 1.*nx*ny)#e3
422        e33 = reform(e33, nx, ny, nz, /over)
423        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
424          IF keyword_set(wdepth) THEN $
425            e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
426          ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
427        ENDIF
428        if keyword_set(integration) then divi = 1 else begin
429          divi = e33*mask
430          if msknan[0] NE -1 then divi = temporary(divi)*msknan
431          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
432          divi = total(divi, 3)
433        ENDELSE
434        res = res*e33*mask
435        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
436        res = total(res, 3, nan = nan)/(divi > 1.)
437        e33 = 1
438        if msknan[0] NE -1 then begin
439          testnan  = msknan*mask
440          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
441          testnan = total(testnan, 3)+(total(mask, 3) EQ 0)
442        endif
443      end
444      (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin
445        e123 = (e1*e2)[*]#replicate(1., nz)
446        e123 = reform(e123, nx, ny, nz, /over)
447        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
448          AND nx*ny NE 1 THEN BEGIN
449          IF msknan[0] EQ -1 THEN BEGIN
450            msknan = replicate(1b, nx, ny, nz)
451            nan = 1
452          endif
453          msknan[bottom] = 0
454          res[bottom] = !values.f_nan
455        ENDIF
456        if keyword_set(integration) then divi = 1 else BEGIN
457          divi = e123*mask
458          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
459          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
460          divi = total(total(divi, 1), 1)
461        ENDELSE
462        res = res*e123*mask
463        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
464        res = total(total(res, 1, nan = nan), 1, nan = nan) / (divi > 1.)
465        e123 = 1
466        if msknan[0] NE -1 then begin
467          testnan  = msknan*mask
468          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
469          testnan = total(total(testnan, 1), 1)+(total(total(mask, 1), 1) EQ 0)
470        endif
471      end
472      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin
473        e133 = e1[*]#e3
474        e133 = reform(e133, nx, ny, nz, /over)
475        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
476          IF keyword_set(wdepth) THEN $
477            e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
478          ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
479        ENDIF
480        if keyword_set(integration) then divi = 1 else BEGIN
481          divi = e133*mask
482          if msknan[0] NE -1 then divi = temporary(divi)*msknan
483          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
484          divi = total(total(divi, 1), 2)
485        ENDELSE
486        res = res*e133*mask
487        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
488        res = total(total(res, 1, nan = nan), 2, nan = nan) / (divi > 1.)
489        e133 = 1
490        if msknan[0] NE -1 then begin
491          testnan  = msknan*mask
492          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
493          testnan = total(total(testnan, 1), 2)+(total(total(mask, 1), 2) EQ 0)
494        endif
495      end
496      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin
497        e233 = e2[*]#e3
498        e233 = reform(e233, nx, ny, nz, /over)
499        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
500          IF keyword_set(wdepth) THEN $
501            e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
502          ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
503        ENDIF
504        if keyword_set(integration) then divi = 1 else BEGIN
505          divi = e233*mask
506          if msknan[0] NE -1 then divi = temporary(divi)*msknan
507          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
508          divi = total(total(divi, 2), 2)
509        ENDELSE
510        res = res*e233*mask
511        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
512        res = total(total(res, 2, nan = nan), 2, nan = nan) / (divi > 1.)
513        e233 = 1
514        if msknan[0] NE -1 then begin
515          testnan  = msknan*mask
516          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
517          testnan = total(total(testnan, 2), 2)+(total(total(mask, 2), 2) EQ 0)
518        endif
519      end
520      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin
521        e1233 = (e1*e2)[*]#e3
522        e1233 = reform(e1233, nx, ny, nz, /over)
523        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
524          IF keyword_set(wdepth) THEN $
525            e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
526          ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
527        ENDIF
528        if keyword_set(integration) then divi = 1 else BEGIN
529          if msknan[0] NE -1 then divi = total(e1233*mask*msknan) $
530          ELSE divi = total(e1233*mask)
531        ENDELSE
532        res = total(res*e1233*mask, nan = nan) / (divi > 1.)
533        e1233 = 1
534        if msknan[0] NE -1 then begin
535          testnan  = msknan*mask
536          testnan = total(testnan)+(total(mask) EQ 0)
537        endif
538      end
539    endcase
540  endif
541;------------------------------------------------------------
542;------------------------------------------------------------
543;IV ) finishing
544;------------------------------------------------------------
545;------------------------------------------------------------
546; IV.1) We mask land by a value equal to 1.e+20
547;------------------------------------------------------------
548  valmask = 1e+20
549  terre = where(divi EQ 0)
550  IF terre[0] NE -1 THEN BEGIN
551    res[terre] = 1e+20
552  ENDIF 
553;------------------------------------------------------------
554; IV.2) We replace, when nan equal 1, !values.f_nan by nan
555;------------------------------------------------------------
556  if keyword_set(nan) NE 0 then BEGIN
557    puttonan = where(testnan EQ 0)
558    if puttonan[0] NE -1 then res[puttonan] = !values.f_nan
559    if nan NE 1 then BEGIN
560      notanumber = where(finite(res) eq 0)
561      if notanumber[0] NE -1 then res[notanumber] = nan
562    ENDIF
563  ENDIF
564;------------------------------------------------------------
565; IV.3) We replace in the domain whch was defined at the entry of average
566;------------------------------------------------------------
567  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
568;------------------------------------------------------------
569  if keyword_set(key_performance) THEN print, 'temps moyenne', systime(1)-tempsun
570  return, res
571;------------------------------------------------------------
572;------------------------------------------------------------
573end
Note: See TracBrowser for help on using the repository browser.