source: trunk/MESOSCALE/PLOT/RESERVE/wind_analysis.pro @ 163

Last change on this file since 163 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: 2.3 KB
Line 
1pro wind_analysis
2
3;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
5retrieve='true'
6;retrieve='false'
7file='zefolder/wrfout_d01_2024-03-21_00:00:00'
8start_time=0.
9interv=1.
10landing_site=[353.87-360.,-1.88]
11;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
12;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13
14;
15; graphics definition
16;
17PS_Start, filename='wind.ps'
18!P.Charsize = 1.2
19
20;!p.charthick = 2.0
21;!p.thick = 3.0
22;!x.thick = 2.0
23;!y.thick = 2.0
24
25if (retrieve eq 'true') then begin
26
27;
28; load WRF coordinates
29;
30getcdf, file=file, charvar='XLAT', invar=lat
31getcdf, file=file, charvar='XLONG', invar=lon
32getcdf, file=file, charvar='HGT', invar=elevation
33;
34; find simulated lander grid point
35;
36findxy, $
37        x=lon, $
38        y=lat, $
39        point=landing_site, $
40        ind=indp, $
41        flag=2
42;
43; get field
44;
45getcdf, file=file, charvar='U', invar=u
46getcdf, file=file, charvar='V', invar=v
47;;field=sqrt(u^2+v^2) & u=0. & v=0.
48;;what_I_plot=reform(field(indp[0],indp[1],*,*))
49  ;;; BEWARE: arrays not the same size... because staggering... u^2+v^2 is wrong and IDL does not complain...
50field=sqrt( u(indp[0],indp[1],*,*)^2 + v(indp[0],indp[1],*,*)^2 ) & u=0. & v=0.
51what_I_plot=reform(field)
52
53;
54; local time
55;
56localtime = float(start_time) + float(interv)*findgen(n_elements(field(0,0,0,*))) ;+ lon(indp[0],indp[1])/15.
57;
58; get model height
59;
60getcdf, file=file, charvar='PHTOT', invar=ph
61height=reform(ph(indp[0],indp[1],*,*))
62height=total(height,2)/n_elements(height(0,*))
63height=height/1000./3.72
64heightp=height(0:n_elements(height(*))-2)-elevation(indp[0],indp[1])/1000.
65;
66; save
67;
68save, what_I_plot, localtime, heightp, filename='yeahyeah'
69
70endif else begin
71
72restore, filename='yeahyeah'   
73       
74endelse
75
76;yeah=where(localtime lt 0) & if (yeah(0) ne -1) then localtime[yeah]=24+localtime[yeah]
77;;; 24h issue
78;localtime=localtime MOD 24 & yeah=where(localtime eq 0) & if (yeah(0) ne -1) then localtime[yeah]=24
79;localtime=24.+localtime
80
81map_latlon, $
82        transpose(what_I_plot), $      ; 2D field
83        localtime, $                   ; 1D latitude
84        heightp, $                      ; 1D longitude
85        overcontour=0
86;        overcontour=transpose(what_I_plot)     ; another 2D field to overplot with contour lines (=0: no)
87
88PS_End, /PNG
89
90end
Note: See TracBrowser for help on using the repository browser.