source: trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro

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

fix thanks to coding rules; typo; dupe empty lines; trailing blanks

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