Context Navigation

getturb.pro @ 1306

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

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

Download in other formats:

Original Format