source: trunk/UTIL/IDLplot/SPEC/GW/gravitwaveprof.pro @ 1185

Last change on this file since 1185 was 324, checked in by aslmd, 13 years ago

MESOSCALE : Preparatory commit for the ultimate option mars=42 which

would allow mesoscale modeling with photochemistry.

[see r76 method to add 'mars' options]
Modified module_lmd_driver.F and Registry.EM and runmeso
Modified readmeteo.F90 and introduced an option -DPHOTOCHEM
...
Transparent to the casual user
Option mars=42 not yet finished though -- so do not use!
  • Property svn:executable set to *
File size: 9.9 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='/home/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'
30;experiment='GW_MARS_highwind_3D_loweropacity_morepoints_LT15.161544'
31;experiment='GW_MARS_highwind_3D_loweropacity_morepoints_LTcst.168440'
32;experiment='GW_MARS_highwind_3D_loweropacity_morepoints_widehill.168468'
33        file=path+experiment+'_wrfout_d01_9999-09-09_09:00:00_z'
34file=path+experiment+'/wrfout_d01_9999-09-09_09:00:00_z'
35;file=path+experiment+'/wrfout_d01_9999-09-09_10:00:00_z'
36;file=path+experiment+'/wrfout_d01_9999-09-09_11:00:00_z'
37charvar='W' & charvarc='W'
38charvar='tk' & charvarc='W'
39charvar='tk' & charvarc='W' & cond=1
40charvar='tpot' & charvarc='tpot'
41;charvar='PTOT' & charvarc='PTOT'
42charvar='tk' & charvarc='tk'
43;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
44nts=1  & nte=15   & nxs=0   & nxe=59
45nts=3  & nte=nts  & nxs=999 & nxe=999 ;; moyenne ;; attention a la montagne !!!
46nts=4  & nte=nts  & nxs=31  & nxe=nxs ;; AU CENTRE ! sinon choisir le point en traçant de 0a 59
47nts=1  & nte=nts  & nxs=0   & nxe=59
48nts=4  & nte=nts  & nxs=0   & nxe=59
49nts=1  & nte=15   & nxs=31  & nxe=31
50nts=1  & nte=15   & nxs=45  & nxe=45
51nts=9  & nte=9    & nxs=0   & nxe=59
52nts=1  & nte=13   & nxs=59  & nxe=59
53nts=8  & nte=8    & nxs=0   & nxe=121
54nts=1  & nte=9    & nxs=59  & nxe=59
55;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
56nts=1  & nte=10  & nxs=50   & nxe=50
57;nts=1  & nte=1  & nxs=6   & nxe=6
58nts=1  & nte=100  & nxs=50   & nxe=50
59nts=1  & nte=100  & nxs=60   & nxe=60
60;nts=1  & nte=100  & nxs=20   & nxe=20
61nts=29 & nte=29   & nxs=60   & nxe=60
62;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
63;
64;
65;
66what_I_plot=0. & overcontour=0.
67;SPAWN, '\rm param_plot.idl ; cp gravitwave_inc2.pro param_plot.idl'
68;
69;
70;
71getcdf, $
72        file=file, $
73        charvar=charvar, $
74        invar=invar
75;getcdf, $
76;        file=file, $
77;        charvar=charvarc, $
78;        invar=invarc
79getcdf, $
80        file=file, $
81        charvar='vert', $
82        invar=vert
83getcdf, $
84        file=file, $
85        charvar='PTOT', $
86        invar=columnp
87;;
88;;
89;;
90        s = size(invar) & middle = floor(s[2]/2) & print, 'PLOT at y subs ', middle
91        nte = min(s[4]-1,nte) & print, 'nte ', nte
92        mean_invar = TOTAL(invar(*,*,*,nts:nte),4) / float(nte-nts+1)
93        ;mean_columnp = TOTAL(columnp(*,*,*,nts:nte),4) / float(nte-nts+1)
94        anomal_invar = invar & for i=0,s[4]-1 do anomal_invar(*,*,*,i) = anomal_invar(*,*,*,i) - mean_invar
95;;
96;;
97;;
98;       goto, perturb
99for nt=nts,nte do begin
100for nx=nxs,nxe do begin
101;
102zefile=experiment+'_'+charvar+'_'+string(nt+100,'(I0)')+'_'+string(nx+100,'(I0)')
103;
104;
105;
106  PS_Start, filename=zefile+'.ps'
107  print, zefile+'.ps'
108  !P.Charsize = 1.2
109  !p.charthick = 2.0
110  !p.thick = 2.0
111  !x.thick = 2.0
112  !y.thick = 2.0
113;
114;
115;
116;;what_I_plot = reform(invar(*,1,*,nt))
117;what_I_plot = reform(invar(*,59,*,nt))
118what_I_plot = reform(invar(*,middle,*,nt))
119
120   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
121   file=path+experiment+'/input_sounding' & header='' & nlines_header=1 & ncol = 5
122   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
123   nlines = FILE_LINES(file)-nlines_header & data=FLTARR(ncol,nlines)
124   OPENR, lun, file, /GET_LUN & READF, lun, header & READF, lun, data & CLOSE, lun & FREE_LUN, lun
125   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
126   mcd_column      = reform(data(0,*))
127   calc_tpot       = reform(data(1,*))
128
129   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
130   file=path+experiment+'/input_therm' & ncol = 5
131   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
132   nlines = FILE_LINES(file) & data=FLTARR(ncol,nlines)
133   OPENR, lun, file, /GET_LUN & READF, lun, data & CLOSE, lun & FREE_LUN, lun
134   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
135   pressure        = reform(data(2,*))
136   mcd_tpot        = reform(data(4,*))   
137
138   if (nx lt 999) then begin
139     what_I_plot = reform(what_I_plot(nx,*))
140   endif else begin
141     what_I_plot = TOTAL(what_I_plot,1) / n_elements(reform(what_I_plot(*,0)))
142   endelse
143
144     ;;;;;;;;;;;
145     ;;;;;;;;;;;
146     ;;;;;;;;;;;
147     column = vert/1000.
148     overplot_column = mcd_column/1000.
149     alt = [20.,150.]
150     alt = [40.,100.]
151     alt = [40.,80.]
152     alt = [0.,120.]
153     alt = [10.,110.]
154     ;;;;;;;;;;;
155     ;;;;;;;;;;;
156     ;mention = ''
157     SPAWN, '\rm param_plot.idl'
158     minfield_init = 85.  ;100.
159     maxfield_init = 185. ;250.
160     minfield_init = 90.
161     maxfield_init = 230.
162     title_axis = ['Temperature (K)','Altitude (km)']
163     overplot = mcd_tpot
164     ;;;;;;;;;;;
165     SPAWN, "echo 'intervaly=10.' >> param_plot.idl"
166     SPAWN, "echo 'intervalx=5.' >> param_plot.idl"
167     SPAWN, "echo 'intervalx=10.' >> param_plot.idl"
168     ;;;;;;;;;;;
169
170        ;overplot = reform(mean_invar(nx,middle,*))
171        ;overplot_column = column
172
173;;;;;;;;;;;;;;;;;;;;;;;;;;;; PLOT PRESSION mettre profile en log
174;column = reform(columnp(nx,middle,*,nt))
175;overplot_column = pressure
176;       ;;overplot = mean_invar
177;       ;;overplot_column = mean_columnp
178;alt = [1.e2,1.e-3]
179;title_axis = ['Temperature (K)','Pressure (Pa)']
180;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PLOT PRESSION
181
182;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
183;; calculer la CAPE ? Tsat - Tenv / Tenv
184
185
186                        ;what_I_plot = pressure
187                        ;       what_I_plot = calc_tpot
188                        ;overplot = alog(overplot) / alog(10.)
189                        ;overplot = overplot[where(overplot lt 35.)] ;; evite les valeurs impossibles
190                        ;what_I_plot = alog(what_I_plot) / alog(10.)
191                        ;minfield_init = -5.
192                        ;maxfield_init = 3.
193                        ;       ;;tpot
194                        ;       minfield_init = 2.
195                        ;       maxfield_init = 4.
196                        ;what_I_plot = 100. * abs(calc_tpot - overplot) / overplot
197                        ;minfield_init = 0.000001
198                        ;maxfield_init = 100.
199
200   profile, $
201        what_I_plot, $                          ; 1D vertical profile
202        column, $                               ; altitudes     
203        alt=alt, $                              ; altitude range [altmin, altmax]
204        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
205        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
206        inprofile=overplot, $                   ; another vertical profile to overplot
207        incolumn=overplot_column, $             ; altitudes of the other vertical profile (in case /= column)
208;        discrete=discrete, $                    ; show the profile points (= type of points in !psym)
209;        title_plot=title_user, $                ; title of the plot ('Profile' is default)
210        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
211        mention=mention                         ; add text precision within the plot window (default is nothing or '')
212
213        ;;; sponge layer
214        ;w=where(FINITE(overplot) ne 0)
215        ;oplot, [minfield_init,maxfield_init], [max(overplot_column[w])-50.,max(overplot_column[w])-50.], linestyle=1
216
217        ;column = columnp
218        ;@tempcond.inc
219        ;yeye = reform(overplot(nx,1,*,nt))
220        ;;w = where(yeye le 0.) & yeye[w] = !VALUES.F_NAN
221        ;;oplot, yeye, vert/1000., linestyle=1
222
223
224loadct, 4
225for i=16,37 do begin  ;10-37
226for j=0,40 do begin
227dec = j - 20 
228overoverplot = reform(invar(nx+dec,middle+dec,*,i))
229oplot, overoverplot, column, color=220  ;psym=3
230endfor
231endfor
232
233!p.thick = 4.0
234oplot, overplot, overplot_column, linestyle=2
235!p.thick = 2.0
236
237!p.thick = 4.0
238oplot, reform(invar(nx,middle,*,28)), column
239!p.thick = 2.0
240
241@tempcond.inc
242overplot = reform(overplot(nx,middle,*,nt))
243w = where(overplot le 0.) & overplot[w] = !VALUES.F_NAN
244oplot, overplot, overplot_column, linestyle=1
245
246        ;calc_tk = calc_tpot * (pressure/610.)^(1/4.4)
247        ;oplot, calc_tk, column, linestyle=3
248
249PS_End, /PNG
250endfor
251endfor
252
253perturb:
254PS_Start, filename=experiment+'_'+charvar+'_perturb_'+string(nxs+100,'(I0)')+'.ps'
255
256     what_I_plot_loop = reform(anomal_invar(nxs,middle,*,*))
257
258     what_I_plot = reform(what_I_plot_loop(*,nts))
259     column = vert/1000.
260     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
261     alt = [0.,120.]
262     minfield_init = -30.
263     maxfield_init = 30.
264     title_axis = ['Temperature anomaly (K)','Altitude (km)']
265                ;; peut pas changer     
266                ;SPAWN, "echo 'intervaly=10.' >> param_plot.idl"
267                ;SPAWN, "echo 'intervalx=5.' >> param_plot.idl"
268     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
269     profile, $
270        what_I_plot, $                          ; 1D vertical profile
271        column, $                               ; altitudes     
272        alt=alt, $                              ; altitude range [altmin, altmax]
273        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
274        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
275        inprofile=overplot, $                   ; another vertical profile to overplot
276        incolumn=overplot_column, $             ; altitudes of the other vertical profile (in case /= column)
277;        discrete=discrete, $                    ; show the profile points (= type of points in !psym)
278;        title_plot=title_user, $                ; title of the plot ('Profile' is default)
279        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
280        mention=mention                         ; add text precision within the plot window (default is nothing or '')
281
282     for iii=nts+1, nte do begin
283        overplot = reform(what_I_plot_loop(*,iii))
284        overplot_column = column
285        oplot, overplot, overplot_column
286     endfor
287
288        ;temppp = reform(mean_invar(nxs,middle,*))
289        ;hache = 191.*temppp/3.72/1000.
290        ;ok = column / ( 2. * hache )
291        ;print, exp(ok)
292        ;yeah = exp(ok)
293        ;oplot, yeah, column
294
295PS_End, /PNG
296
297end
Note: See TracBrowser for help on using the repository browser.