;+
;
; @file_comments
; averages a 3- or 4-d time series field over a selected
; geographical area or along the time axis. For one or more
; selected axes (x, y, z, t)
;
; @categories
;
; @param TAB {in}{required}
; 3 or 4d field
;
; @param DIREC {in}{required}
; 'x' 'y' 'z' 't' 'xy' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt'
; 'xyt' 'xzt' 'yzt' or 'xyzt'
;
; @keyword BOXZOOM
; [xmin,xmax,ymin,ymax (,(zmin,)zmax)] to more details, see domdef
; boxzoom can have 5 forms :
; [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2],
; [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1, lat2, vert1,vert2]
;
; @keyword NAN
; not a number, we activate it if we want to average without considerate some
; masked values of TAB.
; If masked values of TAB are values consecrated by IDL(!values.f_nan), we
; just have to put NAN.
; If masked values of TAB are valued a (a must be different of 1,
; corresponding to nan = !values.f_nan and of 0, which desactivate nan).
; We have to put NAN=a.
; Comment: In output, result points which are NAN will be valued a or
; !values.f_nan.
;
; @keyword NODOMDEF
; We activate it if we do not want to pass in domdef even if the
; keyword boxzoom is present (like when grossemoyenne is called via
; checkfield)
;
; @keyword INTEGRATION
; To make an integral rather than an average
;
; @keyword SPATIALFIRST
; when performing at the same time
; spatial and temporal mean, grossemoyenne is assuming
; that the mask is not changing with the time. In
; consequence, grossemoyenne performs temporal mean
; first and then call moyenne. Activate /SPATIALFIRST if
; you want to perform the spatial mean before the
; temporal mean. Note that if NAN is activated, then
; SPATIALFIRST is activated automatically.
;
; @keyword TEMPORALFIRST
; to force to perform first temporal
; mean even if nan is activated (see SPATIALFIRST explanations...)
;
; @keyword WDEPTH
; to specify that the field is at W depth instead of T
; depth (automatically activated if vargrid eq 'W')
;
; @uses
; result:un tableau
; common
; domdef
;
; @restrictions
; Put values corresponding to land at 1.e20
;
; @history
; Jerome Vialard (jv\@lodyc.jussieu.fr)
; 2/7/98
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; adaptation array containing a temporal dimension
; 14/8/98
; 15/1/98
; 12/3/99 adaptation for NAN and utilization of TEMPORARY
;
; @version
; $Id$
;-
;
FUNCTION grossemoyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $
, NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $
, SPATIALFIRST = spatialfirst, TEMPORALFIRST = temporalfirst $
, _EXTRA = ex
;
compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
ENDIF
;---------------------------------------------------------
tempsun = systime(1) ; For key_performance
;------------------------------------------------------------
;I) preliminaries
;------------------------------------------------------------
dirt = 0
dirx = 0
diry = 0
dirz = 0
dim = 'aa'
;------------------------------------------------------------
; I.1) Directions(s) we follow to integrate
;------------------------------------------------------------
if ( strpos(direc, 't') ge 0 ) then dirt = 1
if ( strpos(direc, 'x') ge 0 ) then dirx = 1
if ( strpos(direc, 'y') ge 0 ) then diry = 1
if ( strpos(direc, 'z') ge 0 ) then dirz = 1
IF keyword_set(NAN) AND (dirx EQ 1 OR diry EQ 1 OR dirz EQ 1) $
THEN spatialfirst = 1
IF keyword_set(temporalfirst) THEN spatialfirst = 0
;------------------------------------------------------------
; I.2) verification of the input array's size
;------------------------------------------------------------
taille = size(tab)
case 1 of
taille[0] eq 1 :return, report('array has only one dimension, not implemented!')
taille[0] eq 2 :return, report('array has only two dimensions, not implemented!')
taille[0] eq 3 :BEGIN
dim = '3d'
if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab
END
taille[0] eq 4 :BEGIN
dim = '4d'
if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab
END
else : return, report('array must have 3 or 4 dimensions if there is not time dimension use moyenne')
endcase
;------------------------------------------------------------
; I.3) Obtainment of scale's factors and of the mask on the subdomain concernedby the average.
; Redefinition of the domain ajusted at boxzoom (at 6 elements)
; This will allowed us to calculate only in the domain concerned by the average.
; Domdef, followed by grid give us all arrays of the grid on the subdomain
;------------------------------------------------------------
if keyword_set(boxzoom) then BEGIN
Case 1 Of
N_Elements(Boxzoom) Eq 1: bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
N_Elements(Boxzoom) Eq 2: bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
N_Elements(Boxzoom) Eq 4: bte = [Boxzoom, vert1, vert2]
N_Elements(Boxzoom) Eq 5: bte = [Boxzoom[0:3], 0, Boxzoom[4]]
N_Elements(Boxzoom) Eq 6: bte = Boxzoom
Else: return, report('Wrong Definition of Boxzoom')
endcase
if NOT keyword_set(nodomdef) then BEGIN
savedbox = 1b
saveboxparam, 'boxparam4grmoyenne.dat'
domdef, bte, GRIDTYPE = vargrid, _extra = ex
ENDIF
ENDIF
;---------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays...
;---------------------------------------------------------------
grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth
;------------------------------------------------------------
; I.4) if dirt equal 1, we make the temporal average and we send it in moyenne
;------------------------------------------------------------
if dirt EQ 1 AND NOT keyword_set(spatialfirst) then begin
if dim EQ 3d then BEGIN
case 1 of
taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $
res = tab[firstx:firstx+nx-1 $
, firsty:firsty+ny-1, *]
taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res = tab
else:BEGIN
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
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))
END
ENDCASE
if keyword_set(integration) then begin
res = total(res, 3, nan = nan)
ENDIF ELSE BEGIN
if keyword_set(nan) then BEGIN
divi = finite(res)
divi = total(temporary(divi), 3)
notanum = where(divi EQ 0)
res = total(res, 3, nan = keyword_set(nan))/ (1 > divi)
if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan
ENDIF ELSE res = total(res, 3)/(1.*taille[3])
ENDELSE
ENDIF ELSE BEGIN
case 1 of
taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $
res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *]
taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
res = tab[firstx:lastx, firsty:lasty, *, *]
taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz and taille[4] eq jpt:res = tab
taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $
res = tab[*, *, firstz:lastz, *]
else:BEGIN
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $
+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $
+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $
+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $
+strtrim(taille[4], 1))
END
endcase
if keyword_set(integration) then begin
res = total(res, 4, nan = nan)
ENDIF ELSE BEGIN
if keyword_set(nan) then begin
divi = finite(res)
divi = total(temporary(divi), 4)
notanum = where(divi EQ 0)
res = total(res, 4, /nan)/(1 > divi)
if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan
ENDIF ELSE res = total(res, 4)/(1.*taille[4])
ENDELSE
ENDELSE
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
return, moyenne(temporary(res), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex)
ENDIF ELSE res = tab
IF jpt EQ 1 THEN BEGIN
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
return, moyenne(reform(res, /over), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex)
END
;------------------------------------------------------------
;------------------------------------------------------------
; II) Case 2d arrays series (tab3d)
;------------------------------------------------------------
;------------------------------------------------------------
if (dim eq '3d') then begin
;---------------------------------------------------------------
; II.1) verification of the coherence of the array's size to average
; Verification of the coherence between the array's size and the domain defined by domdef
; The input array must have either the total domain's size (jpi,jpj,jpt) or
; this one of the reduced domain (nx,ny,jpt)
;---------------------------------------------------------------
case 1 of
taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $
res = tab[firstx:firstx+nx-1 $
, firsty:firsty+ny-1, *]
taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res = tab
else:BEGIN
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
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))
enD
endcase
if keyword_set(nan) NE 0 then BEGIN
if nan NE 1 then BEGIN ; If nan is not !values.f_nan
; we put it at !values.f_nan
if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $
ELSE notanumber = where(abs(res) GT abs(nan)/10.)
if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan
ENDIF
ENDIF
;---------------------------------------------------------------
; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO
; AVERAGE = 1, AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE
; reform(...,nx,ny,...) WHICH CAN LOOK USELESS AT THE BEGINNING
;---------------------------------------------------------------
if nx EQ 1 OR ny EQ 1 then BEGIN
res = reform(res, nx, ny, jpt, /over)
e1 = reform(e1, nx, ny, /over)
e2 = reform(e2, nx, ny, /over)
endif
if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $
mask = reform(mask, nx, ny, nz, /over)
;---------------------------------------------------------------
; II.3) Different types of average
;---------------------------------------------------------------
if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1
mask = mask[*, *, 0]
case 1 of
(dirx eq 1) and (diry eq 0) : begin
e = temporary(e1)*temporary(mask)
echelle = (temporary(e))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE BEGIN
IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $
ELSE divi = total(echelle, 1)
ENDELSE
res = total(temporary(res)*echelle, 1, nan = nan)/(divi > 1.)
if msknan[0] NE -1 then BEGIN
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(temporary(testnan), 1) $
+(total(temporary(echelle), 1) EQ 0)
endif
end
(dirx eq 0) and (diry eq 1) : begin
e = temporary(e2)*temporary(mask)
if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over)
echelle = (temporary(e))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE BEGIN
IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $
ELSE divi = total(echelle, 2)
ENDELSE
res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1.)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(temporary(testnan), 2) $
+(total(temporary(echelle), 2) EQ 0)
endif
end
(dirx eq 1) and (diry eq 1) : begin
echelle = (temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE BEGIN
IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $
ELSE divi = total(total(echelle, 1), 1)
ENDELSE
res = total(temporary(total(temporary(res)*echelle, 1, nan = nan)), 1, nan = nan)/(divi > 1.)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(total(temporary(testnan), 1), 1) $
+(total(total(temporary(echelle), 1), 1) EQ 0)
endif
end
endcase
endif
;------------------------------------------------------------
;------------------------------------------------------------
; III) Case 3d arrays series (tab4d)
;------------------------------------------------------------
if (dim eq '4d') then begin
;---------------------------------------------------------------
; III.1) Verification of the coherence of the array to average size
; Verification of the coherence between the array's size and the domain
; defind by domdef
; The input array must have either the total domain size (jpi,jpj,jpk,jpt)
; or this one of the reduced domain (nx,ny,ny,jpt)
;---------------------------------------------------------------
case 1 of
taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $
res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *]
taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
res = tab[firstx:lastx, firsty:lasty, *, *]
taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz and taille[4] eq jpt:res = tab
taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $
res = tab[*, *, firstz:lastz, *]
else:BEGIN
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $
+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $
+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $
+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $
+strtrim(taille[4], 1))
END
endcase
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)
if keyword_set(nan) NE 0 then BEGIN
if nan NE 1 then BEGIN ; if nan is not !values.f_nan
; we put it at !values.f_nan
if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $
ELSE notanumber = where(abs(res) GT abs(nan)/10.)
if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan
ENDIF
ENDIF
;---------------------------------------------------------------
; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO
; AVERAGE = 1, AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE
; reform(...,nx,ny,...) WHICH CAN LOOK USELESS AT THE BEGINNING
;---------------------------------------------------------------
if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN
res = reform(res, nx, ny, nz, jpt, /over)
mask = reform(mask, nx, ny, nz, /over)
ENDIF
IF keyword_set(key_partialstep) THEN BEGIN
; the top of the ocean floor is
IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $
ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
; we suppress columns with only ocean or land
good = where(bottom NE 0 AND bottom NE nz)
; the bottom of the ocean in 3D index is:
bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny
IF good[0] NE -1 THEN bottom = bottom[good] $
ELSE bottom = -1
ENDIF ELSE bottom = -1
;---------------------------------------------------------------
; III.2) different average types
;---------------------------------------------------------------
IF keyword_set(nan) NE 0 THEN msknan = finite(res) ELSE msknan = -1
case 1 of
(dirx eq 1) and (diry eq 0) and (dirz eq 0) : BEGIN
e13 = (temporary(e1))[*]#replicate(1., nz)
e13 = reform(e13, nx, ny, nz, /over)
echelle = (temporary(e13)*temporary(mask))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
AND nx NE 1 THEN BEGIN
IF msknan[0] EQ -1 THEN BEGIN
msknan = replicate(1b, nx, ny, nz, jpt)
nan = 1
ENDIF
bottom = bottom#replicate(1, jpt) $ ; 4D bottom!
+ replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt))
msknan[bottom] = 0
res[temporary(bottom)] = !values.f_nan
ENDIF
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $
ELSE divi = total(echelle, 1)
endelse
res = temporary(res)*echelle
res = total(temporary(res), 1, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(temporary(testnan), 1) $
+(total(temporary(echelle), 1) EQ 0)
endif
end
(dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin
e23 = temporary(e2[*])#replicate(1., nz)
e23 = reform(e23, nx, ny, nz, /over)
echelle = (temporary(e23)*temporary(mask))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
AND ny NE 1 THEN BEGIN
IF msknan[0] EQ -1 THEN BEGIN
msknan = replicate(1b, nx, ny, nz, jpt)
nan = 1
endif
bottom = bottom#replicate(1, jpt) $ ; 4D bottom!
+ replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt))
msknan[bottom] = 0
res[temporary(bottom)] = !values.f_nan
ENDIF
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $
ELSE divi = total(echelle, 2)
endelse
res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(temporary(testnan), 2) $
+(total(temporary(echelle), 2) EQ 0)
endif
end
(dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin
e33 = replicate(1, 1.*nx*ny)#e3
e33 = reform(e33, nx, ny, nz, /over)
IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
IF keyword_set(wdepth) THEN $
e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
ENDIF
echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $
ELSE divi = total(echelle, 3)
endelse
res = total(temporary(res)*echelle, 3, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(temporary(testnan), 3) $
+(total(temporary(echelle), 3) EQ 0)
endif
end
(dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin
e13 = e1[*]#replicate(1., nz)
e13 = reform(e13, nx, ny, nz, /over)
e23 = e2[*]#replicate(1., nz)
e23 = reform(e23, nx, ny, nz, /over)
echelle = (temporary(e13)*temporary(e23)*temporary(mask))[*]#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
AND nx*ny NE 1 THEN BEGIN
IF msknan[0] EQ -1 THEN BEGIN
msknan = replicate(1b, nx, ny, nz, jpt)
nan = 1
endif
bottom = bottom#replicate(1, jpt) $ ; 4D bottom!
+ replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt))
msknan[bottom] = 0
res[temporary(bottom)] = !values.f_nan
ENDIF
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $
ELSE divi = total(total(echelle, 1), 1)
endelse
res = total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(total(temporary(testnan), 1), 1) $
+(total(total(temporary(echelle), 1), 1) EQ 0)
endif
end
(dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin
e133 = e1[*]#e3
IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
IF keyword_set(wdepth) THEN $
e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
ENDIF
echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $
ELSE divi = total(total(echelle, 1), 2)
endelse
res = total(total(temporary(res)*echelle, 1, nan = nan), 2, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(total(temporary(testnan), 1), 2) $
+(total(total(temporary(echelle), 1), 2) EQ 0)
endif
end
(dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin
e233 = e2[*]#e3
IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
IF keyword_set(wdepth) THEN $
e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
ENDIF
echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $
ELSE divi = total(total(echelle, 2), 2)
endelse
res = total(total(temporary(res)*echelle, 2, nan = nan), 2, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(total(temporary(testnan), 2), 2) $
+(total(total(temporary(echelle), 2), 2) EQ 0)
endif
end
(dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin
e1233 = (e1[*]*e2[*])#e3
IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
IF keyword_set(wdepth) THEN $
e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
ENDIF
echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt)
echelle = reform(echelle, nx, ny, nz, jpt, /over)
if keyword_set(integration) then divi = 1 ELSE begin
IF msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $
ELSE divi = total(total(total(echelle, 1), 1), 1)
endelse
res = total(total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan), 1, nan = nan)/(divi > 1)
if msknan[0] NE -1 then begin
echelle = temporary(echelle) NE 0
testnan = temporary(msknan)*echelle
testnan = total(total(total(temporary(testnan), 1), 1), 1) $
+(total(total(total(temporary(echelle), 1), 1), 1) EQ 0)
endif
end
endcase
endif
;------------------------------------------------------------
if dirt EQ 1 AND keyword_set(spatialfirst) then BEGIN
IF (reverse(size(res, /dimension)))[0] NE jpt THEN BEGIN
ras = report('the last dimension of res is not equal to jpt: '+strtrim(jpt, 2))
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
return, -1
ENDIF
tdim = size(res, /n_dimensions)
if keyword_set(integration) then res = total(res, tdim, nan = nan) $
ELSE BEGIN
if keyword_set(nan) then BEGIN
testnan = testnan < 1
testnan = total(temporary(testnan), tdim)
divi = testnan
ENDIF ELSE divi = jpt
res = total(res, tdim, nan = nan)/(1 > divi)
ENDELSE
ENDIF
;------------------------------------------------------------
;------------------------------------------------------------
;IV ) finishing
;------------------------------------------------------------
;------------------------------------------------------------
; IV.1) We mask land by a value equal to 1.e+20
;------------------------------------------------------------
valmask = 1e+20
terre = where(divi EQ 0)
IF terre[0] NE -1 THEN BEGIN
res[temporary(terre)] = 1e+20
ENDIF
;------------------------------------------------------------
; IV.2) We replace, when nan equal 1, !values.f_nan by nan
;------------------------------------------------------------
if keyword_set(nan) NE 0 then BEGIN
puttonan = where(temporary(testnan) EQ 0)
if puttonan[0] NE -1 then res[temporary(puttonan)] = !values.f_nan
if nan NE 1 then BEGIN
notanumber = where(finite(res) eq 0)
if notanumber[0] NE -1 then res[temporary(notanumber)] = nan
ENDIF
ENDIF
;------------------------------------------------------------
; IV.3) We replace in the domain whch was defined at the entry of average
;------------------------------------------------------------
if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
;------------------------------------------------------------
if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun
return, res
;------------------------------------------------------------
;------------------------------------------------------------
end