Context Navigation

grossemoyenne.pro @ 370

Last change on this file since 370 was 370, checked in by pinsard, 16 years ago
improvemnts of headers (typo, links)
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 27.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format