source: trunk/MESOSCALE/PLOT/SPEC/MAP/map_uvt.pro @ 205

Last change on this file since 205 was 160, checked in by aslmd, 14 years ago

MESOSCALE: modifications to use a tweaked polar topography corresponding to paleo conditions. plus purely cosmetic changes.
File size: 8.5 KB
Line 
1pro map_uvt
2;
3;
4;
5what_I_plot=0.
6overcontour=0
7no3d='true'
8@map_uvt_inc.pro
9SPAWN, '\rm param_plot.idl ; cp map_uvt_inc.pro param_plot.idl ; cp -f map_uvt_inc.pro '+save_ps+'.map_uvt_inc.pro'
10if (n_elements(coord2d) eq 0) then coord2d='false'
11if (n_elements(overvector_x) eq 0) then overvector_x = 1  ;; initialisation (!! map_uvt_inc.pro utilise dans map_uvt et map_latlon)
12;
13;
14;
15if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield1
16if (n_elements(field2) ne 0) then getcdf, file=filename, charvar=field2, invar=cfield2
17if ( ( overvector_x ne 0 ) or ( n_elements(field1) eq 0 ) ) then begin
18   u = getget(filename, 'Um', count=[0,0,1,0], offset=[0,0,nlevel,0])
19   v = getget(filename, 'Vm', count=[0,0,1,0], offset=[0,0,nlevel,0])
20   help, u
21endif
22if (n_elements(fieldscoord) eq 0) then begin
23 getcdf, file=filename, charvar='XLONG', invar=longi
24 getcdf, file=filename, charvar='XLAT', invar=lati
25 getcdf, file=filename, charvar='HGT', invar=hgt
26endif else begin
27 print, 'coordinates are ', fieldscoord
28 getcdf, file=filename, charvar=fieldscoord(0), invar=longi
29 getcdf, file=filename, charvar=fieldscoord(1), invar=lati
30 getcdf, file=filename, charvar=fieldscoord(2), invar=hgt
31endelse
32
33;
34;
35;
36sp = 5 ;; relaxation width
37nx = n_elements(longi(0,*))
38ny = n_elements(longi(*,0))
39if (n_elements(field1) ne 0) then begin
40        if (no3d ne 'true') then cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*,*) else cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*)
41endif
42if (n_elements(field2) ne 0) then begin
43        if (no3d ne 'true') then cfield2 = cfield2 (sp:nx-sp-1,sp:ny-sp-1,*,*) else cfield2 = cfield2 (sp:nx-sp-1,sp:ny-sp-1,*)
44endif
45if ( ( overvector_x ne 0 ) or ( n_elements(field1) eq 0 ) ) then begin
46  u             = reform (u             (sp:nx-sp-1,sp:ny-sp-1,0,*) )
47  v             = reform (v             (sp:nx-sp-1,sp:ny-sp-1,0,*) )
48endif
49longi           = longi         (sp:nx-sp-1,sp:ny-sp-1,*)
50lati            = lati          (sp:nx-sp-1,sp:ny-sp-1,*)
51hgt             = hgt           (sp:nx-sp-1,sp:ny-sp-1,*)
52nx = n_elements(longi(*,0,0))
53ny = n_elements(longi(0,*,0))
54nt = n_elements(longi(0,0,*))
55if ( overvector_x ne 0 ) then begin 
56  overvector_x = u
57  overvector_y = v
58endif
59;
60;
61;
62                if (ntime eq 99) then begin
63                        PRINT, '-- ALL TIME STEPS', nt-1
64                        ntstart = 0  & ntend = nt-1
65                endif else begin
66                        PRINT, '-- ONLY TIME STEP ', string(ntime,'(I0)')
67                        ntstart = ntime & ntend = ntime
68                endelse
69                for ntime = ntstart,ntend do begin
70help, u
71;
72;
73;
74  zefile=save_ps+string(100+ntime,'(I0)')
75  PS_Start, filename=zefile+'.ps'
76  print, zefile+'.ps'
77  !P.Charsize = 1.2
78  !p.charthick = 2.0
79  !p.thick = 2.0
80  !x.thick = 2.0
81  !y.thick = 2.0
82;
83;
84;
85lon = reform(longi(*,*,ntime))
86lat = reform(lati(*,*,ntime))
87if (n_elements(field1) eq 0) then begin
88        print, 'field1: horizontal velocity'
89                ;;;attention attention overvector_x et overvector_y sont vides
90                ;zevel = overvector_x^2 + overvector_y^2  ;; attention il faut que les tableaux soient de la meme taille, OK si uvmet
91        zevel = (reform(u(*,*,ntime)))^2 + (reform(v(*,*,ntime)))^2 ;; attention il faut que les tableaux soient de la meme taille, OK si uvmet
92        what_I_plot = sqrt(zevel)
93endif else begin
94        if (no3d ne 'true') then what_I_plot = reform(cfield1(*,*,nlevel,ntime)) else what_I_plot = reform(cfield1(*,*,ntime))
95        if (n_elements(u) ne 0) then overvector_x = reform(u(*,*,ntime))   ;; ne pas utiliser test overvector_x a cause de la boucle temps
96        if (n_elements(v) ne 0) then overvector_y = reform(v(*,*,ntime))   ;; ne pas utiliser test overvector_y a cause de la boucle temps
97endelse
98if (n_elements(field2) eq 0) then begin
99        print, 'field2: topography'
100        overcontour = reform(hgt(*,*,ntime))
101endif else begin
102        if (no3d ne 'true') then overcontour = reform(cfield2(*,*,nlevel,ntime)) else overcontour = reform(cfield2(*,*,ntime))
103endelse
104;
105;
106;
107if ( coord2d eq 'polar' ) then begin
108   print, 'OK YOU USE MAP_SET with POLAR PROJECTION. VECTORS ARE NOT SUPPORTED. USE polar_uv OR ADAPT THIS SCRIPT.'
109   overvector_x = 0
110   overvector_y = 0
111   if (n_elements(windowx) ne 0) then begin
112       latmin = windowy(0) & latmax = windowy(1) & lonmin = windowx(0) & lonmax = windowx(1)
113   endif else begin
114       latmin = 65. & latmax = 90. & lonmin = -180. & lonmax = 180.
115   endelse
116   print, 'latmin,lonmin,latmax,lonmax', latmin,lonmin,latmax,lonmax
117   map_set, latmax, 0., /isotropic, /azimuthal, /noborder, limit=[latmin,lonmin,latmax,lonmax], title=title_user, /advance
118endif else begin
119  minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon)
120  if (coord2d eq 'true') then begin
121        ;;;; CECI EST DESORMAIS FAIT DANS MAP_LATLON 
122        ;;npoints=n_elements(lon(*,0)) + n_elements(lon(0,*))  ;; trop de points, mais au moins on ne perd rien
123        ;;TRIANGULATE, lon, lat, tr
124        ;;what_I_plot  = GRIDDATA( lon, lat, what_I_plot,  /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
125        ;;if ( overvector_x ne 0 ) then begin
126        ;;  overvector_x = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
127        ;;  overvector_y = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
128        ;;endif
129        ;;overcontour  = GRIDDATA( lon, lat, overcontour,  /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
130        ;;                if ( overvector_x ne 0 ) then begin
131        ;;                ; sale sale sale
132        ;;                if (minlat lt min(lat(*,0))) then overvector_y=-overvector_y
133        ;;                if (minlon lt min(lon(0,*))) then overvector_x=-overvector_x
134        ;;                endif
135        ;;lon =  minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1)
136        ;;lat =  minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1)
137  endif else begin
138        lon = reform(lon(*,0))
139        lat = reform(lat(0,*))
140        ;;npoints=n_elements(lon(*,0)) + n_elements(lon(0,*))
141        ;;what_I_plot = REBIN( what_I_plot, npoints, npoints )
142        ;;overvector_x = REBIN( overvector_x, npoints, npoints )
143        ;;overvector_y = REBIN( overvector_y, npoints, npoints )
144        ;;overcontour = REBIN( overcontour, npoints, npoints )
145        ;;lon =  minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1)
146        ;;lat =  minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1)
147                ;lon = findgen( n_elements(what_I_plot(*,0)) )
148                ;lat = findgen( n_elements(what_I_plot(0,*)) )
149                ;print, 'ATTENTION!!!! NOUS SOMMES BIEN D ACCORD QUE VOUS NE CHARGEZ PAS UMET et VMET ??'
150  endelse
151                        ;;;;; trouve dans polar_uv.pro
152                        ;;if (coord2d eq 'regular') then begin ;; carte avec mercator
153                        ;;      lon = reform(lon(*,0))
154                        ;;      lat = reform(lat(0,*))
155                        ;;endif
156endelse
157help, what_I_plot, lon, lat
158print, 'COLOR PLOT min/max ', min(what_I_plot), max(what_I_plot)
159print, 'CONTOUR PLOT min/max ', min(overcontour), max(overcontour)
160print, 'VECTOR PLOT min/max ', min(overvector_x), min(overvector_y)
161;
162;
163;
164map_latlon, $
165        what_I_plot, $                          ; 2D field
166        lon, $                                  ; 1D latitude
167        lat, $                                  ; 1D longitude
168        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
169        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
170        overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
171        overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
172        overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
173        ct=pal, $                               ; color table (33-rainbow is default)
174        colors=colors, $                        ; number of colors/levels (32 is default)
175        title=title_user, $                     ; title of the plot ('' is default)
176        format=format                           ; format of colorbar annotations ('(F6.2)' is default)
177;
178;
179;
180if ( coord2d eq 'polar' ) then begin    ;;; pourrait aller dans map_latlon
181  loadct, 0
182  MAP_GRID, $
183          CHARSIZE = 1., $
184          COLOR    = 0,   $
185          ;lats=-60, $
186          LABEL    = 1,   $   ;; /LABEL or LABEL=2 (one label any 2 grid lines)
187          LATDEL   = intervaly, $   ;;5 10
188          LONDEL   = intervalx, $   ;;15
189          ;LONLAB   = latmin + intervaly/2., $ ;5. + (latmin+latmax)/2., $ ;0.
190          LONLAB   = (latmin+latmax)/2., $
191          LATLAB   = -0.001, $
192          GLINESTYLE = 2, $
193          GLINETHICK = 0.3
194          ;LONALIGN = 0., $
195          ;LATALIGN = 1.
196endif
197;
198;
199;
200PS_End, /PNG
201                endfor
202end
Note: See TracBrowser for help on using the repository browser.