source: trunk/mesoscale/LMD_MM_MARS/SRC/ARWpost/idl/graphwrf.pro @ 13

Last change on this file since 13 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 18.2 KB
Line 
1pro graphwrf, typeplot
2
3
4;-----------------------------------------
5; Read files generated by ARWpost & plot
6;
7; typeplot is s for see         (.ps are not saved, .def saved)
8;             t for trace       (.ps are saved, .def saved)
9;             r for report      (same as t + tex report)
10;             m for movie       (same as t + gif movie)
11;             w for web page    (same as t + html animation)
12;
13; default is 's'
14;
15; use:
16;  1. set parameters in graphwrf.def
17;  2. @graphwrf.comp
18;  3. graphwrf, 's' (or 't' or 'r' or 'm' or 'w')
19;
20;-----------------------------------------
21; A. Spiga, Mars 2008
22;-----------------------------------------
23
24
25nomovie='yes'
26
27if (n_elements(typeplot) eq 0) then typeplot='s'
28
29print, '---------------------'
30print, 'init. please wait ...'
31print, '---------------------'
32
33;
34; init.
35;
36field1='' & field2='' & winds=['','']
37backup_data='no' & already_data='no'
38datafolder='./' & plotfolder='./'
39textitle='plot'
40texcomments='plot'
41topo=0
42extract='no'
43
44
45;
46; defining fields, coordinates, type of plot ... etc ...
47;
48@user.idl
49
50
51;
52; extracting/interpolating fields from WRF NETCDF files
53;
54if (extract eq 'yes') then begin
55        print, 'EXTRACT DATA. OK ?'
56        read, dumb, prompt='type 1 to proceed, 0 to keep the same input >>> '
57        if (dumb ne 0) then begin
58                print, 'WE WILL EXTRACT DATA USING ARWPOST.EXE'
59                if (n_elements(nest) eq 0) then nest='d01'
60                call_arwpost, nam1, nam2, nam3, nam4, nam5, tabnam, nest, datafolder
61        endif
62endif
63if (already_data eq 'yes') then begin
64        if (nam5 eq " interp_method = 0 ") then coord='model_level'
65        if ((nam5 eq " interp_method = 1 ") and (tabnam(0) lt tabnam(1))) then coord='height'
66        SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.ctl input.ctl'
67        SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.dat input.dat'
68endif
69
70                               
71;-----------------------------;
72; STARTING TO READ INPUT.CTL  ;
73;-----------------------------;
74
75
76;
77; get the local time
78;
79interv=0.
80openr, 99, 'timefil'
81readf,99,interv
82close, 99
83
84
85;
86; get preliminar info on coordinates
87;
88wnx=1
89wny=1
90wnz=1
91wnt=1
92txt=''
93param=[0.,0.,0.,0.,0.,0.,0.,0.]
94
95openr, 1, 'input.ctl'
96while not eof(1) do begin
97    readf,1,txt
98    match=STRPOS(txt,'linear')
99        if (match ne -1) then begin
100                wnx=(STRPOS(txt,'xdef') eq -1)*wnx  ;0 if linear
101                wny=(STRPOS(txt,'ydef') eq -1)*wny  ;0 if linear
102                wnz=(STRPOS(txt,'zdef') eq -1)*wnz  ;0 if linear
103                wnt=(STRPOS(txt,'tdef') eq -1)*wnt  ;0 if linear
104                paramstart=STRSPLIT(txt, 'r', LENGTH=paramlength)
105                paramchar=STRMID(txt,paramstart(1),paramlength(1))
106                READS, paramchar, start, step
107                param=param+[start*(STRPOS(txt,'xdef') ne -1),$
108                        step*(STRPOS(txt,'xdef') ne -1),$
109                        start*(STRPOS(txt,'ydef') ne -1),$
110                        step*(STRPOS(txt,'ydef') ne -1),$
111                        start*(STRPOS(txt,'zdef') ne -1),$
112                        step*(STRPOS(txt,'zdef') ne -1),$
113                        start*(STRPOS(txt,'tdef') ne -1),$
114                        step*(STRPOS(txt,'tdef') ne -1)]
115        endif   
116endwhile
117close, 1
118timebegin=param(6)
119       
120
121;---------------------------------
122; READ .CTL FILE INFORMATIONS
123;---------------------------------
124
125info=read_ascii('input.ctl', missing_value=1e37)
126infodat=info.field1
127
128;
129; second column : dimensions
130;
131
132w=where((infodat(0,*) eq 1e37) and (infodat(1,*) lt 1e37))
133dimensions=reform(infodat(1,w))
134
135nx=round(dimensions(0))
136ny=round(dimensions(1))
137nz=round(dimensions(2))
138nt=round(dimensions(3))
139fields=round(dimensions(4))
140vertical=intarr(fields)
141for i=0,fields-1 do begin
142        vertical(i)=round(dimensions(5+i))
143endfor
144
145
146;
147; first column : coordinates
148;
149
150w=where((infodat(1,*) eq 1e37) and (infodat(0,*) lt 1e37))
151if (w(0) ne -1) then coordinates=reform(infodat(0,w))
152
153
154if (wnx ne 0) then begin
155
156        wnx=nx
157        x=coordinates(0:nx-1)
158
159endif else begin
160
161        xmin=param(0)
162        xmax=param(0)+param(1)*(nx-1) MOD 360
163        ;inc=param(0)+param(1)*(nx-1)   ;cas mercator
164        ;x=start+inc*findgen(nx)
165        x=xmin+(xmax-xmin)*findgen(nx)/(nx-1)
166       
167;param(1)=param(1)-param(0) ;;bricolo plot GCM
168;xmin=param(0)
169;xmax=param(0)+param(1)*(nx-1)
170;x=xmin+(xmax-xmin)*findgen(nx)/(nx-1)
171;;x = x - (LONG(x)/180)*360.0
172;print, x
173
174        ;stop
175
176endelse
177
178
179if (wny ne 0) then begin
180        wny=ny
181        y=coordinates(wnx:wnx+wny-1)
182endif else begin
183        start=param(2)
184        step=param(3)
185        y=start+step*findgen(ny)
186endelse
187if (wnz ne 0) then begin
188        wnz=nz
189        z=coordinates(wnx+wny:wnx+wny+wnz-1)
190endif else begin
191        start=param(4)
192        step=param(5)
193        z=start+step*findgen(nz)
194endelse
195t=interv*findgen(nt)/3700. + timebegin + mean(x)/15.
196lct = round(t+24) MOD 24       
197       
198;
199; get info on fields
200;
201openr, 1, 'input.ctl'
202toto=5+wnx+1+wny+1+wnz+2+fields
203infofields=strarr(toto)
204readf, 1, infofields
205close, 1
206desc_field=transpose(infofields(toto-fields:toto-1))
207
208;
209; vert coord
210;
211if (z(0) lt z(1)) then coord='height'
212if (z(0) gt z(1)) then coord='pressure'
213if ((z(0) eq 1.) and (z(nz-1) eq (z(nz-2)+1.))) then coord='model_level'
214
215
216
217;---------------------------------
218; READ .DAT FILE DATA
219;---------------------------------
220; GrADS views gridded data sets as 5-dimensional arrays
221; varying in longitude, latitude, vertical level, variable, and time.
222;---------------------------------
223
224;-----------------------------;
225; STARTING TO READ INPUT.DATA ;
226;-----------------------------;
227print, '----------------'
228print, 'reading data ...'
229print, '----------------'
230print, 'longitudes: ', min(x), max(x), nx, ' points'
231print, 'latitudes: ', min(y), max(y), ny, ' points'
232print, 'altitudes: ', min(z), max(z), nz, ' points'
233print, 'local times: ', min(lct), max(lct), nt, ' points'
234print, fields, ' recorded fields :'
235print, desc_field
236print, '----------------'
237
238; see http://www.dfanning.com/tips/ascii_column_data.html       
239OPENR, lun, 'input.dat', /GET_LUN
240
241
242;
243; titles and plot name
244;
245denom=plot+'_'+coord+'_'+field1
246if ((field2 ne '') and (topo eq 0)) then denom=denom+'_'+field2
247if (topo eq 1) then denom=denom+'_HGT'
248if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
249        case coord of
250        'height': charlev='Altitude (km)'
251        ;'pressure': begin
252        ;               H=10.
253        ;               p0=6.1
254        ;               z=H*alog(p0/z)
255        ;               charlev='Log-Pressure altitude (km)'
256        ;end
257        'model_level': charlev='Model level'
258        endcase
259                                                                                                                                                                                       
260
261;
262; graphics settings
263;
264!p.charthick = 2.0
265!p.thick = 3.0
266!x.thick = 2.0
267!y.thick = 2.0
268
269        ; plot with topography
270        if (topo eq 1) then field2='HGT'
271
272save_ps='true'
273if (save_ps ne 'false') then begin
274;       set_plot, 'ps'
275endif else begin
276        set_plot, 'x'
277        !P.MULTI = [0, 2, 3]
278        window, 0, xsize=600, ysize=900
279endelse
280
281no3D='false'   
282print, 'GENERATE PLOTS ',min([nt-1,num-1])+1
283
284for l=0,min([nt-1,num-1]) do begin                                              ;time loop
285print, 'time loop ...', l+1, ' / ', string(nt,'(I0)') 
286T = SYSTIME(1)
287  for f=0, fields-1 do begin ;fields (whether 2D or 3D) ;variable loop
288  nzf=vertical(f)
289  print, 'scanning field ...', f+1, ' / ', string(fields,'(I0)')
290
291        data = FLTARR(1,nzf*ny*nx)
292        print, 'read data ...'
293        READF, lun, data        ;; the IDL way, two times faster than 3 FOR loops
294        print, 'read data done'
295       
296        recognize=STRSPLIT(desc_field(f), /EXTRACT)
297        recognize=recognize[0]
298       
299;       IF (STRPOS(desc_field(f),field1) ne -1) then begin
300;       IF ( ( STRPOS(desc_field(f),field1) eq 0 ) and ( STRLEN(desc_field(f)) eq STRLEN(field1) ) ) then begin
301        IF (recognize eq field1) then begin
302                data_save=data
303                print, 'found field1 ', field1
304                data=REFORM(data,nx,ny,nzf,/OVERWRITE)
305                        w=where(data eq -9999.) ; missing values: -9999 (before, was 1e20)
306                        if (w(0) ne -1) then data[w]=1e37
307                CASE (plot) OF
308                'map': what_I_plot=reform(data(*,*,-1+min([level,nzf])))
309                'meridional': what_I_plot=reform(data(-1+min([nlon,nx]),*,*))
310                'zonal': what_I_plot=reform(data(*,-1+min([nlat,ny]),*))
311                ENDCASE
312                if (nzf eq 1) then no3D='true'
313                data=0. & data=data_save & data_save=0.
314        ENDIF
315;       IF ((STRPOS(desc_field(f),field2) ne -1) and (field2 ne '')) then begin
316;       IF ((STRPOS(desc_field(f),field2) eq 0) and (field2 ne '')) then begin
317;       IF ( ( STRPOS(desc_field(f),field2) eq 0 ) and ( field2 ne '' ) and ( STRLEN(desc_field(f)) eq STRLEN(field2) ) ) then begin
318        IF (recognize eq field2) then begin
319                data_save=data
320                print, 'found field2 ', field2
321                if (topo eq 1) then begin
322                        data=data/100.
323                        data=round(data)
324                        data=float(data)/10.
325                endif   
326                data=REFORM(data,nx,ny,nzf,/OVERWRITE)
327                        w=where(data eq -9999.)   ; missing values
328                        if (w(0) ne -1) then data[w]=1e37
329                CASE (plot) OF
330                'map': overcontour=reform(data(*,*,-1+min([level,nzf])))
331                'meridional': overcontour=reform(data(-1+min([nlon,nx]),*,*))
332                'zonal': overcontour=reform(data(*,-1+min([nlat,ny]),*))
333                ENDCASE
334                data=0. & data=data_save & data_save=0.
335        ENDIF   
336;       IF (STRPOS(desc_field(f),winds(0)) ne -1 and (winds(0) ne '')) then begin
337;       IF (STRPOS(desc_field(f),winds(0)) eq 0 and (winds(0) ne '')) then begin
338        IF (recognize eq winds(0)) then begin
339                data_save=data
340                print, 'found wind1 ', winds(0)
341                data=REFORM(data,nx,ny,nzf,/OVERWRITE)
342                        w=where(data eq -9999.)   ; missing values
343                        if (w(0) ne -1) then data[w]=!Values.F_NAN
344                CASE (plot) OF
345                'map': overvector_x=reform(data(*,*,-1+min([level,nzf])))
346                'meridional': overvector_x=reform(data(-1+min([nlon,nx]),*,*))
347                'zonal': overvector_x=reform(data(*,-1+min([nlat,ny]),*))
348                ENDCASE
349                data=0. & data=data_save & data_save=0.
350        ENDIF
351;       IF (STRPOS(desc_field(f),winds(1)) ne -1 and (winds(1) ne '')) then begin
352;       IF (STRPOS(desc_field(f),winds(1)) eq 0 and (winds(1) ne '')) then begin
353        IF (recognize eq winds(1)) then begin
354                data_save=data
355                print, 'found wind2 ', winds(1)
356                data=REFORM(data,nx,ny,nzf,/OVERWRITE)
357                        w=where(data eq -9999.)   ; missing values
358                        if (w(0) ne -1) then data[w]=!Values.F_NAN
359                CASE (plot) OF
360                'map': overvector_y=reform(data(*,*,-1+min([level,nzf])))
361                'meridional': overvector_y=reform(data(-1+min([nlon,nx]),*,*))
362                'zonal': overvector_y=reform(data(*,-1+min([nlat,ny]),*))
363                ENDCASE
364                data=0. & data=data_save & data_save=0.
365        ENDIF
366        data=0.         ;; free memory
367        data_save=0.         ;; free memory
368       
369  endfor
370
371        if (n_elements(overcontour) eq 0) then overcontour=0.
372        if (n_elements(overvector_x) eq 0) then overvector_x=0.
373        if (n_elements(overvector_y) eq 0) then overvector_y=0.
374       
375        @graphwrf_name.idl
376
377;;--------------------------------------------------------
378;; NB: you can use the winds fields to combine variables
379;;--------------------------------------------------------
380@wind_stress.idl        ;** WIND STRESS
381@sea_level.idl          ;** HYDROSTATIC REDUCTION TO LEVEL OF REFERENCE
382@ertel_vorticity.idl    ;** ERTEL POTENTIAL VORTICITY
383@tau_cloud.idl          ;** CLOUD OPACITY
384@wind_velocity.idl      ;** 3D WIND VELOCITY
385@hwind_velocity.idl     ;** HORIZONTAL WIND VELOCITY
386;;--------------------------------------------------------
387;; NB: you can use the winds fields to combine variables
388;;--------------------------------------------------------
389
390        title=title+', LT '+string(lct(l),'(I0)')
391
392        if (save_ps ne 'false') then begin
393        save_ps=denom+string(1000+l,'(I0)')
394;        device, filename=save_ps+'.ps', /color
395PS_Start, filename=save_ps+'.ps'
396        endif
397
398        CASE (plot) OF
399        'map': begin
400                case coord of
401                   'height': charlev=string(round(z(level-1)*1000.),'(I0)')+' m'
402                   'pressure': charlev=string(round(z(level-1)*100.),'(I0)')+' Pa'
403                   'model_level': charlev=string(round(z(level-1)),'(I0)')
404                endcase
405                if (no3D ne 'true') then title=title+', at level '+charlev
406                print, title
407        map_latlon, $
408                what_I_plot, $                          ; 2D field
409                x, $                                  ; 1D latitude
410                y, $                                  ; 1D longitude
411                minfield=0., $               ; minimum value of plotted field (=0: calculate)
412                maxfield=0., $               ; maximum value of plotted field (=0: calculate)
413                overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
414                overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
415                overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
416                ct=pal, $                               ; color table (33-rainbow is default)
417                colors=colors, $                        ; number of colors/levels (32 is default)
418                title=title, $                     ; title of the plot ('' is default)
419                format=format                           ; format of colorbar annotations ('(F6.2)' is default)
420        end
421        'zonal': begin
422                displaylat=float(round(y(nlat-1)*10))/10.
423                title=title+ ', section at latitude '+string(displaylat,format='(1F6.1)')+ '!Uo!N '
424                print, title
425        section, $
426                what_I_plot, $                          ; 2D field
427                x, $                                ; horizontal coordinate
428                z, $                             ; altitude coordinate
429                minfield=0., $               ; minimum value of plotted field (=0: calculate)
430                maxfield=0., $               ; maximum value of plotted field (=0: calculate)
431                minspace=0., $                    ; minimum value of space window (=0: calculate)
432                maxspace=0., $                    ; maximum value of space window (=0: calculate)
433                overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
434                overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
435                overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
436                colors=colors, $                        ; number of colors/levels (32 is default)
437                title_plot=title, $                ; title of the plot ('Profile' is default)
438                title_axis=['Longitude',charlev], $             ; title of the [x,y] axis (['Field','Altitude'] is default)
439                ct=pal, $                               ; color table (33-rainbow is default)
440                format=format
441        end
442        'meridional': begin
443                displaylon=float(round(x(nlon-1)*10))/10.
444                title=title+', section at longitude '+string(displaylon,format='(1F6.1)')+ '!Uo!N '
445                print, title
446        section, $
447                what_I_plot, $                          ; 2D field
448                y, $                                ; horizontal coordinate
449                z, $                             ; altitude coordinate
450                minfield=0., $               ; minimum value of plotted field (=0: calculate)
451                maxfield=0., $               ; maximum value of plotted field (=0: calculate)
452                minspace=0., $                    ; minimum value of space window (=0: calculate)
453                maxspace=0., $                    ; maximum value of space window (=0: calculate)
454                overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
455                overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
456                overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
457                colors=colors, $                        ; number of colors/levels (32 is default)
458                title_plot=title, $                ; title of the plot ('Profile' is default)
459                title_axis=['Latitude',charlev], $             ; title of the [x,y] axis (['Field','Altitude'] is default)
460                ct=pal, $                               ; color table (33-rainbow is default)
461                format=format
462        end
463        endcase 
464
465         if (save_ps ne 'false') then device, /close
466         if (typeplot eq 's') then SPAWN, 'display '+save_ps+'.ps &'
467
468;PRINT, 1000*(SYSTIME(1) - T), ' milliSeconds'
469PRINT, SYSTIME(1) - T, ' seconds to process'
470endfor
471
472
473
474CLOSE, lun
475print, '...done !'
476help, /memory
477print, '---------'
478
479
480CASE (typeplot) OF
481's': begin
482        ;SPAWN, 'find '+denom+'????.ps -exec display {} \;'
483        SPAWN, 'cp -f graphwrf.def '+plotfolder+'/see_'+denom+'.def'
484        ;SPAWN, 'mv -f ./'+denom+'????.ps '+plotfolder+'/'
485        SPAWN, '\rm ./'+denom+'????.ps'
486end
487't': begin
488        SPAWN, 'mv -f ./'+denom+'????.ps '+plotfolder+'/'
489        SPAWN, 'cp -f graphwrf.def '+plotfolder+'/'+denom+'.def'
490end
491'r': begin
492        SPAWN, 'touch report.comments ; \rm report.comments ; echo '+textitle+' > report.comments ; echo '+texcomments+' >> report.comments'
493        SPAWN, 'rm -f '+denom+'_.ps ; echo executing ... base '+denom+'????.ps < report.comments'
494        SPAWN, 'base '+denom+'????.ps < report.comments ; mv -f base.ps '+denom+'_.ps'
495
496        SPAWN, 'mv -f ./'+denom+'????.ps '+plotfolder+'/'
497        SPAWN, 'mv -f ./'+denom+'_.ps '+plotfolder+'/'
498        SPAWN, 'cp -f graphwrf.def '+plotfolder+'/'+denom+'.def'
499end
500'm': begin
501        print, 'generating movie ...'
502        SPAWN, 'rm -rf temp ; mkdir temp ; mv -f '+denom+'????.ps temp ; cp bigconvert200 temp ; cd temp'
503        SPAWN, 'bigconvert200 *.ps ; convert -delay 60 *.png movie.gif ; cd ..'
504        if (nomovie ne 'yes') then SPAWN, 'animate temp/movie.gif &'
505        print, 'movie done'
506
507        SPAWN, 'mv -f temp/'+denom+'????.ps '+plotfolder+'/'
508        SPAWN, 'mv -f temp/movie.gif '+plotfolder+'/'+denom+'.gif'
509        SPAWN, 'mv -f temp/'+denom+'????.png '+plotfolder+'/'
510        SPAWN, 'cp -f graphwrf.def '+plotfolder+'/'+denom+'.def'
511end
512'w': begin
513;
514;
515netfolder='/u/aslmd/WWW/ANIMATIONS/'
516;
517; set also textitle
518;
519        print, 'generating frames ...'
520        SPAWN, 'rm -rf temp ; mkdir temp ; mv -f '+denom+'????.ps temp/ ; ls temp/'
521        SPAWN, "find temp/* -exec convert -density 100x100 -crop 0x0 {} '{}.png' \;"
522        SPAWN, 'ls temp/'
523        ;SPAWN, 'mv -f temp/'+denom+'????.ps '+plotfolder+'/'
524        SPAWN, 'cp -f graphwrf.def '+plotfolder+'/web_'+denom+'.def'
525        ;; webpage
526        ;SPAWN, 'rm -rf '+netfolder+textitle+'_'+denom+' ; mkdir '+netfolder+textitle+'_'+denom
527        SPAWN, 'mkdir '+netfolder+textitle
528        SPAWN, 'mv -f temp/'+denom+'????.ps.png '+netfolder+textitle+'/'
529        SPAWN, 'chmod 777 -R '+netfolder+textitle
530        SPAWN, 'rm -rf '+denom+'.html ; more header.html > '+textitle+'_'+denom+'.html'
531        genanim, num, textitle, denom
532        SPAWN, 'more templist >> '+textitle+'_'+denom+'.html ; rm -rf templist'
533        SPAWN, 'more body.html >> '+textitle+'_'+denom+'.html'
534        SPAWN, 'mv -f '+textitle+'_'+denom+'.html '+netfolder+'/'
535end     
536ENDCASE
537
538if (backup_data eq 'yes') then SPAWN, 'cp -f namelist.ARWpost '+datafolder+'/'+coord+'.namelist.ARWpost'
539if (backup_data eq 'yes') then SPAWN, 'cp -f input.ctl '+datafolder+'/'+coord+'.input.ctl'
540if (backup_data eq 'yes') then SPAWN, 'cp -f input.dat '+datafolder+'/'+coord+'.input.dat'
541
542
543
544
545end
Note: See TracBrowser for help on using the repository browser.