Context Navigation

getturb.pro @ 117

Last change on this file since 117 was 113, checked in by aslmd, 14 years ago
LMD_LES_MARS: qq modifications pour les routines graphiques
File size: 46.9 KB

Rev	Line
[85]	1	pro getturb, saveps=saveps, noplot=noplot
	2
	3	;----------------------------------
	4	; USE: getturb
	5	; getturb, saveps='false'
	6	;----------------------------------
	7
	8	history_interval_s = 400.
	9	history_interval_s = 100.
	10	smoothampl=3700/history_interval_s
	11	;smoothampl=0. ;; no smooth
	12	;smoothampl=10.
	13	smoothampl=2.
	14	;smoothampl=5.
	15
[96]	16	;;
	17	;; constantes
	18	;;
[113]	19	p0=610. & t0=220. & r_cp=1.0/4.4 & grav=3.72 & R=192.
	20	;p0=610. & t0=220. & r_cp=1.0/3.9 & grav = 3.72 & R=191.
	21	;print, 'ATTENTION ATTENTION R/cp !!!!', r_cp
[85]	22
	23	;
	24	; graphics definition
	25	;
	26	if (n_elements(saveps) eq 0) then saveps='true'
	27	if (n_elements(noplot) eq 0) then noplot='false'
	28	if (saveps eq 'false') then begin
	29	;!p.multi=[0,3,2]
	30	!P.CHARSIZE=2.
	31	WINDOW, /PIXMAP & WDELETE & DEVICE,BYPASS_TRANSLATION=0,DECOMPOSED=0,RETAIN=2
	32	endif else begin
[96]	33	;PREF_SET, 'IDL_PATH', '/home/spiga/Save/SOURCES/IDL/fsc_psconfig:<IDL_DEFAULT>', /COMMIT
	34	PREF_SET, 'IDL_PATH', '/home/spiga/SVN/trunk/mesoscale/PLOT/MINIMAL:<IDL_DEFAULT>', /COMMIT
[85]	35	endelse
	36
	37	;
	38	; retrieve fields
	39	;
	40	openr,unit,'getturb.dat',/get_lun,error=err
	41	IF (err ne 0) THEN BEGIN
	42
	43	;
	44	; input files
	45	;
	46	OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22
	47	OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
	48
	49	;
	50	; get fields
	51	;
	52	domain='d01' & filesWRF = FindFile('wrfout_'+domain+'_????-??-??_??:??:??') & nf=n_elements(filesWRF)
	53
	54
	55	; get dimensions
	56	;
	57	id=ncdf_open(filesWRF(0))
	58	NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east' ), toto, nx & NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north' ), toto, ny
	59	NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top' ), toto, nz & NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, nt
	60	NCDF_CLOSE, id
	61	id=ncdf_open(filesWRF(nf-1)) ;; for interrupted runs
	62	NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, ntlast
	63	NCDF_CLOSE, id
	64
	65	;
	66	; prepare loop
	67	;
	68	;nloop1 = nf-1 & nloop2 = nt & yeye = 0
	69	yeye = 0 & nttot = (nf-1)*nt + ntlast
	70
	71	localtime = lctu + history_interval_s*findgen(nttot)/3700.
	72
	73	wt = fltarr(nz,nttot) & tke = fltarr(nz,nttot) & ztke = fltarr(nz,nttot) & t = fltarr(nz,nttot)
	74	p = fltarr(nz) & ph = fltarr(nz) & pht = fltarr(nz,nttot) & pt = fltarr(nz,nttot) & stst = fltarr(nz,nttot)
	75	xtke = fltarr(nz,nttot) & ytke = fltarr(nz,nttot) & wmax = fltarr(nz,nttot) & wmin = fltarr(nz,nttot)
	76	depressions = fltarr(nttot) & psmin = fltarr(nttot)
	77
	78	;
	79	; loop loop
	80	;
	81	for loop = 0, nf-1 do begin
	82	timetime = SYSTIME(1)
	83	if (loop ne nf-1) then nloop2=nt else nloop2=ntlast
	84
	85	for loop2 = 0, nloop2-1 do begin
	86
	87	anomalt = 1. & anomalu = 1. & anomalv = 1.
	88	;print, loop, loop2
	89
	90	; ------------
	91	; t' = t - <t>
	92	; ------------
	93	tprime = getget(filesWRF(loop), 'T', anomaly=anomalt, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; t' = t - <t>
	94	t(*,yeye) = t0 + TEMPORARY(anomalt)
	95	; ------
	96	; w' = w
	97	; ------
	98	wprime = getget(filesWRF(loop), 'W', count=[0,0,0,1], offset=[0,0,0,loop2])
	99	for toto = 0, nz-1 do begin
	100	wmax(toto,yeye) = max(reform(wprime(,,toto,0)),min=tutu) & wmin(toto,yeye) = tutu
	101	endfor
	102	; ; ------
	103	; ; u' = u and v' = v (car PAS de background wind !)
	104	; ; ------
	105	; uprime = getget(filesWRF(loop), 'U', count=[0,0,0,1], offset=[0,0,0,loop2])
	106	; vprime = getget(filesWRF(loop), 'V', count=[0,0,0,1], offset=[0,0,0,loop2])
	107	; --------------------------------------------------------
	108	; tke = 0.5 ( <u'^2> + <v'^2> + <w'^2> ) ; u' = u ; v' = v
	109	; --------------------------------------------------------
	110	; ztke is 0.5 * <w'^2>/2 or 0.5 * sigma_w^2
	111	ztke(,yeye) = 0.5 TOTAL(TOTAL(wprime^2,1),1) / float(nx) / float(ny)
	112	xtke(,yeye) = 0.5 TOTAL(TOTAL(getget(filesWRF(loop), 'U', anomaly=anomalu, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny)
	113	ytke(,yeye) = 0.5 TOTAL(TOTAL(getget(filesWRF(loop), 'V', anomaly=anomalv, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny)
	114	; xtke(,yeye) = 0.5 TOTAL(TOTAL(uprime^2,1),1) / float(nx) / float(ny)
	115	; ytke(,yeye) = 0.5 TOTAL(TOTAL(vprime^2,1),1) / float(nx) / float(ny)
	116	tke(,yeye) = xtke(,yeye) + $
	117	ytke(*,yeye) + $
	118	ztke(*,yeye)
	119	; ------
	120	; <w't'>
	121	; ------
	122	wt(,yeye) = TOTAL(TOTAL(TEMPORARY(tprime) TEMPORARY(wprime),1),1) / float(nx) / float(ny)
	123	; ------
	124	; p & ph
	125	; ------
	126	;if (loop + loop2 eq 0) then nloopbeware = nloop2-1 else nloopbeware = nloop2 ; 1ere valeur vaut 0
	127	nttotbeware = nttot-1 ; 1ere valeur vaut 0
	128	pht(*,yeye) = TOTAL(TOTAL(getget(filesWRF(loop), 'PHTOT', count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny) / 1000. / 3.72
	129	pt(*,yeye) = TOTAL(TOTAL(getget(filesWRF(loop), 'PTOT' , count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny)
	130	ph = TEMPORARY(ph) + pht(*,yeye) / nttotbeware ;/ nloopbeware / nloop1
	131	p = TEMPORARY(p ) + pt(*,yeye) / nttot ;/ nloop2 / nloop1
	132	; ----------------
	133	; static stability
	134	; ----------------
	135	stst(,yeye) = DERIV( reform(pht(,yeye)) - hgtu/1000. , reform(t(,yeye)) ( reform(pt(,yeye)) /p0 )^r_cp ) + 1000.grav / (R / r_cp) ;; 4.9 dans Hinson
	136	; ----------------
	137	; surface pressure
	138	; ----------------
	139	!QUIET=1
	140	;psfc = getget(filesWRF(loop), 'PSFC', anomaly=1, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; plus rapide que autre routine! pas grave avertissement
	141	psfc = getget(filesWRF(loop), 'PSFC', count=[0,0,0,1], offset=[0,0,0,loop2])
	142	!QUIET=0
	143	; psfc = getget2d(filesWRF(loop), 'PSFC', anomaly=1, count=[0,0,1], offset=[0,0,loop2])
	144	psmin(yeye) = min(psfc) ;& print, psmin(yeye)
	145	w = where(TEMPORARY(psfc) lt -0.5) ;& print, n_elements(w)
	146	if (w(0) ne -1) then depressions(yeye) = 1000. * n_elements(w) / float(nx) / float(ny) else depressions(yeye) = 0.
	147	; ----------
	148	; loop count
	149	; ----------
	150	yeye = TEMPORARY(yeye) + 1 ;& print, yeye & print, SYSTIME(1) - timetime, ' s'
	151
	152	endfor
	153	print, 'file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s'
	154	endfor
	155	h = TEMPORARY(ph) - hgtu/1000. ;; altitude above ground
	156	ht = TEMPORARY(pht) - hgtu/1000.
	157	;
	158	; save
	159	;
	160	save, wt, tke, ztke, h, ht, t, p, pt, stst, localtime, xtke, ytke, wmax, wmin, depressions, psmin, filename='getturb.dat'
	161	nz = n_elements(h)
	162	nnn = n_elements(tke(0,*))
	163
	164	ENDIF ELSE BEGIN
	165
	166	print, 'OK, file is here'
	167	restore, filename='getturb.dat'
	168	nz = n_elements(h)
	169	nnn = n_elements(tke(0,*))
	170
	171	OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
	172
[113]	173	;@include_ustar.pro
[85]	174
	175	;@include_hfx.pro
	176
	177	;@include_w.pro
	178
[96]	179	;@include_tprime.pro
	180
[85]	181	ENDELSE
	182
	183	;
	184	; pbl height
	185	;
	186	pbl_height = fltarr(nnn) & hfsurf = fltarr(nnn) & w_star = fltarr(nnn) & hfbot = fltarr(nnn) & w_star_bot = fltarr(nnn) & hfbotmax = fltarr(nnn)
	187	a_vel = fltarr(nz,nnn) & a_h = fltarr(nz,nnn) & a_wt = fltarr(nz,nnn) & a_vel_bot = fltarr(nz,nnn) & a_wt_bot = fltarr(nz,nnn)
	188	tt_top = fltarr(nnn) & z_top = fltarr(nnn) & t_bot = fltarr(nnn) & p_bot = fltarr(nnn)
	189	tt_bot = fltarr(nnn)
	190	wmaxmax = fltarr(nnn)
	191	wminmin = fltarr(nnn)
	192	tkemax = fltarr(nnn)
	193
	194	for i=1, nnn-1 do begin ;; ne pas traiter l'init
	195
	196	;
	197	; pbl height
	198	;
	199	z_in = reform(ht(*,i))
	200	profile_in = reform(stst(*,i)) ;; wt pas suffisamment regulier
	201	resol = 1000 ;200 un peu juste
	202	;********************************************************************
	203	; ---- numerical recipes in C - section 3.3 ----
	204	;
	205	; Calculate interpolating cubic spline
	206	yspline = SPL_INIT(z_in,profile_in)
	207	; Define the X values P at which we desire interpolated Y values
	208	xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol
	209	; Calculate the interpolated Y values
	210	result=spl_interp(z_in,profile_in,yspline,xspline)
	211	;********************************************************************
	212	w = where( (result ge 1.5) and (xspline gt 0.01) ) ;0.3) ) ;0.1 suffit si HR
	213	;if ( localtime(i) gt 16. ) then w = where( (result ge 1.5) and (xspline gt 2.) )
	214
	215	if ( localtime(i) gt 15. ) then w = where( (result ge 1.5) and (xspline gt 1.) )
	216
	217	;;;; special tau = 5.
	218	;if ( localtime(i) ge 13. ) then w = where( (result ge 1.5) and (xspline gt 0.35) )
	219	;if ( localtime(i) ge 16.05 ) then w = [0]
	220
	221	pbl_height(i) = xspline(w(0)) ;& print, 'PBL height - mean profile', pbl_height(i)
	222	z_top(i) = xspline(w(0)) + hgtu/1000.
	223
	224	;
	225	; temperature @ pbl height
	226	;
	227	z_in = reform(ht(*,i))
	228	profile_in = reform(t(*,i))
	229	;profile_in = reform(t(,i))(reform(pt(*,i))/p0)^r_cp
	230	profile_in2 = reform(pt(*,i))
	231	resol = 1000 ;200 un peu juste
	232	;********************************************************************
	233	; ---- numerical recipes in C - section 3.3 ----
	234	;
	235	; Calculate interpolating cubic spline
	236	yspline = SPL_INIT(z_in,profile_in)
	237	yspline2 = SPL_INIT(z_in,profile_in2)
	238	; Define the X values P at which we desire interpolated Y values
	239	xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol
	240	; Calculate the interpolated Y values
	241	result=spl_interp(z_in,profile_in,yspline,xspline)
	242	result2=spl_interp(z_in,profile_in2,yspline2,xspline)
	243	;********************************************************************
	244	w = where( abs(xspline-pbl_height(i)) eq min(abs(xspline-pbl_height(i))) ) & tt_top(i) = result[w(0)]*(result2[w(0)]/p0)^r_cp ;& print, z_top(i), t_top(i)
	245	hhh = 1000. / 1000. ;100. / 1000.
	246	w = where( abs(xspline-hhh) eq min(abs(xspline-hhh)) ) & t_bot(i) = result[w(0)] & p_bot(i) = result2[w(0)] ;& print, p_bot(i), t_bot(i)
	247	tt_bot(i) = result[w(0)]*(result2[w(0)]/p0)^r_cp
	248
	249	;
	250	; heat flux at the "surface"
	251	;
	252	hfsurf(i) = wt(1,i) ;; attention vaut 0 a la hauteur 0
	253
	254	;
	255	; heat flux at the bottom of mixed layer = top of radiative layer
	256	;
	257	z_in = reform(ht(*,i))
	258	profile_in = reform(wt(*,i)) ;; wt pas suffisamment regulier
	259	resol = 1000 ;200 un peu juste
	260	;********************************************************************
	261	; ---- numerical recipes in C - section 3.3 ----
	262	;
	263	; Calculate interpolating cubic spline
	264	yspline = SPL_INIT(z_in,profile_in)
	265	; Define the X values P at which we desire interpolated Y values
	266	xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol
	267	; Calculate the interpolated Y values
	268	result=spl_interp(z_in,profile_in,yspline,xspline)
	269	;********************************************************************
	270	hhh = 300./1000. & yeah=where(abs(xspline-hhh) eq (min(abs(xspline-hhh)))) & nnzz=yeah(0) & hfbot(i)=result(nnzz) ;; 300m semble OK, 500m super
	271	w=where(xspline lt 1.5) & result=result[w] & xspline=xspline[w] & yeah=where(result eq max(result)) & nnzz=yeah(0) & hfbotmax(i)=result(nnzz)
	272	;print, 'surf, 300m, max ', hfsurf(i), hfbot(i), hfbotmax(i) ;& plot, result, xspline
	273
	274	result = reform(wmax(*,i)) & yeah=where(result eq max(result)) & nnzz=yeah(0) & wmaxmax(i)=result(nnzz)
	275	result = reform(wmin(*,i)) & yeah=where(result eq min(result)) & nnzz=yeah(0) & wminmin(i)=result(nnzz)
	276	result = reform(tke(*,i)) & yeah=where(result eq max(result)) & nnzz=yeah(0) & tkemax(i )=result(nnzz)
	277
	278	;
	279	; free convection velocity scale
	280	;
	281	; w_star_bot(i) = ( 1000.* grav * pbl_height(i) * hfbot(i) / t(0,i) ) ^ (1./3.)
	282	w_star_bot(i) = ( 1000.* grav * pbl_height(i) * hfbotmax(i) / t(0,i) ) ^ (1./3.)
	283	w_star(i) = ( 1000.* grav * pbl_height(i) * hfsurf(i) / t(0,i) ) ^ (1./3.)
	284
	285	;
	286	; dimensional quantities
	287	;
	288	a_vel(,i) = reform( 2.ztke(*,i) / w_star(i)^2 )
	289	a_vel_bot(,i) = reform( 2.ztke(*,i) / w_star_bot(i)^2 )
	290	a_h(,i) = reform( ht(,i) / pbl_height(i) )
	291	a_wt(,i) = reform( wt(,i) / hfsurf(i) ) ;reform( wt(*,i) / w_star(i) )
	292	a_wt_bot(,i) = reform( wt(,i) / hfbotmax(i) ) ;reform( wt(*,i) / w_star_bot(i) )
	293
	294	endfor
	295	;
	296	; free convection time scale
	297	;
	298	fcts = 1000. * pbl_height / w_star / 3700.
	299	;
	300	; mixed layer temperature scale
	301	;
	302	mlts = hfsurf / w_star
	303
	304	;pbl_height = SMOOTH(TEMPORARY(pbl_height), [smoothampl], /EDGE_TRUNCATE)
	305	;hfsurf = SMOOTH(TEMPORARY(hfsurf), [smoothampl], /EDGE_TRUNCATE)
	306	;w_star = SMOOTH(TEMPORARY(w_star), [smoothampl], /EDGE_TRUNCATE)
	307	;hfbot = SMOOTH(TEMPORARY(hfbot), [smoothampl], /EDGE_TRUNCATE)
	308	;w_star_bot = SMOOTH(TEMPORARY(w_star_bot), [smoothampl], /EDGE_TRUNCATE)
	309	;hfbotmax = SMOOTH(TEMPORARY(hfbotmax), [smoothampl], /EDGE_TRUNCATE)
	310	;a_vel = SMOOTH(TEMPORARY(a_vel), [0,smoothampl], /EDGE_TRUNCATE)
	311	;a_h = SMOOTH(TEMPORARY(a_h), [0,smoothampl], /EDGE_TRUNCATE)
	312	;a_wt = SMOOTH(TEMPORARY(a_wt), [0,smoothampl], /EDGE_TRUNCATE)
	313	;a_vel_bot = SMOOTH(TEMPORARY(a_vel_bot), [0,smoothampl], /EDGE_TRUNCATE)
	314	;a_wt_bot = SMOOTH(TEMPORARY(a_wt_bot), [0,smoothampl], /EDGE_TRUNCATE)
	315	;tt_top = SMOOTH(TEMPORARY(tt_top), [smoothampl], /EDGE_TRUNCATE)
	316	;z_top = SMOOTH(TEMPORARY(z_top), [smoothampl], /EDGE_TRUNCATE)
	317	;t_bot = SMOOTH(TEMPORARY(t_bot), [smoothampl], /EDGE_TRUNCATE)
	318	;p_bot = SMOOTH(TEMPORARY(p_bot), [smoothampl], /EDGE_TRUNCATE)
	319
	320	;
	321	; save for comparison sake
	322	;
	323	save, pbl_height, hfbot, hfsurf, w_star, w_star_bot, a_vel_bot, a_vel, a_h, a_wt_bot, a_wt, fcts, mlts, hfbotmax, tt_top, z_top, t_bot, p_bot, tt_bot, wmaxmax, wminmin, tkemax, filename='compturb.dat'
	324	if (noplot ne 'false') then stop
	325
	326	;
	327	; smooth smooth
	328	;
	329	wt = SMOOTH(TEMPORARY(wt), [0,smoothampl], /EDGE_TRUNCATE)
	330	tke = SMOOTH(TEMPORARY(tke), [0,smoothampl], /EDGE_TRUNCATE)
	331	ztke = SMOOTH(TEMPORARY(ztke), [0,smoothampl], /EDGE_TRUNCATE)
	332	xtke = SMOOTH(TEMPORARY(xtke), [0,smoothampl], /EDGE_TRUNCATE)
	333	ytke = SMOOTH(TEMPORARY(ytke), [0,smoothampl], /EDGE_TRUNCATE)
	334	wmax = SMOOTH(TEMPORARY(wmax), [0,smoothampl], /EDGE_TRUNCATE)
	335	wmin = SMOOTH(TEMPORARY(wmin), [0,smoothampl], /EDGE_TRUNCATE)
	336	depressions = SMOOTH(TEMPORARY(depressions), [smoothampl], /EDGE_TRUNCATE)
	337	psmin = SMOOTH(TEMPORARY(psmin), [smoothampl], /EDGE_TRUNCATE)
	338
	339
	340	;*******************;
	341	;*******************;
	342	; PLOTS PLOTS PLOTS ;
	343	;*******************;
	344	;*******************;
	345
[113]	346	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	347	;;zefield = a_vel_bot
	348	;;zey = a_h
	349	;;set_name = 'A_VEL.ps'
	350	;;set_title = "Scaled vertical velocity variance (m.s!U-1!N)"
	351	;;set_titlex = 'Scaled vertical velocity variance (m.s!U-1!N)'
	352	;;set_titley = 'Scaled height (km)'
	353	;;set_subtitle = ''
	354	;;set_xrange = [0.,1.2]
	355	;;set_yrange = [0.,1.2]
	356	;;set_tickx = 0.1
	357	;;set_ticky = 0.1
	358	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	359	;;zesim = 1.8 * a_h^(2./3.) * ( 1. - 0.8 * a_h )^2
	360	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	361	;;localtimes = localtime[where( (localtime ge 11) and (localtime le 16) )]
	362	;;;localtimes = localtime
	363	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	364	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	365	;;!P.Charsize = 1.2
	366	;;altlin=0
	367	;;user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0)
	368	;;plot, zefield(,nntt), zey(,nntt), xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin, /ISOTROPIC
	369	;;for ll = 1, n_elements(localtimes)-1 do begin
	370	;; user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0)
	371	;; if (nntt ne -1) then oplot, zefield(,nntt), zey(,nntt), linestyle=altlin
	372	;;endfor
	373	;;loadct, 33 & oplot, zesim, zey, psym=5, color=255 & loadct, 0
	374	;;if (saveps eq 'true') then PS_End, /PNG
	375	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	376	;
[85]	377	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	378	;zefield = localtime
[113]	379	;zey = pbl_height
	380	;set_name = 'HEIGHT.ps'
	381	;set_title = 'Boundary layer height (km)'
[85]	382	;set_titlex = 'Local Time (h)'
[113]	383	;set_titley = 'Boundary layer height (km)'
	384	;set_subtitle = '' ;'Criterion is : static stability > 1.5 K.m!U-1!N'
	385	;set_xrange = [11.,17.] ;[8.,17.]
	386	;set_yrange = [0.,8.]
[85]	387	;set_tickx = 1.
[113]	388	;set_ticky = 1.
[85]	389	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	390	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	391	;!P.Charsize = 1.2
	392	;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	393	;;oplot, zefield, zey, psym=5
	394	;if (saveps eq 'true') then PS_End, /PNG
	395	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	396	;
	397	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	398	;;zefield = localtime
	399	;;zey = w_star_bot
	400	;;set_name = 'W_STAR.ps'
	401	;;set_title = "Free convection velocity scale (m.s!U-1!N)"
	402	;;set_titlex = 'Local Time (h)'
	403	;;set_titley = 'Free convection velocity scale (m.s!U-1!N)'
	404	;;set_subtitle = ''
	405	;;set_xrange = [8.,18.]
	406	;;set_yrange = [0.,6.]
	407	;;set_tickx = 1.
	408	;;set_ticky = 0.5
	409	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	410	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	411	;;!P.Charsize = 1.2
	412	;;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	413	;;;oplot, zefield, zey, psym=5
	414	;;if (saveps eq 'true') then PS_End, /PNG
	415	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	416	;
	417	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	418	;;zefield = localtime
	419	;;zey = hfbot
	420	;;zey2 = hfsurf
	421	;;zey3 = hfbotmax
	422	;;set_name = 'HFBOT.ps'
	423	;;set_title = "Heat flux at the top of the radiative layer (K.m.s!U-1!N)"
	424	;;set_titlex = 'Local Time (h)'
	425	;;set_titley = 'Heat flux at the top of the radiative layer (K.m.s!U-1!N)'
	426	;;set_subtitle = ''
	427	;;set_xrange = [8.,18.]
	428	;;set_yrange = [0.,2.5]
	429	;;set_tickx = 1.
	430	;;set_ticky = 0.5
	431	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	432	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	433	;;!P.Charsize = 1.2
	434	;;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	435	;;oplot, zefield, zey2, linestyle=1
	436	;;oplot, zefield, zey3, linestyle=2
	437	;;if (saveps eq 'true') then PS_End, /PNG
	438	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	439	;
[85]	440	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	441	;zefield = transpose(tke)
	442	;zex = localtime
	443	;zey = h
	444	;set_name = 'TKE.ps'
	445	;set_title = "Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<u'!U2!N>+<v'!U2!N>+<w'!U2!N>]
[85]	446	;set_titlex = 'Local Time (h)'
[113]	447	;set_titley = 'Altitude above surface (km)'
	448	;set_subtitle = '' ;'Mean over the simulation domain'
	449	;set_xrange = [11.,17.] ;[8.,17.]
	450	;set_yrange = [0.,8.] ;; [0.,200.]
[85]	451	;set_tickx = 1.
[113]	452	;set_ticky = 1. ;; 50.
	453	;minval = 0.
	454	;maxval = 20.
	455	;nlev = maxval-minval
	456	;pal = 22
	457	;rrr = 'no'
[85]	458	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	459	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	460	;!P.Charsize = 1.2
[113]	461	;;; 0. levels
	462	;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	463	;;;; 1. background
	464	;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	465	; ;; 2. color field
	466	; loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	467	; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	468	;;; 3. contour field
	469	;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	470	;;; 4. choose output
[85]	471	;if (saveps eq 'true') then PS_End, /PNG
	472	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	473	;
[85]	474	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	475	;;zefield = 2.*transpose(ztke)
	476	;;set_title = "Vertical wind variance (m!U2!N.s!U-2!N)"
	477	;;minval = -6.
	478	;;maxval = 12.
	479	;;pal = 0
	480	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	481	;zefield = transpose(ztke)
[85]	482	;zex = localtime
	483	;zey = h
[113]	484	;set_name = 'zTKE.ps'
	485	;set_title = "Vertical Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<w'!U2!N>]
[85]	486	;set_titlex = 'Local Time (h)'
	487	;set_titley = 'Altitude above surface (km)'
	488	;set_subtitle = '' ;'Mean over the simulation domain'
[113]	489	;set_xrange = [8.,17.]
[85]	490	;set_yrange = [0.,7.]
	491	;set_tickx = 1.
	492	;set_ticky = 1.
	493	;minval = 0.
[113]	494	;maxval = 8.
[85]	495	;nlev = maxval-minval
	496	;pal = 22
	497	;rrr = 'no'
	498	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	499	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	500	;!P.Charsize = 1.2
	501	;;; 0. levels
	502	;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	503	;;;; 1. background
	504	;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	505	;;; 2. color field
	506	;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	507	; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	508	;;; 3. contour field
	509	;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	510	;;; 4. choose output
	511	;if (saveps eq 'true') then PS_End, /PNG
	512	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	513	;
	514	;;restore, filename='addturb.dat'
	515	;;modvar = SMOOTH(TEMPORARY(modvar), [0,smoothampl], /EDGE_TRUNCATE)
	516	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	517	;;zefield = transpose(modvar)
	518	;;zex = localtime
	519	;;zey = h
	520	;;set_name = 'modvar.ps'
	521	;;set_title = "Horizontal wind speed variance (m!U2!N.s!U-2!N)"
	522	;;set_titlex = 'Local Time (h)'
	523	;;set_titley = 'Altitude above surface (km)'
	524	;;set_subtitle = '' ;'Mean over the simulation domain'
	525	;;set_xrange = [7.,17.]
	526	;;set_yrange = [0.,7.]
	527	;;set_tickx = 1.
	528	;;set_ticky = 1.
	529	;;minval = 0.
	530	;;maxval = 10.
	531	;;nlev = maxval-minval
	532	;;pal = 22
	533	;;rrr = 'no'
	534	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	535	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	536	;;!P.Charsize = 1.2
	537	;;;; 0. levels
	538	;;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	539	;;;;; 1. background
	540	;;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	541	;;;; 2. color field
	542	;;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	543	;; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	544	;;;; 3. contour field
	545	;;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	546	;;;; 4. choose output
	547	;;if (saveps eq 'true') then PS_End, /PNG
	548	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	549	;
	550	;
[85]	551	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	552	;zefield = transpose(wt)
	553	;zex = localtime
	554	;zey = h
	555	;set_name = 'HF.ps'
	556	;set_title = "Turbulent Heat Flux (K.m.s!U-1!N)" ;"Vertical Eddy Heat Flux <w'T'>" ;<w'!7h!3'>
	557	;set_titlex = 'Local Time (h)'
[85]	558	;set_titley = 'Altitude above surface (km)'
	559	;set_subtitle = '' ;'Mean over the simulation domain'
[113]	560	;set_xrange = [8.,17.]
	561	;set_yrange = [0.,8.] ;; [0.,200.]
	562	;set_tickx = 1.
	563	;set_ticky = 1. ;; 50.
	564	;minval = -1.5 ;-2.
	565	;maxval = 2.5 ;2.
	566	;nlev = floor(maxval-minval)*10
	567	;pal = 33
	568	;rrr = 'no'
[85]	569	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	570	;;minval = -0.8
	571	;;maxval = 1.2
	572	;;pal = 0
[85]	573	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	574	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	575	;!P.Charsize = 1.2
[113]	576	;;; 0. levels
	577	;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	578	;;;; 1. background
	579	;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	580	;;; 2. color field
	581	;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	582	; ;;;--------------------------------------------------------------------------------------------------------------------------------
	583	; ;;; WHITE ZONE - 1. get location of interval in the CT - 2. change the CT to have a white zone
	584	; ulim=0.09 & dlim=-0.09 & w=where(lev le dlim) & n1=w[n_elements(w)-1] & w=where(lev ge ulim) & n2=w[0] & yy=BYTSCL(lev) & nd=yy[n1] & nu=yy[n2]-5
	585	; nu = nd + (nu-nd)/2 ;; otherwise the interval is too large (because we removed 0)
	586	; TVLCT, r, g, b, /Get & r[nd:nu]=255 & g[nd:nu]=255 & b[nd:nu]=255 & TVLCT, r, g, b
	587	; ;;;--------------------------------------------------------------------------------------------------------------------------------
	588	; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	589	;;; 3. contour field
	590	;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	591	;;; 4. choose output
[85]	592	;if (saveps eq 'true') then PS_End, /PNG
	593	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	594	;
	595	;;;restore, filename='tpot_profB'
[85]	596	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	597	;;zefield = t
	598	;;zey = ht
	599	;;set_name = 'T.ps'
	600	;;set_title = "Potential Temperature (K)"
	601	;;set_titlex = 'Potential Temperature (K)'
	602	;;set_titley = 'Altitude above surface (km)'
	603	;;set_subtitle = '' ;'Mean over the simulation domain'
	604	;;set_xrange = [190.,230.] ;[min(t),max(t)]
	605	;;set_yrange = [0.,7.]
	606	;;set_tickx = 5.
	607	;;set_ticky = 1.
[85]	608	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	609	;;;localtimes = [9,10,11,12,13,14,15,16,17,18]
	610	;;localtimes = [7,9,11,13,15,17]
	611	;;;localtimes = [17.0,17.1,17.2,17.3,17.4,17.5,17.6,17.7,17.8,17.9,18.0]
	612	;;;localtimes = [15.0]
[85]	613	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	614	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	615	;;!P.Charsize = 1.2
	616	;;;altlin=0 & loadct, 0
	617	;;altlin=1 & loadct, 0
	618	;;user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0)
	619	;;plot, zefield(,nntt), zey(,nntt), xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	620	;;;plot, zefield(*,nntt), les_column, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	621	;;for ll = 1, n_elements(localtimes)-1 do begin
	622	;; CASE altlin OF
	623	;; 0: altlin=1
	624	;; 1: altlin=0
	625	;; ENDCASE
	626	;; user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0)
	627	;; if (nntt ne -1) then oplot, zefield(,nntt), zey(,nntt), linestyle=altlin
	628	;;; if (nntt ne -1) then oplot, zefield(*,nntt), les_column, linestyle=altlin
	629	;;endfor
	630	;;;oplot, ro_tpot, ro_column, psym=7
	631	;;;xyouts, 192.0, 0.25, '09:00'
	632	;;;xyouts, 205.5, 0.25, '11:00'
	633	;;;xyouts, 214.5, 0.25, '13:00'
	634	;;;xyouts, 219.5, 0.25, '17:00'
	635	;;;xyouts, 220.5, 2.30, '17:00 [RO]'
	636	;;if (saveps eq 'true') then PS_End, /PNG
	637	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	638	;
	639	;
	640	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	641	;;zefield = stst
	642	;;zey = ht
	643	;;set_name = 'STST.ps'
	644	;;set_title = "Static stability (K.m!U-1!N)"
	645	;;set_titlex = 'Static stability (K.m!U-1!N)'
	646	;;set_titley = 'Altitude above surface (km)'
	647	;;set_subtitle = 'Mean over the simulation domain'
	648	;;set_xrange = [-5.,5.]
	649	;;set_yrange = [0.,7.]
	650	;;set_tickx = 1.
	651	;;set_ticky = 1.
	652	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	653	;;localtimes = [9,11,13,15,17]
	654	;;localtimes = [12,13,14,15,16]
	655	;;;localtimes = [17.0,17.1,17.2,17.3,17.4,17.5,17.6,17.7,17.8,17.9,18.0]
	656	;;;localtimes = [15.0]
	657	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	658	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	659	;;!P.Charsize = 1.2
	660	;;altlin=0 & loadct, 0
	661	;;user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0)
	662	;;plot, zefield(,nntt), zey(,nntt), xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	663	;;oplot, zefield(,nntt), zey(,nntt), psym=5
	664	;;for ll = 1, n_elements(localtimes)-1 do begin
	665	;; CASE altlin OF
	666	;; 0: altlin=1
	667	;; 1: altlin=0
	668	;; ENDCASE
	669	;; user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0)
	670	;; if (nntt ne -1) then oplot, zefield(,nntt), zey(,nntt), linestyle=altlin
	671	;; if (nntt ne -1) then oplot, zefield(,nntt), zey(,nntt), psym=5
	672	;;endfor
	673	;;oplot, 0.zefield(,nntt) + 1.5, zey(*,nntt), linestyle=2
	674	;;if (saveps eq 'true') then PS_End, /PNG
	675	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	676	;
	677	;
	678	;
	679	;if (n_elements(xtke) eq 0) then stop
	680	;
	681	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	682	;;;zefield = 2.*transpose(xtke)
	683	;;;minval = -4
	684	;;;maxval = 8. ;12.
	685	;;;pal = 0
	686	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	687	;;zex = localtime
	688	;;zey = h
	689	;;set_name = 'xTKE.ps'
	690	;;set_title = "Horizontal Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<u'!U2!N>]
	691	;;set_titlex = 'Local Time (h)'
	692	;;set_titley = 'Altitude above surface (km)'
	693	;;set_subtitle = '';'Mean over the simulation domain'
	694	;;set_xrange = [8.,17.]
	695	;;set_yrange = [0.,7.]
	696	;;set_tickx = 1.
	697	;;set_ticky = 1.
	698	;;minval = 0.
	699	;;maxval = 8.
	700	;;nlev = maxval-minval
	701	;;pal = 22
	702	;;rrr = 'no'
	703	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	704	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	705	;;!P.Charsize = 1.2
	706	;;;; 0. levels
	707	;;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	708	;;;;; 1. background
	709	;;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	710	;;;; 2. color field
	711	;;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	712	;; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	713	;;;; 3. contour field
	714	;;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	715	;;;; 4. choose output
	716	;;if (saveps eq 'true') then PS_End, /PNG
	717	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	718	;
	719	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	720	;zefield = transpose(ytke) + transpose(xtke)
[85]	721	;zex = localtime
	722	;zey = h
[113]	723	;set_name = 'hTKE.ps'
	724	;set_title = "Horizontal Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<v'!U2!N>]
[85]	725	;set_titlex = 'Local Time (h)'
	726	;set_titley = 'Altitude above surface (km)'
[113]	727	;set_subtitle = '' ;'Mean over the simulation domain'
[85]	728	;set_xrange = [8.,17.]
[113]	729	;set_yrange = [0.,1.] ;[0.,7.] ;; [0.,200.]
[85]	730	;set_tickx = 1.
[113]	731	;set_ticky = 0.1 ;1. ;50.
[85]	732	;minval = 0.
[113]	733	;maxval = 12. ;10.
[85]	734	;nlev = maxval-minval
	735	;pal = 22
	736	;rrr = 'no'
	737	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	738	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	739	;!P.Charsize = 1.2
	740	;;; 0. levels
	741	;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	742	;;;; 1. background
	743	;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	744	;;; 2. color field
	745	;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	746	; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	747	;;; 3. contour field
	748	;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	749	;;; 4. choose output
	750	;if (saveps eq 'true') then PS_End, /PNG
	751	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	752	;
	753	;if (n_elements(wmax) eq 0) then stop
[85]	754
	755	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	756	zefield = transpose(wmax)
	757	zex = localtime
	758	zey = h
	759	set_name = 'Wmax.ps'
	760	set_title = "Maximum updraft speed (m.s!U-1!N)"
	761	set_titlex = 'Local Time (h)'
	762	set_titley = 'Altitude above surface (km)'
	763	set_subtitle = ''
[113]	764	set_xrange = [11.,17.] ;[8.,18.]
	765	set_yrange = [0.,8.] ;[0.,8.]
[85]	766	set_tickx = 1.
	767	set_ticky = 1.
	768	minval = 0.
[113]	769	maxval = 18. ;16. ;15.;14.;13.
[85]	770	nlev = maxval-minval
	771	pal = 22
	772	rrr = 'no'
	773	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	774	if (saveps eq 'true') then PS_Start, FILENAME=set_name
	775	!P.Charsize = 1.2
	776	;; 0. levels
	777	lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
[113]	778	levhalf = minval + (maxval-minval)*findgen((nlev/2)+1)/float(nlev/2) & if (minval ne 0.) then levhalf = levhalf[where(levhalf ne 0.)]
[85]	779	;;; 1. background
	780	loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	781	;; 2. color field
	782	loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	783	contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	784	;; 3. contour field
[113]	785	loadct, 0 & contour, zefield, zex, zey, levels=levhalf, c_labels=findgen(n_elements(levhalf))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (levhalf LT 0.0)
[85]	786	;; 4. choose output
	787	if (saveps eq 'true') then PS_End, /PNG
	788	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	789
	790	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	791	zefield = abs(transpose(wmin))
	792	zex = localtime
	793	zey = h
	794	set_name = 'Wmin.ps'
	795	set_title = "Maximum downdraft speed (m.s!U-1!N)"
	796	set_titlex = 'Local Time (h)'
	797	set_titley = 'Altitude above surface (km)'
	798	set_subtitle = ''
[113]	799	set_xrange = [11.,17.] ;[8.,18.]
	800	set_yrange = [0.,8.]
[85]	801	set_tickx = 1.
	802	set_ticky = 1.
	803	minval = 0.
[113]	804	maxval = 10.;12.
[85]	805	nlev = maxval-minval
	806	pal = 22
	807	rrr = 'no' ;'yes'
	808	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	809	if (saveps eq 'true') then PS_Start, FILENAME=set_name
	810	!P.Charsize = 1.2
	811	;; 0. levels
	812	lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
[113]	813	levhalf = minval + (maxval-minval)*findgen((nlev/2)+1)/float(nlev/2) & if (minval ne 0.) then levhalf = levhalf[where(levhalf ne 0.)]
[85]	814	;;; 1. background
	815	loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	816	;; 2. color field
	817	loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	818	contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	819	;; 3. contour field
[113]	820	loadct, 0 & contour, zefield, zex, zey, levels=levhalf, c_labels=findgen(n_elements(levhalf))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (levhalf LT 0.0)
[85]	821	;; 4. choose output
	822	if (saveps eq 'true') then PS_End, /PNG
	823	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	824	;
	825	;
	826	;restore, filename='addturb.dat'
	827	;velmax = SMOOTH(TEMPORARY(velmax), [0,smoothampl], /EDGE_TRUNCATE)
[85]	828	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	829	;zefield = transpose(velmax)
	830	;zex = localtime
	831	;zey = h
	832	;set_name = 'velmax.ps'
	833	;set_title = "Maximum horizontal wind speed (m.s!U-1!N)"
[85]	834	;set_titlex = 'Local Time (h)'
[113]	835	;set_titley = 'Altitude above surface (km)'
	836	;set_subtitle = '' ;'Mean over the simulation domain'
	837	;set_xrange = [11.,17.] ;[8.,17.]
	838	;set_yrange = [0.,4.]
[85]	839	;set_tickx = 1.
	840	;set_ticky = 1.
[113]	841	;minval = 0.
	842	;maxval = 30. ;12.
	843	;nlev = maxval-minval
	844	;pal = 22
	845	;rrr = 'no'
[85]	846	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	847	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	848	;!P.Charsize = 1.2
[113]	849	;;; 0. levels
	850	;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
	851	;;;; 1. background
	852	;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0
	853	;;;; 2. color field
	854	;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
	855	; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
	856	;;; 3. contour field
	857	;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0)
	858	;;; 4. choose output
[85]	859	;if (saveps eq 'true') then PS_End, /PNG
	860	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[113]	861	;
	862	;
	863	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	864	;;zefield = localtime
	865	;;zey = depressions
	866	;;set_name = 'DEPR.ps'
	867	;;set_title = "Percentage of depressions in the domain"
	868	;;set_titlex = 'Local Time (h)'
	869	;;set_titley = 'Percentage'
	870	;;set_subtitle = ''
	871	;;set_xrange = [8.,18.]
	872	;;set_yrange = [0.,3.]
	873	;;set_tickx = 1.
	874	;;set_ticky = 1.
	875	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	876	;;if (saveps eq 'true') then PS_Start, FILENAME=set_name
	877	;;!P.Charsize = 1.2
	878	;;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin
	879	;;oplot, zefield, zey, psym=5
	880	;; oplot, localtime, abs(psmin), linestyle=1
	881	;;if (saveps eq 'true') then PS_End, /PNG
	882	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[85]	883
	884	stop
	885
	886
	887
	888
	889	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	890	goto, no_staticstab
	891	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	892	if (saveps eq 'true') then begin
	893	device, /close
	894	set_plot, 'ps' & device, filename='plot/staticstab.ps'
	895	endif
	896
	897	user_lt=17. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
	898
	899	;;; recompute height @ given local time
	900	caca=heightp+hgtu/1000.
	901	height=reform(ph(,0,,*))
	902	height=total(height,1)/n_elements(height(*,0,0))
	903	height=reform(height)
	904	height=reform(height(*,yeah(0)))
	905	height=height/1000./3.72
	906	heightp=height(0:n_elements(height(*))-2)
	907	print, 'new minus old', heightp - caca
	908	;;; recompute height @ given local time
	909
	910	press=reform(p(,0,,*))
	911	press=total(press,1)/n_elements(press(*,0,0))
	912	press=reform(press)
	913	press=reform(press(*,yeah(0)))
	914	press=press(0:n_elements(press(*))-2)
	915
	916	staticstab = DERIV( heightp , reform(what_I_plot5(,yeah(0))) ) + 1000.3.72/844.6 ;; 4.9 dans Hinson
	917	staticstab = staticstab(0:n_elements(staticstab)-5)
	918	heightp = heightp(0:n_elements(heightp)-5)
	919
	920	plot, $
	921	staticstab, $
	922	heightp,$
	923	xtitle='Static stability (K/km)',$
	924	xrange=[-1.,5.], $
	925	xtickinterval=0.5, $
	926	ytitle='Altitude above surface (km)', $
	927	yrange=[min(heightp),max(heightp)], $
	928	ytickinterval=1., $
	929	title="LMD LES Static stability @ LT 17h (K/km)", $
	930	subtitle='zonal average at lat. '+latwrite
	931	oplot, $
	932	staticstab, $
	933	heightp,$
	934	psym=5
	935	oplot, $
	936	findgen(n_elements(heightp))*0. + 1.,$
	937	heightp,$
	938	linestyle=2
	939	oplot, $
	940	findgen(n_elements(heightp))*0. + 2.,$
	941	heightp,$
	942	linestyle=2
	943
	944	;;;;;;;;;;
	945	w = where((staticstab gt 1.5) and ((heightp- hgtu/1000.) gt 1.))
	946	t_top_plus = what_I_plot5(w(0),yeah(0))
	947	z_top_plus = heightp(w(0))
	948	p_top_plus = press(w(0))
	949	t_top_moins = what_I_plot5(w(0)-1,yeah(0))
	950	z_top_moins = heightp(w(0)-1)
	951	p_top_moins = press(w(0)-1)
	952	pbl_depth = (z_top_plusalog(p_top_plus) + z_top_moinsalog(p_top_moins))/(alog(p_top_plus) + alog(p_top_moins)) - hgtu/1000.
	953	xyouts, 3., 1.5 + (max(heightp) + min(heightp)) / 3., 'Ls = '+string(lsu,'(F5.1)')+'!Uo!N', CHARSIZE=1
	954	xyouts, 3., 1. + (max(heightp) + min(heightp)) / 3., 'Lat = '+string(latu,'(F5.1)')+'!Uo!N', CHARSIZE=1
	955	xyouts, 3., 0.5 + (max(heightp) + min(heightp)) / 3., 'LonE = '+string(lonu,'(F6.1)')+'!Uo!N', CHARSIZE=1
	956	xyouts, 3., (max(heightp) + min(heightp)) / 3., 'T!Dt!N = '+string(t_top_plus,'(I0)')+'/'+string(t_top_moins,'(I0)')+' K ', CHARSIZE=1
	957	xyouts, 3., -0.5 + (max(heightp) + min(heightp)) / 3., 'p!Dt!N = '+string(p_top_plus,'(I0)')+'/'+string(p_top_moins,'(I0)')+' Pa', CHARSIZE=1
	958	xyouts, 3., -1. + (max(heightp) + min(heightp)) / 3., 'z!Dt!N = '+string(z_top_plus,'(F4.1)')+'/'+string(z_top_moins,'(F4.1)')+' km', CHARSIZE=1
	959	xyouts, 3., -1.5 + (max(heightp) + min(heightp)) / 3., 'z!Ds!N = '+string(hgtu/1000.,'(F4.1)')+' km', CHARSIZE=1
	960	xyouts, 3., -2. + (max(heightp) + min(heightp)) / 3., 'D = '+string(pbl_depth,'(F4.1)')+' km', CHARSIZE=1
	961	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	962	no_staticstab:
	963	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	964
	965
	966
	967	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	968	if (saveps eq 'true') then begin
	969	device, /close
	970	endif
	971	stop
	972	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	973
	974
	975	user_lt=13.
	976	yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
	977
	978	user_h=0.1
	979	walt=where(abs(heightp-user_h) eq (min(abs(heightp-user_h))))
	980
	981	;mapfield=reform(t[,,walt(0),w(0)])
	982	;help, mapfield
	983	;contour, mapfield
	984
	985	section=reform(w[,0,,yeah(0)])
	986
	987	section=reform(w[,0,,160])
	988
	989	lev=[-12.,-8.,-4.,4.,8.,12.]
	990	lev=[-5.,-4.,-3.,-2.,-1.,1.,2.,3.,4.,5.]
	991	contour, $
	992	section, $
	993	(lon(,0)-lon(0,0))59., $
	994	heightp, $
	995	levels=lev , $
	996	C_LINESTYLE = (lev LT 0.0)
	997
	998	device, /close
	999
	1000	end
	1001

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

Download in other formats:

Original Format