Context Navigation

getturb.pro @ 105

Last change on this file since 105 was 96, checked in by aslmd, 14 years ago
LMD_MM_MARS et LMD_LES_MARS: ajouts mineurs: commentaires et modifications routines de plot
File size: 46.0 KB

Line
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
16	;;
17	;; constantes
18	;;
19	;p0=610. & t0=220. & r_cp=1/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
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
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
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
173	@include_ustar.pro
174
175	;@include_hfx.pro
176
177	;@include_w.pro
178
179	;@include_tprime.pro
180
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
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
377	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
378	zefield = localtime
379	zey = pbl_height
380	set_name = 'HEIGHT.ps'
381	set_title = 'Boundary layer height (km)'
382	set_titlex = 'Local Time (h)'
383	set_titley = 'Boundary layer height (km)'
384	set_subtitle = '' ;'Criterion is : static stability > 1.5 K.m!U-1!N'
385	set_xrange = [8.,17.]
386	set_yrange = [0.,7.]
387	set_tickx = 1.
388	set_ticky = 1.
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	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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
440	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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>]
446	set_titlex = 'Local Time (h)'
447	set_titley = 'Altitude above surface (km)'
448	set_subtitle = '' ;'Mean over the simulation domain'
449	set_xrange = [8.,17.]
450	set_yrange = [0.,7.] ;; [0.,200.]
451	set_tickx = 1.
452	set_ticky = 1. ;; 50.
453	minval = 0.
454	maxval = 15.
455	nlev = maxval-minval
456	pal = 22
457	rrr = 'no'
458	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
459	if (saveps eq 'true') then PS_Start, FILENAME=set_name
460	!P.Charsize = 1.2
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
471	if (saveps eq 'true') then PS_End, /PNG
472	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
473
474	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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)
482	zex = localtime
483	zey = h
484	set_name = 'zTKE.ps'
485	set_title = "Vertical Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<w'!U2!N>]
486	set_titlex = 'Local Time (h)'
487	set_titley = 'Altitude above surface (km)'
488	set_subtitle = '' ;'Mean over the simulation domain'
489	set_xrange = [8.,17.]
490	set_yrange = [0.,7.]
491	set_tickx = 1.
492	set_ticky = 1.
493	minval = 0.
494	maxval = 8.
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	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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
551	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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)'
558	set_titley = 'Altitude above surface (km)'
559	set_subtitle = '' ;'Mean over the simulation domain'
560	set_xrange = [8.,17.]
561	set_yrange = [0.,7.] ;; [0.,200.]
562	set_tickx = 1.
563	set_ticky = 1. ;; 50.
564	minval = -2.
565	maxval = 2.
566	nlev = floor(maxval-minval)*10
567	pal = 33
568	rrr = 'no'
569	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
570	;minval = -0.8
571	;maxval = 1.2
572	;pal = 0
573	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
574	if (saveps eq 'true') then PS_Start, FILENAME=set_name
575	!P.Charsize = 1.2
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
592	if (saveps eq 'true') then PS_End, /PNG
593	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
594
595	;;restore, filename='tpot_profB'
596	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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.
608	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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]
613	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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)
721	zex = localtime
722	zey = h
723	set_name = 'hTKE.ps'
724	set_title = "Horizontal Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<v'!U2!N>]
725	set_titlex = 'Local Time (h)'
726	set_titley = 'Altitude above surface (km)'
727	set_subtitle = '' ;'Mean over the simulation domain'
728	set_xrange = [8.,17.]
729	set_yrange = [0.,1.] ;[0.,7.] ;; [0.,200.]
730	set_tickx = 1.
731	set_ticky = 0.1 ;1. ;50.
732	minval = 0.
733	maxval = 12. ;10.
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	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
752
753	if (n_elements(wmax) eq 0) then stop
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 = ''
764	set_xrange = [8.,18.]
765	set_yrange = [0.,7.]
766	set_tickx = 1.
767	set_ticky = 1.
768	minval = 0.
769	maxval = 16. ;15.;14.;13.
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.)]
778	;;; 1. background
779	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
780	;; 2. color field
781	loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
782	contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
783	;; 3. contour field
784	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)
785	;; 4. choose output
786	if (saveps eq 'true') then PS_End, /PNG
787	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
788
789	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
790	zefield = abs(transpose(wmin))
791	zex = localtime
792	zey = h
793	set_name = 'Wmin.ps'
794	set_title = "Maximum downdraft speed (m.s!U-1!N)"
795	set_titlex = 'Local Time (h)'
796	set_titley = 'Altitude above surface (km)'
797	set_subtitle = ''
798	set_xrange = [8.,18.]
799	set_yrange = [0.,7.]
800	set_tickx = 1.
801	set_ticky = 1.
802	minval = 0.
803	maxval = 9.;12.
804	nlev = maxval-minval
805	pal = 22
806	rrr = 'no' ;'yes'
807	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
808	if (saveps eq 'true') then PS_Start, FILENAME=set_name
809	!P.Charsize = 1.2
810	;; 0. levels
811	lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
812	;;; 1. background
813	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
814	;; 2. color field
815	loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
816	contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
817	;; 3. contour field
818	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)
819	;; 4. choose output
820	if (saveps eq 'true') then PS_End, /PNG
821	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
822
823
824	restore, filename='addturb.dat'
825	velmax = SMOOTH(TEMPORARY(velmax), [0,smoothampl], /EDGE_TRUNCATE)
826	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
827	zefield = transpose(velmax)
828	zex = localtime
829	zey = h
830	set_name = 'velmax.ps'
831	set_title = "Maximum horizontal wind speed (m.s!U-1!N)"
832	set_titlex = 'Local Time (h)'
833	set_titley = 'Altitude above surface (km)'
834	set_subtitle = '' ;'Mean over the simulation domain'
835	set_xrange = [8.,17.]
836	set_yrange = [0.,4.]
837	set_tickx = 1.
838	set_ticky = 1.
839	minval = 0.
840	maxval = 30. ;12.
841	nlev = maxval-minval
842	pal = 22
843	rrr = 'no'
844	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
845	if (saveps eq 'true') then PS_Start, FILENAME=set_name
846	!P.Charsize = 1.2
847	;; 0. levels
848	lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)]
849	;;; 1. background
850	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
851	;;; 2. color field
852	loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b)
853	contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill
854	;; 3. contour field
855	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)
856	;; 4. choose output
857	if (saveps eq 'true') then PS_End, /PNG
858	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
859
860
861	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
862	;zefield = localtime
863	;zey = depressions
864	;set_name = 'DEPR.ps'
865	;set_title = "Percentage of depressions in the domain"
866	;set_titlex = 'Local Time (h)'
867	;set_titley = 'Percentage'
868	;set_subtitle = ''
869	;set_xrange = [8.,18.]
870	;set_yrange = [0.,3.]
871	;set_tickx = 1.
872	;set_ticky = 1.
873	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
874	;if (saveps eq 'true') then PS_Start, FILENAME=set_name
875	;!P.Charsize = 1.2
876	;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
877	;oplot, zefield, zey, psym=5
878	; oplot, localtime, abs(psmin), linestyle=1
879	;if (saveps eq 'true') then PS_End, /PNG
880	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
881
882	stop
883
884
885
886
887	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
888	goto, no_staticstab
889	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
890	if (saveps eq 'true') then begin
891	device, /close
892	set_plot, 'ps' & device, filename='plot/staticstab.ps'
893	endif
894
895	user_lt=17. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
896
897	;;; recompute height @ given local time
898	caca=heightp+hgtu/1000.
899	height=reform(ph(,0,,*))
900	height=total(height,1)/n_elements(height(*,0,0))
901	height=reform(height)
902	height=reform(height(*,yeah(0)))
903	height=height/1000./3.72
904	heightp=height(0:n_elements(height(*))-2)
905	print, 'new minus old', heightp - caca
906	;;; recompute height @ given local time
907
908	press=reform(p(,0,,*))
909	press=total(press,1)/n_elements(press(*,0,0))
910	press=reform(press)
911	press=reform(press(*,yeah(0)))
912	press=press(0:n_elements(press(*))-2)
913
914	staticstab = DERIV( heightp , reform(what_I_plot5(,yeah(0))) ) + 1000.3.72/844.6 ;; 4.9 dans Hinson
915	staticstab = staticstab(0:n_elements(staticstab)-5)
916	heightp = heightp(0:n_elements(heightp)-5)
917
918	plot, $
919	staticstab, $
920	heightp,$
921	xtitle='Static stability (K/km)',$
922	xrange=[-1.,5.], $
923	xtickinterval=0.5, $
924	ytitle='Altitude above surface (km)', $
925	yrange=[min(heightp),max(heightp)], $
926	ytickinterval=1., $
927	title="LMD LES Static stability @ LT 17h (K/km)", $
928	subtitle='zonal average at lat. '+latwrite
929	oplot, $
930	staticstab, $
931	heightp,$
932	psym=5
933	oplot, $
934	findgen(n_elements(heightp))*0. + 1.,$
935	heightp,$
936	linestyle=2
937	oplot, $
938	findgen(n_elements(heightp))*0. + 2.,$
939	heightp,$
940	linestyle=2
941
942	;;;;;;;;;;
943	w = where((staticstab gt 1.5) and ((heightp- hgtu/1000.) gt 1.))
944	t_top_plus = what_I_plot5(w(0),yeah(0))
945	z_top_plus = heightp(w(0))
946	p_top_plus = press(w(0))
947	t_top_moins = what_I_plot5(w(0)-1,yeah(0))
948	z_top_moins = heightp(w(0)-1)
949	p_top_moins = press(w(0)-1)
950	pbl_depth = (z_top_plusalog(p_top_plus) + z_top_moinsalog(p_top_moins))/(alog(p_top_plus) + alog(p_top_moins)) - hgtu/1000.
951	xyouts, 3., 1.5 + (max(heightp) + min(heightp)) / 3., 'Ls = '+string(lsu,'(F5.1)')+'!Uo!N', CHARSIZE=1
952	xyouts, 3., 1. + (max(heightp) + min(heightp)) / 3., 'Lat = '+string(latu,'(F5.1)')+'!Uo!N', CHARSIZE=1
953	xyouts, 3., 0.5 + (max(heightp) + min(heightp)) / 3., 'LonE = '+string(lonu,'(F6.1)')+'!Uo!N', CHARSIZE=1
954	xyouts, 3., (max(heightp) + min(heightp)) / 3., 'T!Dt!N = '+string(t_top_plus,'(I0)')+'/'+string(t_top_moins,'(I0)')+' K ', CHARSIZE=1
955	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
956	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
957	xyouts, 3., -1.5 + (max(heightp) + min(heightp)) / 3., 'z!Ds!N = '+string(hgtu/1000.,'(F4.1)')+' km', CHARSIZE=1
958	xyouts, 3., -2. + (max(heightp) + min(heightp)) / 3., 'D = '+string(pbl_depth,'(F4.1)')+' km', CHARSIZE=1
959	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
960	no_staticstab:
961	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
962
963
964
965	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
966	if (saveps eq 'true') then begin
967	device, /close
968	endif
969	stop
970	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
971
972
973	user_lt=13.
974	yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
975
976	user_h=0.1
977	walt=where(abs(heightp-user_h) eq (min(abs(heightp-user_h))))
978
979	;mapfield=reform(t[,,walt(0),w(0)])
980	;help, mapfield
981	;contour, mapfield
982
983	section=reform(w[,0,,yeah(0)])
984
985	section=reform(w[,0,,160])
986
987	lev=[-12.,-8.,-4.,4.,8.,12.]
988	lev=[-5.,-4.,-3.,-2.,-1.,1.,2.,3.,4.,5.]
989	contour, $
990	section, $
991	(lon(,0)-lon(0,0))59., $
992	heightp, $
993	levels=lev , $
994	C_LINESTYLE = (lev LT 0.0)
995
996	device, /close
997
998	end
999

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

Download in other formats:

Original Format