source: trunk/mesoscale/PLOT/SPEC/GW/gravitwaveprof.pro @ 114

Last change on this file since 114 was 114, checked in by aslmd, 15 years ago

LMD_MM_MARS: update graphic tools for GW, water cycle + generic graphic stuff (ps_start, map_latlon) + update notes
  • Property svn:executable set to *
File size: 9.0 KB
Line 
1pro gravitwaveprof
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3
4;
5; faire avant api avec kappa 3.9
6;
7
8;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
9        path='../TESTGW/'
10        path='/tmp7/aslmd/'
11path='/d5/aslmd/GRAVITWAVE/'
12        experiment='mcd_lt16_rcpcst'
13        ;experiment='mcd_lt16'
14        ;experiment='short_mcd_lt16_rcpcst'
15        ;experiment='short_mcd_lt16'
16        ;experiment='short_mcd_lt16_calcprho'
17        ;experiment='mcd_lt16_calcprho'
18        ;experiment='short_mcd_lt16_calcprho_tot'
19        ;experiment='colder_gcm_add_short'
20        ;experiment='colder_gcm_add'
21        experiment=''
22        ;experiment='wind15'
23        ;experiment='wind20c'
24        ;experiment='colder_gcm_add_3.9_fine_rad_mountain'
25        experiment='test'
26experiment='GW_MARS_highwind_3D.157077'
27experiment='GW_MARS_highwind_3D_lower_opacity.161436'
28experiment='GW_MARS_highwind_3D_loweropacity_lotspoints.161468'
29experiment='GW_MARS_highwind_3D_loweropacity_morepoints.161542'
30experiment='GW_MARS_highwind_3D_loweropacity_morepoints_LT15.161544'
31        file=path+experiment+'_wrfout_d01_9999-09-09_09:00:00_z'
32file=path+experiment+'/wrfout_d01_9999-09-09_09:00:00_z'
33;file=path+experiment+'/wrfout_d01_9999-09-09_10:00:00_z'
34;file=path+experiment+'/wrfout_d01_9999-09-09_11:00:00_z'
35charvar='W' & charvarc='W'
36charvar='tk' & charvarc='W'
37charvar='tk' & charvarc='W' & cond=1
38charvar='tpot' & charvarc='tpot'
39;charvar='PTOT' & charvarc='PTOT'
40charvar='tk' & charvarc='tk'
41;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
42nts=1  & nte=15   & nxs=0   & nxe=59
43nts=3  & nte=nts  & nxs=999 & nxe=999 ;; moyenne ;; attention a la montagne !!!
44nts=4  & nte=nts  & nxs=31  & nxe=nxs ;; AU CENTRE ! sinon choisir le point en traçant de 0a 59
45nts=1  & nte=nts  & nxs=0   & nxe=59
46nts=4  & nte=nts  & nxs=0   & nxe=59
47nts=1  & nte=15   & nxs=31  & nxe=31
48nts=1  & nte=15   & nxs=45  & nxe=45
49nts=9  & nte=9    & nxs=0   & nxe=59
50nts=1  & nte=13   & nxs=59  & nxe=59
51nts=8  & nte=8    & nxs=0   & nxe=121
52nts=1  & nte=9    & nxs=59  & nxe=59
53;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
54nts=1  & nte=10  & nxs=50   & nxe=50
55;nts=1  & nte=1  & nxs=6   & nxe=6
56nts=1  & nte=100  & nxs=50   & nxe=50
57;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
58;
59;
60;
61what_I_plot=0. & overcontour=0.
62;SPAWN, '\rm param_plot.idl ; cp gravitwave_inc2.pro param_plot.idl'
63;
64;
65;
66getcdf, $
67        file=file, $
68        charvar=charvar, $
69        invar=invar
70;getcdf, $
71;        file=file, $
72;        charvar=charvarc, $
73;        invar=invarc
74getcdf, $
75        file=file, $
76        charvar='vert', $
77        invar=vert
78getcdf, $
79        file=file, $
80        charvar='PTOT', $
81        invar=columnp
82;;
83;;
84;;
85        s = size(invar) & middle = floor(s[2]/2) & print, 'PLOT at y subs ', middle
86        nte = min(s[4]-1,nte) & print, 'nte ', nte
87        mean_invar = TOTAL(invar(*,*,*,nts:nte),4) / float(nte-nts+1)
88        ;mean_columnp = TOTAL(columnp(*,*,*,nts:nte),4) / float(nte-nts+1)
89        anomal_invar = invar & for i=0,s[4]-1 do anomal_invar(*,*,*,i) = anomal_invar(*,*,*,i) - mean_invar
90;
91;
92;
93        ;goto, perturb
94for nt=nts,nte do begin
95for nx=nxs,nxe do begin
96;
97zefile=experiment+'_'+charvar+'_'+string(nt+100,'(I0)')+'_'+string(nx+100,'(I0)')
98;
99;
100;
101  PS_Start, filename=zefile+'.ps'
102  print, zefile+'.ps'
103  !P.Charsize = 1.2
104  !p.charthick = 2.0
105  !p.thick = 2.0
106  !x.thick = 2.0
107  !y.thick = 2.0
108;
109;
110;
111;;what_I_plot = reform(invar(*,1,*,nt))
112;what_I_plot = reform(invar(*,59,*,nt))
113what_I_plot = reform(invar(*,middle,*,nt))
114
115   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
116   file=path+experiment+'/input_sounding' & header='' & nlines_header=1 & ncol = 5
117   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
118   nlines = FILE_LINES(file)-nlines_header & data=FLTARR(ncol,nlines)
119   OPENR, lun, file, /GET_LUN & READF, lun, header & READF, lun, data & CLOSE, lun & FREE_LUN, lun
120   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
121   mcd_column      = reform(data(0,*))
122   calc_tpot       = reform(data(1,*))
123
124   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
125   file=path+experiment+'/input_therm' & ncol = 5
126   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
127   nlines = FILE_LINES(file) & data=FLTARR(ncol,nlines)
128   OPENR, lun, file, /GET_LUN & READF, lun, data & CLOSE, lun & FREE_LUN, lun
129   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
130   pressure        = reform(data(2,*))
131   mcd_tpot        = reform(data(4,*))   
132
133   if (nx lt 999) then begin
134     what_I_plot = reform(what_I_plot(nx,*))
135   endif else begin
136     what_I_plot = TOTAL(what_I_plot,1) / n_elements(reform(what_I_plot(*,0)))
137   endelse
138
139     ;;;;;;;;;;;
140     ;;;;;;;;;;;
141     ;;;;;;;;;;;
142     column = vert/1000.
143     overplot_column = mcd_column/1000.
144     alt = [20.,150.]
145     alt = [40.,100.]
146     alt = [40.,80.]
147     alt = [0.,120.]
148     ;;;;;;;;;;;
149     ;;;;;;;;;;;
150     ;mention = ''
151     SPAWN, '\rm param_plot.idl'
152     minfield_init = 85.  ;100.
153     maxfield_init = 185. ;250.
154     title_axis = ['Temperature (K)','Altitude (km)']
155     overplot = mcd_tpot
156     ;;;;;;;;;;;
157     SPAWN, "echo 'intervaly=10.' >> param_plot.idl"
158     SPAWN, "echo 'intervalx=5.' >> param_plot.idl"
159     ;;;;;;;;;;;
160
161        ;overplot = reform(mean_invar(nx,middle,*))
162        ;overplot_column = column
163
164;;;;;;;;;;;;;;;;;;;;;;;;;;;; PLOT PRESSION mettre profile en log
165;column = reform(columnp(nx,middle,*,nt))
166;overplot_column = pressure
167;       ;;overplot = mean_invar
168;       ;;overplot_column = mean_columnp
169;alt = [1.e2,1.e-3]
170;title_axis = ['Temperature (K)','Pressure (Pa)']
171;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PLOT PRESSION
172
173;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
174;; calculer la CAPE ? Tsat - Tenv / Tenv
175
176
177                        ;what_I_plot = pressure
178                        ;       what_I_plot = calc_tpot
179                        ;overplot = alog(overplot) / alog(10.)
180                        ;overplot = overplot[where(overplot lt 35.)] ;; evite les valeurs impossibles
181                        ;what_I_plot = alog(what_I_plot) / alog(10.)
182                        ;minfield_init = -5.
183                        ;maxfield_init = 3.
184                        ;       ;;tpot
185                        ;       minfield_init = 2.
186                        ;       maxfield_init = 4.
187                        ;what_I_plot = 100. * abs(calc_tpot - overplot) / overplot
188                        ;minfield_init = 0.000001
189                        ;maxfield_init = 100.
190
191   profile, $
192        what_I_plot, $                          ; 1D vertical profile
193        column, $                               ; altitudes     
194        alt=alt, $                              ; altitude range [altmin, altmax]
195        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
196        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
197        inprofile=overplot, $                   ; another vertical profile to overplot
198        incolumn=overplot_column, $             ; altitudes of the other vertical profile (in case /= column)
199;        discrete=discrete, $                    ; show the profile points (= type of points in !psym)
200;        title_plot=title_user, $                ; title of the plot ('Profile' is default)
201        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
202        mention=mention                         ; add text precision within the plot window (default is nothing or '')
203
204        ;;; sponge layer
205        ;w=where(FINITE(overplot) ne 0)
206        ;oplot, [minfield_init,maxfield_init], [max(overplot_column[w])-50.,max(overplot_column[w])-50.], linestyle=1
207
208        ;column = columnp
209        ;@tempcond.inc
210        ;yeye = reform(overplot(nx,1,*,nt))
211        ;;w = where(yeye le 0.) & yeye[w] = !VALUES.F_NAN
212        ;;oplot, yeye, vert/1000., linestyle=1
213
214@tempcond.inc
215overplot = reform(overplot(nx,middle,*,nt))
216oplot, overplot, overplot_column, linestyle=1
217
218        ;calc_tk = calc_tpot * (pressure/610.)^(1/4.4)
219        ;oplot, calc_tk, column, linestyle=3
220
221PS_End, /PNG
222endfor
223endfor
224
225perturb:
226PS_Start, filename=experiment+'_'+charvar+'_perturb_'+string(nxs+100,'(I0)')+'.ps'
227
228     what_I_plot_loop = reform(anomal_invar(nxs,middle,*,*))
229
230     what_I_plot = reform(what_I_plot_loop(*,nts))
231     column = vert/1000.
232     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
233     alt = [0.,120.]
234     minfield_init = -30.
235     maxfield_init = 30.
236     title_axis = ['Temperature anomaly (K)','Altitude (km)']
237                ;; peut pas changer     
238                ;SPAWN, "echo 'intervaly=10.' >> param_plot.idl"
239                ;SPAWN, "echo 'intervalx=5.' >> param_plot.idl"
240     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
241     profile, $
242        what_I_plot, $                          ; 1D vertical profile
243        column, $                               ; altitudes     
244        alt=alt, $                              ; altitude range [altmin, altmax]
245        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
246        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
247        inprofile=overplot, $                   ; another vertical profile to overplot
248        incolumn=overplot_column, $             ; altitudes of the other vertical profile (in case /= column)
249;        discrete=discrete, $                    ; show the profile points (= type of points in !psym)
250;        title_plot=title_user, $                ; title of the plot ('Profile' is default)
251        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
252        mention=mention                         ; add text precision within the plot window (default is nothing or '')
253
254     for iii=nts+1, nte do begin
255        overplot = reform(what_I_plot_loop(*,iii))
256        overplot_column = column
257        oplot, overplot, overplot_column
258     endfor
259PS_End, /PNG
260
261end
Note: See TracBrowser for help on using the repository browser.