| 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 |
|---|
