Context Navigation

getturb.pro @ 86

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

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

Download in other formats:

Original Format