source: trunk/mesoscale/PLOT/SPEC/LES/turbulence.pro @ 85

Last change on this file since 85 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: 15.9 KB
Line 
1pro turbulence
2
3;;----------------------------------------------------
4;;----------------------------------------------------
5;;----------------------------------------------------
6
7
8retrieve='true'
9;retrieve='false'
10
11
12
13;; CALCUL DE LA MOYENNE
14;; 18 pas assez 74 trop
15meansmooth=37
16;meansmooth=18
17;meansmooth=74
18
19;; CALCUL DE LA MOYENNE DES FLUX
20;smoothampl=74
21;smoothampl=2
22;smoothampl=10
23;smoothampl=18
24smoothampl=meansmooth
25
26saveps='false'
27;saveps='true'
28
29yu=20
30
31;;----------------------------------------------------
32;;----------------------------------------------------
33;;----------------------------------------------------
34
35
36p0=610. & t0=220. & r_cp=1/4.4 & grav=3.72 & R=192.
37
38;
39; graphics definition
40;
41if (saveps eq 'false') then begin
42   !p.multi=[0,3,2]
43   !P.CHARSIZE=2.
44endif else begin
45   !p.charthick = 2.0
46   !p.thick = 3.0
47   !x.thick = 2.0
48   !y.thick = 2.0
49endelse
50
51;
52; input files
53;
54OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22
55OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
56print, lonu, latu, lsu, lctu, hgtu, tsurfu
57
58
59;
60; retrieve
61;
62if (retrieve eq 'true') then begin
63
64
65;
66; get fields
67;
68getcdf, file='u'+string(yu,'(I0)')+'.nc', charvar='U', invar=u
69getcdf, file='v'+string(yu,'(I0)')+'.nc', charvar='V', invar=v
70getcdf, file='t'+string(yu,'(I0)')+'.nc', charvar='T', invar=t
71getcdf, file='w'+string(yu,'(I0)')+'.nc', charvar='W', invar=w
72getcdf, file='p'+string(yu,'(I0)')+'.nc', charvar='PTOT', invar=p
73        ;
74        ; no vertical staggering for w
75        ;
76        wplus=shift(w,0,0,-1,0)
77        w=(w+wplus)/2.
78        w=w(*,*,0:n_elements(w(0,0,*,0))-2,*)
79
80
81ndim=n_elements(t(*,0,0,0))
82lat=fltarr(ndim,ndim)
83lon=fltarr(ndim,ndim)
84  for i=0,ndim-1 do begin
85  for j=0,ndim-1 do begin
86    lat(i,j)=latu+float(j)*100./59000.
87    lon(i,j)=lonu+float(i)*100./59000.
88  endfor
89  endfor
90what_I_plot=0. & what_I_plot2=0. & what_I_plot3=0. & what_I_plot4=0. & what_I_plot5=0. & what_I_plot6=0.
91maxwhat_I_plot=0. & maxwhat_I_plot2=0. & maxwhat_I_plot3=0.
92minwhat_I_plot=0.
93
94
95;
96; loop on west-east coordinate
97;
98nloop=n_elements(t(*,0,0,0))
99for i=0,nloop-1 do begin
100
101        indp=[i,0]      ;;indp=[75,0]
102
103        ;
104        ; 2D field for plot
105        ;
106        t2=reform(t(indp[0],indp[1],*,*))
107t2_save=t0+reform(t(indp[0],indp[1],*,*))
108treal2_save=t2_save*(reform(p(indp[0],indp[1],*,*))/p0)^r_cp
109        u2=reform(u(indp[0],indp[1],*,*))
110        v2=reform(v(indp[0],indp[1],*,*))
111        w2=reform(w(indp[0],indp[1],*,*))
112        nz=n_elements(t2(*,0))
113        nt=n_elements(t2(0,*))
114
115nt=nt-2
116t2=t2[*,0:nt-1]
117u2=u2[*,0:nt-1]
118v2=v2[*,0:nt-1]
119w2=w2[*,0:nt-1]
120
121                ;;
122                ;; local time
123                ;;
124                ;localtime=21.+100.*findgen(nt)/3700.+lon(indp[0],indp[1])/15.-24.
125                localtime=lctu+100.*findgen(nt)/3700.
126                ;;print, localtime
127               
128        ;
129        ; 1. perturbations
130        ;
131        for k=0,nz-1 do begin
132                t2(k,*)=t2(k,*)-SMOOTH(reform(t2(k,*)),meansmooth,/EDGE_TRUNCATE)
133                u2(k,*)=u2(k,*)-SMOOTH(reform(u2(k,*)),meansmooth,/EDGE_TRUNCATE)
134                v2(k,*)=v2(k,*)-SMOOTH(reform(v2(k,*)),meansmooth,/EDGE_TRUNCATE)
135                w2(k,*)=w2(k,*);-SMOOTH(reform(w2(k,*)),meansmooth,/EDGE_TRUNCATE)
136                ;treal2_save(k,*)=treal2_save(k,*)-SMOOTH(reform(treal2_save(k,*)),meansmooth,/EDGE_TRUNCATE)
137        endfor
138        ;
139        ; 2. nonlinear products
140        ;
141        vertical_eddy_heat_flux=w2*t2
142        uprime2=u2^2
143        vprime2=v2^2
144        wprime2=w2^2
145        ;
146        ; 3. Reynolds averaging
147        ;
148        for k=0,nz-1 do begin
149                vertical_eddy_heat_flux(k,*)=SMOOTH(reform(vertical_eddy_heat_flux(k,*)),smoothampl,/EDGE_TRUNCATE)
150                uprime2(k,*)=SMOOTH(reform(uprime2(k,*)),smoothampl,/EDGE_TRUNCATE)
151                vprime2(k,*)=SMOOTH(reform(vprime2(k,*)),smoothampl,/EDGE_TRUNCATE)
152                wprime2(k,*)=SMOOTH(reform(wprime2(k,*)),smoothampl,/EDGE_TRUNCATE)
153        endfor
154        tke=0.5*(uprime2+vprime2+wprime2)
155
156
157what_I_plot=what_I_plot+vertical_eddy_heat_flux/nloop
158what_I_plot2=what_I_plot2+tke/nloop
159what_I_plot3=what_I_plot3+wprime2/nloop
160what_I_plot4=what_I_plot4+t2_save/nloop
161what_I_plot5=what_I_plot5+treal2_save/nloop
162if (i eq nloop-1) then  what_I_plot6=t2
163
164;;help, t2
165;;help, treal2_save
166;treal2_save=treal2_save(0:n_elements(t2(*,0))-1,0:n_elements(t2(0,*))-1)
167;if (i eq nloop-1) then  what_I_plot6=treal2_save
168
169yeah=max(vertical_eddy_heat_flux)
170yeah2=max(tke)
171yeah3=max(wprime2)
172yeah4=min(vertical_eddy_heat_flux)
173
174if (yeah gt maxwhat_I_plot) then maxwhat_I_plot=yeah
175if (yeah2 gt maxwhat_I_plot2) then maxwhat_I_plot2=yeah2
176if (yeah3 gt maxwhat_I_plot3) then maxwhat_I_plot3=yeah3
177if (yeah4 lt minwhat_I_plot) then minwhat_I_plot=yeah4
178
179endfor
180
181
182
183
184
185;
186; get model height
187;
188getcdf, file='ph'+string(yu,'(I0)')+'.nc', charvar='PHTOT', invar=ph
189height=reform(ph(*,0,*,*))   
190height=total(height,1)/n_elements(height(*,0,0))
191height=reform(height)
192height=total(height(*,1:nt-1),2)/(nt-1)
193height=height/1000./3.72
194heightp=height(0:n_elements(height(*))-2)
195        ;
196        ; altitude above ground
197        ;       
198        heightp=heightp-hgtu/1000.
199
200
201;
202; save
203;
204save, what_I_plot, $
205        what_I_plot2, $
206        what_I_plot3, $
207        what_I_plot4, $
208        what_I_plot5, $
209        what_I_plot6, $
210        localtime, $
211        heightp, $
212        lat, $
213        lon, $
214        hgtu, $
215        filename='turbulence.dat'
216
217endif else begin
218        restore, filename='turbulence.dat'     
219endelse
220
221help, what_I_plot
222help, localtime
223help, heightp
224
225
226latwrite=string(lat(0,yu),'(F5.1)')
227
228
229
230
231;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
232if (saveps eq 'true') then begin
233  set_plot, 'ps' & device, filename='plot/TKE.ps'
234endif
235
236lev=findgen(15)+1.
237contour, $
238        transpose(what_I_plot2), $
239        localtime, $
240        heightp, $
241        xtitle='Local Time (h)', $
242        xrange=[8.,18.], $
243        xtickinterval=1., $
244        ytitle='Altitude above surface (km)', $
245        yrange=[0.,10.], $
246        ytickinterval=1., $
247        title="Turbulent Kinetic Energy 0.5[<u'!U2!N>+<v'!U2!N>+<w'!U2!N>] (m!U2!N.s!U-2!N)", $
248        levels=lev, $
249        c_labels=findgen(n_elements(lev))*0.+1., $
250        subtitle='zonal average at lat. '+latwrite+' (max value is '+string(maxwhat_I_plot2,'(F4.1)')+' m!U2!N.s!U-2!N)'
251;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
252
253
254;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
255if (saveps eq 'true') then begin
256  device, /close
257  set_plot, 'ps' & device, filename='plot/zTKE.ps'
258endif
259
260lev=findgen(15)+1.
261contour, $
262        0.5*transpose(what_I_plot3), $
263        localtime, $
264        heightp, $
265        xtitle='Local Time (h)', $
266        xrange=[8.,18.], $
267        xtickinterval=1., $
268        ytitle='Altitude above surface (km)', $
269        yrange=[0.,10.], $
270        ytickinterval=1., $
271        title="Vertical TKE 0.5<w'!U2!N> (m!U2!N.s!U-2!N)", $
272        levels=lev, $
273        c_labels=findgen(n_elements(lev))*0.+1., $
274        subtitle='zonal average at lat. '+latwrite+' (max value is '+string(0.5*maxwhat_I_plot3,'(F4.1)')+' m!U2!N.s!U-2!N)'
275;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
276
277
278
279
280if (saveps eq 'true') then begin
281  device, /close
282  set_plot, 'ps' & device, filename='plot/HF.ps'
283endif
284
285;lev=[-0.15,-0.05,-0.02,0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95]
286lev = -1. + findgen(20)/10.
287contour, $
288        transpose(what_I_plot), $
289        localtime, $
290        heightp, $
291        xtitle='Local Time (h)', $
292        xrange=[8.,18.], $
293        xtickinterval=1., $
294        ytitle='Altitude above surface (km)', $
295        yrange=[0.,10.], $
296        ytickinterval=1., $
297        title="Vertical Eddy Heat Flux <w'!7h!3'> (K.m.s!U-1!N)", $
298        levels=lev, $   
299        C_LINESTYLE = (lev LT 0.0), $
300        c_labels=findgen(n_elements(lev))*0.+1., $
301        subtitle='zonal average at lat. '+latwrite+' (min/max values are '+string(minwhat_I_plot,'(F5.1)')+'/'+string(maxwhat_I_plot,'(F3.1)')+' K.m.s!U-1!N)'
302
303
304
305
306
307;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
308goto, no_pert_temp
309;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
310if (saveps eq 'true') then begin
311  device, /close
312  set_plot, 'ps' & device, filename='plot/pert_temp.ps'
313endif
314
315lev=[-4.,-2.,-1.,-0.5,0.5,1.,2.,4.]
316contour, $
317        transpose(what_I_plot6), $
318        localtime, $
319        heightp, $
320        xtitle='Local Time (h)', $
321        xrange=[13.5,14.], $
322        xtickinterval=0.1, $
323        ytitle='Altitude above surface (km)', $
324        yrange=[0.,0.5], $
325        ytickinterval=0.1, $
326        title="Potential Temperature Perturbation !7h!3' (K)", $
327        levels=lev, $
328;nlevels=20, $
329        C_LINESTYLE = (lev LT 0.0), $
330        c_labels=findgen(n_elements(lev))*0.+1., $
331;/cell_fill, $
332        subtitle='latitude '+latwrite+' longitude '+string(lon(nloop-1,0),'(F6.1)')
333;xyouts, 14.007, heightp(0), '1'
334;xyouts, 14.016, heightp(1), '2'
335;xyouts, 14.007, heightp(2), '3'
336;xyouts, 14.007, heightp(3), '4'
337;xyouts, 14.007, heightp(4), '5'
338;xyouts, 14.007, heightp(5), '6'
339;xyouts, 14.007, heightp(6), '7'
340;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
341no_pert_temp:
342;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
343
344
345
346
347
348if (saveps eq 'true') then begin
349  device, /close
350  set_plot, 'ps' & device, filename='plot/temp_pot.ps'
351endif
352
353user_lt=10. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
354plot, $
355        what_I_plot4(*,yeah(0)), $
356        heightp,$
357        xtitle='Potential Temperature (K)',$
358;       xrange=[215,240], $
359        xtickinterval=5., $
360        ytitle='Altitude above surface (km)', $
361        yrange=[0.,10.], $
362        ytickinterval=1., $
363        title="Potential Temperature !7h!3 profile (K)", $
364        subtitle='zonal average at lat. '+latwrite
365user_lt=12. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
366        oplot, what_I_plot4(*,yeah(0)), heightp, linestyle=1
367user_lt=14. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
368        oplot, what_I_plot4(*,yeah(0)), heightp
369user_lt=16. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
370        oplot, what_I_plot4(*,yeah(0)), heightp, linestyle=1
371user_lt=18. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
372        oplot, what_I_plot4(*,yeah(0)), heightp
373
374;xyouts, 217.3, 0.5, '10:00'
375;xyouts, 225.8, 1.5, '12:00'
376;xyouts, 231.8, 3.0, '14:00'
377;xyouts, 234.8, 0.7, '16:00'
378;xyouts, 232.3, 2.0, '18:00'
379
380
381if (saveps eq 'true') then begin
382  device, /close
383  set_plot, 'ps' & device, filename='plot/temp.ps'
384endif
385
386user_lt=14.+55./60. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
387plot, $
388        what_I_plot5(*,yeah(0)), $
389        heightp,$
390        xtitle='Temperature (K)',$
391        xrange=[210,250], $
392        xtickinterval=5., $
393        ytitle='Altitude above surface (km)', $
394        yrange=[0.,2.], $
395        ytickinterval=0.2, $
396        title="LMD LES Near-surface T profile (K)", $
397        subtitle='zonal average at lat. '+latwrite, $
398        linestyle=2     
399user_lt=10.+24./60. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
400        oplot, what_I_plot5(*,yeah(0)), heightp
401user_lt=10.+59./60. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
402        oplot, what_I_plot5(*,yeah(0)), heightp, linestyle=2
403user_lt=11.+59./60. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
404        oplot, what_I_plot5(*,yeah(0)), heightp
405user_lt=15.+22./60. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
406        oplot, what_I_plot5(*,yeah(0)), heightp
407user_lt=17.+13./60. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
408        oplot, what_I_plot5(*,yeah(0)), heightp, linestyle=2
409
410;xyouts, 221,   0.2,  '10:24'
411;xyouts, 226.5,   0.4,  '10:59'
412;xyouts, 231.8, 0.2,  '11:59'
413;xyouts, 233,   0.8,  '14:55'
414;xyouts, 241,   0.1,  '15:22'
415;xyouts, 236.5, 0.05, '17:13'
416 
417
418;;
419;xyouts, 250.05, heightp(0), '1'
420;xyouts, 250.45, heightp(1), '2'
421;xyouts, 251.00, heightp(2), '3'
422;xyouts, 250.45, heightp(3), '4'
423;xyouts, 250.45, heightp(4), '5'
424;xyouts, 250.45, heightp(5), '6'
425;xyouts, 250.45, heightp(6), '7'
426;xyouts, 250.45, heightp(7), '8'
427;xyouts, 250.45, heightp(8), '9'
428;xyouts, 250.25, heightp(9), '10'
429;xyouts, 250.25, heightp(10), '11'
430;xyouts, 250.25, heightp(11), '12'
431;xyouts, 250.25, heightp(12), '13'
432;xyouts, 250.25, heightp(13), '14'
433;xyouts, 250.25, heightp(14), '15'
434;xyouts, 250.25, heightp(15), '16'
435
436
437
438;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
439;goto, no_staticstab
440;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
441if (saveps eq 'true') then begin
442  device, /close
443  set_plot, 'ps' & device, filename='plot/staticstab.ps'
444endif
445
446user_lt=17. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
447
448   ;;; recompute height @ given local time
449   caca=heightp+hgtu/1000.
450   height=reform(ph(*,0,*,*))
451   height=total(height,1)/n_elements(height(*,0,0))
452   height=reform(height)
453   height=reform(height(*,yeah(0)))
454   height=height/1000./3.72
455   heightp=height(0:n_elements(height(*))-2)
456   print, 'new minus old', heightp - caca
457   ;;; recompute height @ given local time
458
459        press=reform(p(*,0,*,*))
460        press=total(press,1)/n_elements(press(*,0,0))
461        press=reform(press)
462        press=reform(press(*,yeah(0)))
463        press=press(0:n_elements(press(*))-2)
464
465staticstab = DERIV( heightp , reform(what_I_plot5(*,yeah(0))) ) + 1000.*3.72/844.6  ;; 4.9 dans Hinson
466   staticstab = staticstab(0:n_elements(staticstab)-5)
467   heightp = heightp(0:n_elements(heightp)-5)
468
469plot, $
470        staticstab, $
471        heightp,$
472        xtitle='Static stability (K/km)',$
473        xrange=[-1.,5.], $
474        xtickinterval=0.5, $
475        ytitle='Altitude above surface (km)', $
476        yrange=[min(heightp),max(heightp)], $
477        ytickinterval=1., $
478        title="LMD LES Static stability @ LT 17h (K/km)", $
479        subtitle='zonal average at lat. '+latwrite
480oplot, $
481        staticstab, $
482        heightp,$
483        psym=5
484oplot, $
485       findgen(n_elements(heightp))*0. + 1.,$
486       heightp,$
487       linestyle=2
488oplot, $
489       findgen(n_elements(heightp))*0. + 2.,$
490       heightp,$
491       linestyle=2
492
493;;;;;;;;;;
494w = where((staticstab gt 1.5) and ((heightp- hgtu/1000.) gt 1.))
495t_top_plus = what_I_plot5(w(0),yeah(0))
496z_top_plus = heightp(w(0))
497p_top_plus = press(w(0))
498t_top_moins = what_I_plot5(w(0)-1,yeah(0))
499z_top_moins = heightp(w(0)-1)
500p_top_moins = press(w(0)-1)
501pbl_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.
502xyouts, 3., 1.5 + (max(heightp) + min(heightp)) / 3., 'Ls = '+string(lsu,'(F5.1)')+'!Uo!N', CHARSIZE=1
503xyouts, 3., 1. + (max(heightp) + min(heightp)) / 3., 'Lat = '+string(latu,'(F5.1)')+'!Uo!N', CHARSIZE=1
504xyouts, 3., 0.5 + (max(heightp) + min(heightp)) / 3., 'LonE = '+string(lonu,'(F6.1)')+'!Uo!N', CHARSIZE=1
505xyouts, 3., (max(heightp) + min(heightp)) / 3., 'T!Dt!N = '+string(t_top_plus,'(I0)')+'/'+string(t_top_moins,'(I0)')+' K ', CHARSIZE=1
506xyouts, 3., -0.5 + (max(heightp) + min(heightp)) / 3., 'p!Dt!N = '+string(p_top_plus,'(I0)')+'/'+string(p_top_moins,'(I0)')+' Pa', CHARSIZE=1
507xyouts, 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
508xyouts, 3., -1.5 + (max(heightp) + min(heightp)) / 3., 'z!Ds!N = '+string(hgtu/1000.,'(F4.1)')+' km', CHARSIZE=1
509xyouts, 3., -2. + (max(heightp) + min(heightp)) / 3., 'D = '+string(pbl_depth,'(F4.1)')+' km', CHARSIZE=1
510;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
511no_staticstab:
512;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
513
514
515
516;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
517if (saveps eq 'true') then begin
518  device, /close
519endif
520stop
521;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
522
523
524user_lt=13.
525yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
526
527user_h=0.1
528walt=where(abs(heightp-user_h) eq (min(abs(heightp-user_h))))
529
530;mapfield=reform(t[*,*,walt(0),w(0)])
531;help, mapfield
532;contour, mapfield
533
534section=reform(w[*,0,*,yeah(0)])
535
536section=reform(w[*,0,*,160])
537
538lev=[-12.,-8.,-4.,4.,8.,12.]
539lev=[-5.,-4.,-3.,-2.,-1.,1.,2.,3.,4.,5.]
540contour, $
541        section, $
542        (lon(*,0)-lon(0,0))*59., $
543        heightp, $
544        levels=lev , $
545        C_LINESTYLE = (lev LT 0.0)
546
547device, /close
548
549end
Note: See TracBrowser for help on using the repository browser.