Context Navigation

norme.pro @ 157

Last change on this file since 157 was 157, checked in by navarro, 18 years ago
header improvements + xxx doc
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 22.1 KB

Rev	Line
[2]	1	;------------------------------------------------------------
	2	;------------------------------------------------------------
	3	;------------------------------------------------------------
	4	;+
	5	;
[142]	6	; @file_comments
	7	; calculate the norme of a field of vectors, then make a possible average.
	8	; Comment 1: The field of vector can be, 2d:xy, 3d: xyz or xyt,
	9	; 4d: xyzt
	10	; Comment 2:
	11	; The calculation of the norme is made before the possible spatial or
	12	; temporal average because the average of the norme is not equal to the
	13	; norme of averages
	14
[2]	15	;
[142]	16	; @categories
[157]	17	; Calculation
[2]	18	;
[142]	19	; @param COMPOSANTEU {in}{required}
	20	; an 2d, 3d or 4d array
[2]	21	;
[142]	22	; @param COMPOSANTEV {in}{required}
	23	; an 2d, 3d or 4d array
[2]	24	;
[142]	25	; @keyword BOXZOOM
	26	; boxzoom on which do the average (by default the domain selected
	27	; by the last domdef done)
[2]	28	;
[142]	29	; @keyword DIREC
	30	; 't' 'x' 'y' 'z' 'xys' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
	31	; 'xzt' 'yzt' 'xyzt' Direction on which do averages
[2]	32	;
[142]	33	; @returns
	34	; Array to trace with plt, pltz or pltt.
[2]	35	;
[142]	36	; @uses
	37	; common.pro
[2]	38	;
[142]	39	; @restrictions
	40	; The norme is calculated on points TTo do this calculation, we average
	41	; field U and Von points T before calculate the norme. At the edge of
	42	; coast and of domain, we can not calculate fields U and V at points T,
	43	; that is why these points are at value !values.f_nan.
[2]	44	;
[142]	45	; When we calculate on a reduce geographic domain, field U and V have not
	46	; necessarily the same number of point. In this case, we recut U and V to
	47	; keep only common points. We profit of this to redo a domdef which redefine
	48	; a geographic domain on which fields U and V are extracted on same points
[2]	49	;
[142]	50	; @restrictions
	51	; To know what type of array we work with, we test its size and dates
	52	; gave by time[0] and time[jpt-1] to know if thee is a temporal dimension.
	53	; Before to start norme, make sure that time and jpt are defined how
	54	; they have to!
[2]	55	;
[142]	56	; @examples
	57	; To calculate the average of the norme of streams on all the domain
	58	; between 0 et 50:
	59	; IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz')
[2]	60	;
[142]	61	; @history
[157]	62	; Sebastien Masson (smasson\@lodyc.jussieu.fr)
[142]	63	; 9/6/1999
[2]	64	;
[142]	65	; @version
	66	; $Id$
[2]	67	;
	68	;-
	69	;------------------------------------------------------------
	70	;------------------------------------------------------------
	71	;------------------------------------------------------------
[25]	72	FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _extra = ex
	73	;---------------------------------------------------------
[114]	74	;
	75	compile_opt idl2, strictarrsubs
	76	;
[25]	77	@cm_4mesh
	78	@cm_4data
	79	@cm_4cal
	80	IF NOT keyword_set(key_forgetold) THEN BEGIN
	81	@updatenew
	82	@updatekwd
	83	ENDIF
	84	;---------------------------------------------------------
[142]	85	tempsun = systime(1) ; To key_performance
[25]	86	;
	87	IF finite(glamu[0])finite(gphiu[0])finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
	88	return, report(['This version of norme is based on Arakawa C-grid.' $
	89	, 'U and V grids must therefore be defined'])
	90	;
[2]	91	;------------------------------------------------------------
[25]	92	if keyword_set(boxzoom) then BEGIN
	93	Case 1 Of
	94	N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
	95	N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
	96	N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
	97	N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
	98	N_Elements(Boxzoom) Eq 6:bte = Boxzoom
	99	Else: return, report('Mauvaise Definition de Boxzoom')
	100	ENDCASE
	101	domdef, boxzoom
	102	ENDIF
	103	;------------------------------------------------------------
[2]	104	if NOT keyword_set(direc) then direc = 0
[142]	105	; construction of u and v at points T
[2]	106	u = litchamp(composanteu)
	107	v = litchamp(composantev)
	108	date1 = time[0]
	109	if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
	110
	111	if (size(u))[0] NE (size(v))[0] then return, -1
	112
	113	vargrid='T'
	114	varname = 'norme '
	115	valmask = 1e20
	116	;
	117	grilleu = litchamp(composanteu, /grid)
	118	if grilleu EQ '' then grilleu = 'U'
	119	grillev = litchamp(composantev, /grid)
	120	if grillev EQ '' then grillev = 'V'
	121	IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
	122	IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN
	123	interpolle = 0
	124	return, report('cas non code mais facile a faire!')
	125	ENDIF ELSE interpolle = 1
	126	if keyword_set(inverse) then begin
	127	rien = u
	128	u = v
	129	v = rien
	130	endif
	131
	132
	133	;------------------------------------------------------------
[142]	134	; We find common points between u and v
[2]	135	;------------------------------------------------------------
[25]	136	indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
	137	indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
[2]	138	indicex = inter(indicexu, indicexv)
[25]	139	indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
	140	indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
[2]	141	indicey = inter(indiceyu, indiceyv)
	142	nx = n_elements(indicex)
	143	ny = n_elements(indicey)
	144	;----------------------------------------------------------------------------
	145	case 1 of
	146	;----------------------------------------------------------------------------
	147	;----------------------------------------------------------------------------
	148	;xyz
	149	;----------------------------------------------------------------------------
	150	;----------------------------------------------------------------------------
	151	(size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
	152	;----------------------------------------------------------------------------
	153	indice2d = lindgen(jpi, jpj)
	154	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
	155	indice3d = lindgen(jpi, jpj, jpk)
[25]	156	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt]
[2]	157	;------------------------------------------------------------
[142]	158	; extraction of u and v on the appropriated domain
[2]	159	;------------------------------------------------------------
	160	case 1 of
	161	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
	162	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
	163	case (size(u))[3] OF
	164	nzt:BEGIN
	165	if nxu NE nx then $
[25]	166	if indicex[0] EQ firstxu then u = u[0:nx-1,,] ELSE u = u[1: nx, ,]
[2]	167	IF nxv NE nx THEN $
[25]	168	if indicex[0] EQ firstxv then v = v[0:nx-1,,] ELSE v = v[1: nx, ,]
[2]	169	IF nyu NE ny THEN $
[25]	170	if indicey[0] EQ firstyu then u = u[,0:ny-1,] ELSE u = u[, 1: ny,]
[2]	171	IF nyv NE ny THEN $
[25]	172	if indicey[0] EQ firstyv then v = v[,0:ny-1,] ELSE v = v[, 1: ny,]
[2]	173	end
	174	jpk:BEGIN
	175	if nxu NE nx then $
[25]	176	if indicex[0] EQ firstxu then u = u[0:nx-1, ,firstzt:lastzt] ELSE u = u[1: nx, ,firstzt:lastzt]
[2]	177	IF nxv NE nx THEN $
[25]	178	if indicex[0] EQ firstxv then v = v[0:nx-1, ,firstzt:lastzt] ELSE v = v[1: nx, ,firstzt:lastzt]
[2]	179	IF nyu NE ny THEN $
[25]	180	if indicey[0] EQ firstyu then u = u[, 0:ny-1,firstzt:lastzt] ELSE u = u[, 1: ny,firstzt:lastzt]
[2]	181	IF nyv NE ny THEN $
[25]	182	if indicey[0] EQ firstyv then v = v[, 0:ny-1,firstzt:lastzt] ELSE v = v[, 1: ny,firstzt:lastzt]
[2]	183	end
	184	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
	185	endcase
	186	END
	187	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
	188	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
	189	u = u[indice3d]
	190	v = v[indice3d]
	191	END
	192	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
	193	endcase
	194	;------------------------------------------------------------------
[142]	195	; We reshape u and v to make sure that no dimension has been erased
[2]	196	;------------------------------------------------------------------
	197	if nzt EQ 1 then begin
	198	u = reform(u, nx, ny, nzt, /over)
	199	v = reform(v, nx, ny, nzt, /over)
	200	endif
	201	;------------------------------------------------------------------
[142]	202	; construction of u and v at points T
[2]	203	;-----------------------------------------------------------
[114]	204	a=u[0,,]
[2]	205	u=(u+shift(u,1,0,0))/2.
[114]	206	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,,]=a
	207	a=v[,0,]
[2]	208	v=(v+shift(v,0,1,0))/2.
[114]	209	if NOT keyword_set(key_periodic) OR nx NE jpi then v[,0,]=a
[2]	210	;----------------------------------------------------------------------------
[142]	211	; attribution of the mask and of logitude and latitude arrays
[2]	212	;----------------------------------------------------------------------------
	213	mask = tmask[indice3d]
	214	if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
	215	;-----------------------------------------------------------
	216	if n_elements(valmask) EQ 0 THEN valmask = 1e20
	217	landu = where(u GE valmask/10)
	218	if landu[0] NE -1 then u[landu] = 0
	219	landv = where(v GE valmask/10)
	220	if landv[0] NE -1 then v[landv] = 0
	221	res=sqrt(u^2+v^2)
[114]	222	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,, ]=!values.f_nan
	223	res[,0, ]=!values.f_nan
[2]	224	mask = where(mask eq 0)
[114]	225	IF mask[0] NE -1 THEN res[mask] = valmask
[142]	226	; All kind of average
[25]	227	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
	228	if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
[2]	229	;-----------------------------------------------------------
	230	END
	231	;----------------------------------------------------------------------------
	232	;----------------------------------------------------------------------------
	233	;xyt
	234	;----------------------------------------------------------------------------
	235	;----------------------------------------------------------------------------
	236	date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
	237	indice2d = lindgen(jpi, jpj)
	238	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
	239	;------------------------------------------------------------
[142]	240	; extraction of u and v on the appropriated domain
[2]	241	;------------------------------------------------------------
	242	case 1 of
	243	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
	244	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
	245	if nxu NE nx then $
[25]	246	if indicex[0] EQ firstxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
[2]	247	IF nxv NE nx THEN $
[25]	248	if indicex[0] EQ firstxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
[2]	249	IF nyu NE ny THEN $
[25]	250	if indicey[0] EQ firstyu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
[2]	251	IF nyv NE ny THEN $
[25]	252	if indicey[0] EQ firstyv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
[2]	253	END
	254	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
	255	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
	256	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
	257	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
	258	END
	259	ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
	260	endcase
	261	;------------------------------------------------------------------
[142]	262	; construction of u and v at points T
[2]	263	;-----------------------------------------------------------
[114]	264	a=u[0,,]
[2]	265	u=(u+shift(u,1,0,0))/2.
[114]	266	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,,]=a
	267	a=v[,0,]
[2]	268	v=(v+shift(v,0,1,0))/2.
[114]	269	if NOT keyword_set(key_periodic) OR nx NE jpi then v[,0,]=a
[2]	270	;----------------------------------------------------------------------------
[142]	271	; attribution of the mask and of longitude and latitude arrays.
	272	; We recover the complete grid to establish a big mask extent in the four
	273	; direction to cover pointsfor which a land point has been
	274	; considerated (make a small drawing)
[2]	275	;----------------------------------------------------------------------------
[25]	276	mask = tmask[indice2d+jpijpjfirstzt]
[2]	277	if ny EQ 1 then mask = reform(mask, nx, ny, /over)
	278	;-----------------------------------------------------------
[142]	279	; construction of land containing all points to mask
[2]	280	;-----------------------------------------------------------
	281	if n_elements(valmask) EQ 0 THEN valmask = 1e20
	282	landu = where(u GE valmask/10)
	283	if landu[0] NE -1 then u[landu] = 0
	284	landv = where(v GE valmask/10)
	285	if landv[0] NE -1 then v[landv] = 0
	286	res=sqrt(u^2+v^2)
[114]	287	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,, ]=!values.f_nan
	288	res[,0, ]=!values.f_nan
[2]	289	mask = where(mask eq 0)
	290	IF mask[0] NE -1 THEN BEGIN
	291	coeftps = lindgen(jpt)nxny
	292	coeftps = replicate(1, n_elements(mask))#coeftps
	293	mask = (temporary(mask))[*]#replicate(1, jpt)
	294	mask =temporary(mask[]) + temporary(coeftps[])
[114]	295	res[temporary(mask)] = valmask
[2]	296	ENDIF
	297	; moyennes en tous genres
[25]	298	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
	299	if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
[2]	300	END
	301	;----------------------------------------------------------------------------
	302	;----------------------------------------------------------------------------
	303	;xyzt
	304	;----------------------------------------------------------------------------
	305	;----------------------------------------------------------------------------
	306	date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
	307	indice2d = lindgen(jpi, jpj)
	308	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
	309	indice3d = lindgen(jpi, jpj, jpk)
[25]	310	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt]
[2]	311	;------------------------------------------------------------
[142]	312	; extraction of u and v on the appropriated domain
[2]	313	;------------------------------------------------------------
	314	case 1 of
	315	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
	316	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
	317	case (size(u))[3] OF
	318	nzt:BEGIN
	319	if nxu NE nx then $
[25]	320	if indicex[0] EQ firstxu then u = u[0:nx-1,,,] ELSE u = u[1: nx, ,,]
[2]	321	IF nxv NE nx THEN $
[25]	322	if indicex[0] EQ firstxv then v = v[0:nx-1,,,] ELSE v = v[1: nx, ,,]
[2]	323	IF nyu NE ny THEN $
[25]	324	if indicey[0] EQ firstyu then u = u[,0:ny-1,,] ELSE u = u[, 1: ny,,]
[2]	325	IF nyv NE ny THEN $
[25]	326	if indicey[0] EQ firstyv then v = v[,0:ny-1,,] ELSE v = v[, 1: ny,,]
[2]	327	end
	328	jpk:BEGIN
	329	if nxu NE nx then $
[25]	330	if indicex[0] EQ firstxu then u = u[0:nx-1, ,firstzt:lastzt,] ELSE u = u[1: nx, ,firstzt:lastzt,]
[2]	331	IF nxv NE nx THEN $
[25]	332	if indicex[0] EQ firstxv then v = v[0:nx-1, ,firstzt:lastzt,] ELSE v = v[1: nx, ,firstzt:lastzt,]
[2]	333	IF nyu NE ny THEN $
[25]	334	if indicey[0] EQ firstyu then u = u[, 0:ny-1,firstzt:lastzt,] ELSE u = u[, 1: ny,firstzt:lastzt,]
[2]	335	IF nyv NE ny THEN $
[25]	336	if indicey[0] EQ firstyv then v = v[, 0:ny-1,firstzt:lastzt,] ELSE v = v[, 1: ny,firstzt:lastzt,]
[2]	337	end
	338	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
	339	endcase
	340	END
	341	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
	342	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
[25]	343	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *]
	344	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *]
[2]	345	END
	346	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
	347	endcase
	348	;------------------------------------------------------------------
[142]	349	; construction of u and v at points T
[2]	350	;-----------------------------------------------------------
[114]	351	a=u[0,,,*]
[2]	352	u=(u+shift(u,1,0,0,0))/2.
[114]	353	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,,,*]=a
	354	a=v[,0,,*]
[2]	355	v=(v+shift(v,0,1,0,0))/2.
[114]	356	if NOT keyword_set(key_periodic) OR nx NE jpi then v[,0,,*]=a
[2]	357	;----------------------------------------------------------------------------
[142]	358	; attribution of the mask and of logitude and latitude arrays
[2]	359	;----------------------------------------------------------------------------
	360	mask = tmask[indice3d]
	361	if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
	362	;-----------------------------------------------------------
	363	if n_elements(valmask) EQ 0 THEN valmask = 1e20
	364	landu = where(u GE valmask/10)
	365	if landu[0] NE -1 then u[landu] = 0
	366	landv = where(v GE valmask/10)
	367	if landv[0] NE -1 then v[landv] = 0
	368	res=sqrt(u^2+v^2)
[114]	369	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,, , *]=!values.f_nan
	370	res[,0, , *]=!values.f_nan
[2]	371	mask = where(mask eq 0)
	372	IF mask[0] NE -1 THEN BEGIN
	373	coeftps = lindgen(jpt)nxny*nzt
	374	coeftps = replicate(1, n_elements(mask))#coeftps
	375	mask = (temporary(mask))[*]#replicate(1, jpt)
	376	mask =temporary(mask[]) + temporary(coeftps[])
[114]	377	res[temporary(mask)] = valmask
[2]	378	ENDIF
[142]	379	; All kind of average
[25]	380	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
	381	if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
[2]	382	END
	383	;----------------------------------------------------------------------------
	384	;----------------------------------------------------------------------------
	385	;xy
	386	;----------------------------------------------------------------------------
	387	;----------------------------------------------------------------------------
	388	ELSE:BEGIN ;xy
	389	indice2d = lindgen(jpi, jpj)
	390	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
	391	;------------------------------------------------------------
[142]	392	; extraction of u and v on the appropriated domain
[2]	393	;------------------------------------------------------------
	394	case 1 of
	395	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
	396	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
	397	if nxu NE nx then $
[25]	398	if indicex[0] EQ firstxu then u = u[0:nx-1, ] ELSE u = u[1: nx, ]
[2]	399	IF nxv NE nx THEN $
[25]	400	if indicex[0] EQ firstxv then v = v[0:nx-1, ] ELSE v = v[1: nx, ]
[2]	401	IF nyu NE ny THEN $
[25]	402	if indicey[0] EQ firstyu then u = u[, 0:ny-1] ELSE u = u[, 1: ny]
[2]	403	IF nyv NE ny THEN $
[25]	404	if indicey[0] EQ firstyv then v = v[, 0:ny-1] ELSE v = v[, 1: ny]
[2]	405	END
	406	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
	407	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
	408	u = u[indice2d]
	409	v = v[indice2d]
	410	END
	411	ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
	412	endcase
	413	;------------------------------------------------------------------
[142]	414	; We reshape u and v to make sure that no dimension has been erased
[2]	415	;------------------------------------------------------------------
	416	if ny EQ 1 then begin
	417	u = reform(u, nx, ny, /over)
	418	v = reform(v, nx, ny, /over)
	419	endif
	420	;------------------------------------------------------------------
[142]	421	; construction of u and v at points T
[2]	422	;-----------------------------------------------------------
[114]	423	a=u[0,*]
[2]	424	u=(u+shift(u,1,0))/2.
[114]	425	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,*]=a
	426	a=v[*,0]
[2]	427	v=(v+shift(v,0,1))/2.
[114]	428	if NOT keyword_set(key_periodic) OR nx NE jpi then v[*,0]=a
[2]	429	;----------------------------------------------------------------------------
[142]	430	; attribution of the mask and of longitude and latitude arrays.
	431	; We recover the complete grid to establish a big mask extent in the four
	432	; direction to cover pointsfor which a land point has been
	433	; considerated (make a small drawing)
[2]	434	;----------------------------------------------------------------------------
[25]	435	mask = tmask[indice2d+jpijpjfirstzt]
[2]	436	if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
	437	;-----------------------------------------------------------
[142]	438	; construction of land containing all points to mask
[2]	439	;-----------------------------------------------------------
	440	if n_elements(valmask) EQ 0 THEN valmask = 1e20
	441	landu = where(u GE valmask/10)
	442	if landu[0] NE -1 then u[landu] = 0
	443	landv = where(v GE valmask/10)
	444	if landv[0] NE -1 then v[landv] = 0
	445	res=sqrt(u^2+v^2)
[114]	446	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,*]=!values.f_nan
	447	res[*,0]=!values.f_nan
[2]	448	mask = where(mask eq 0)
[114]	449	IF mask[0] NE -1 THEN res[mask] = valmask
[142]	450	; All kind of average
[25]	451	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
	452	if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
[2]	453	END
	454	;----------------------------------------------------------------------------
	455	endcase
	456	;------------------------------------------------------------
	457	if keyword_set(key_performance) THEN print, 'temps norme', systime(1)-tempsun
	458	return, res
	459	end

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

Download in other formats:

Original Format