source: trunk/MESOSCALE/PLOT/SPEC/LES/getturb.pro @ 205

Last change on this file since 205 was 113, checked in by aslmd, 14 years ago

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