1 | pro wind_analysis |
---|
2 | |
---|
3 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
---|
4 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
---|
5 | retrieve='true' |
---|
6 | ;retrieve='false' |
---|
7 | file='zefolder/wrfout_d01_2024-03-21_00:00:00' |
---|
8 | start_time=0. |
---|
9 | interv=1. |
---|
10 | landing_site=[353.87-360.,-1.88] |
---|
11 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
---|
12 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
---|
13 | |
---|
14 | ; |
---|
15 | ; graphics definition |
---|
16 | ; |
---|
17 | PS_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 | |
---|
25 | if (retrieve eq 'true') then begin |
---|
26 | |
---|
27 | ; |
---|
28 | ; load WRF coordinates |
---|
29 | ; |
---|
30 | getcdf, file=file, charvar='XLAT', invar=lat |
---|
31 | getcdf, file=file, charvar='XLONG', invar=lon |
---|
32 | getcdf, file=file, charvar='HGT', invar=elevation |
---|
33 | ; |
---|
34 | ; find simulated lander grid point |
---|
35 | ; |
---|
36 | findxy, $ |
---|
37 | x=lon, $ |
---|
38 | y=lat, $ |
---|
39 | point=landing_site, $ |
---|
40 | ind=indp, $ |
---|
41 | flag=2 |
---|
42 | ; |
---|
43 | ; get field |
---|
44 | ; |
---|
45 | getcdf, file=file, charvar='U', invar=u |
---|
46 | getcdf, 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... |
---|
50 | field=sqrt( u(indp[0],indp[1],*,*)^2 + v(indp[0],indp[1],*,*)^2 ) & u=0. & v=0. |
---|
51 | what_I_plot=reform(field) |
---|
52 | |
---|
53 | ; |
---|
54 | ; local time |
---|
55 | ; |
---|
56 | localtime = float(start_time) + float(interv)*findgen(n_elements(field(0,0,0,*))) ;+ lon(indp[0],indp[1])/15. |
---|
57 | ; |
---|
58 | ; get model height |
---|
59 | ; |
---|
60 | getcdf, file=file, charvar='PHTOT', invar=ph |
---|
61 | height=reform(ph(indp[0],indp[1],*,*)) |
---|
62 | height=total(height,2)/n_elements(height(0,*)) |
---|
63 | height=height/1000./3.72 |
---|
64 | heightp=height(0:n_elements(height(*))-2)-elevation(indp[0],indp[1])/1000. |
---|
65 | ; |
---|
66 | ; save |
---|
67 | ; |
---|
68 | save, what_I_plot, localtime, heightp, filename='yeahyeah' |
---|
69 | |
---|
70 | endif else begin |
---|
71 | |
---|
72 | restore, filename='yeahyeah' |
---|
73 | |
---|
74 | endelse |
---|
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 | |
---|
81 | map_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 | |
---|
88 | PS_End, /PNG |
---|
89 | |
---|
90 | end |
---|