source: trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro @ 91

pro map_uvt
;
;
;
what_I_plot=0.
overcontour=0
no3d='true'
@map_uvt_inc.pro
SPAWN, '\rm param_plot.idl ; cp map_uvt_inc.pro param_plot.idl ; cp -f map_uvt_inc.pro '+save_ps+'.map_uvt_inc.pro'
if (n_elements(coord2d) eq 0) then coord2d='false'
;
;
;
if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield1
if (n_elements(field2) ne 0) then getcdf, file=filename, charvar=field2, invar=cfield2
if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin
   u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,0])
   v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,0])
endif
getcdf, file=filename, charvar='XLONG', invar=longi
getcdf, file=filename, charvar='XLAT', invar=lati
getcdf, file=filename, charvar='HGT', invar=hgt
;
;
;
sp = 5 ;; relaxation width
nx = n_elements(longi(0,*))
ny = n_elements(longi(*,0))
if (n_elements(field1) ne 0) then begin
        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,*)
endif
if (n_elements(field2) ne 0) then begin
        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,*)
endif
if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin
  u             = u             (sp:nx-sp-1,sp:ny-sp-1,*) 
  v             = v             (sp:nx-sp-1,sp:ny-sp-1,*)
endif
longi           = longi         (sp:nx-sp-1,sp:ny-sp-1,*)
lati            = lati          (sp:nx-sp-1,sp:ny-sp-1,*)
hgt             = hgt           (sp:nx-sp-1,sp:ny-sp-1,*)
nx = n_elements(longi(*,0,0))
ny = n_elements(longi(0,*,0))
nt = n_elements(longi(0,0,*))
if ( n_elements( overvector_x ) ne 0 ) then begin  
  overvector_x = u
  overvector_y = v
endif
;
;
;
                if (ntime eq 99) then begin
                        PRINT, '-- ALL TIME STEPS', nt-1 
                        ntstart = 0  & ntend = nt-1 
                endif else begin
                        PRINT, '-- ONLY TIME STEP ', string(ntime,'(I0)')
                        ntstart = ntime & ntend = ntime
                endelse
                for ntime = ntstart,ntend do begin
;
;
;
  zefile=save_ps+string(100+ntime,'(I0)')
  PS_Start, filename=zefile+'.ps'
  print, zefile+'.ps'
  !P.Charsize = 1.2
  !p.charthick = 2.0
  !p.thick = 2.0
  !x.thick = 2.0
  !y.thick = 2.0
;
;
;
lon = reform(longi(*,*,ntime))
lat = reform(lati(*,*,ntime))
if (n_elements(field1) eq 0) then begin
        print, 'field1: horizontal velocity'
        zevel = overvector_x^2 + overvector_y^2  ;; attention il faut que les tableaux soient de la meme taille, OK si uvmet
        what_I_plot = sqrt(zevel)
endif else begin
        if (no3d ne 'true') then what_I_plot = reform(cfield1(*,*,nlevel,ntime)) else what_I_plot = reform(cfield1(*,*,ntime))
        if (n_elements(u) ne 0) then overvector_x = reform(overvector_x(*,*,ntime))   ;; ne pas utiliser test overvector_x a cause de la boucle temps
        if (n_elements(v) ne 0) then overvector_y = reform(overvector_y(*,*,ntime))   ;; ne pas utiliser test overvector_y a cause de la boucle temps
endelse
if (n_elements(field2) eq 0) then begin
        print, 'field2: topography'
        overcontour = reform(hgt(*,*,ntime))
endif else begin
        if (no3d ne 'true') then overcontour = reform(cfield2(*,*,nlevel,ntime)) else overcontour = reform(cfield2(*,*,ntime))
endelse
;
;
;
if ( coord2d eq 'polar' ) then begin
   print, 'OK YOU USE MAP_SET with POLAR PROJECTION. VECTORS ARE NOT SUPPORTED. USE polar_uv OR ADAPT THIS SCRIPT.'
   overvector_x = 0 
   overvector_y = 0
   if (n_elements(windowx) ne 0) then begin
       latmin = windowy(0) & latmax = windowy(1) & lonmin = windowx(0) & lonmax = windowx(1)
   endif else begin
       latmin = 65. & latmax = 90. & lonmin = -180. & lonmax = 180.
   endelse
   print, 'latmin,lonmin,latmax,lonmax', latmin,lonmin,latmax,lonmax
   map_set, latmax, 0., /isotropic, /azimuthal, /noborder, limit=[latmin,lonmin,latmax,lonmax], title=title_user, /advance
endif else begin
  minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon)
  if (coord2d eq 'true') then begin
        npoints=n_elements(lon(*,0)) + n_elements(lon(0,*))  ;; trop de points, mais au moins on ne perd rien 
        TRIANGULATE, lon, lat, tr 
        what_I_plot  = GRIDDATA( lon, lat, what_I_plot,  /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
        if (n_elements(overvector_x) ne 0) then begin
          overvector_x = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
          overvector_y = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
        endif
        overcontour  = GRIDDATA( lon, lat, overcontour,  /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
                        if (n_elements(overvector_x) ne 0) then begin
                          ; sale sale sale
                          if (minlat lt min(lat(*,0))) then overvector_y=-overvector_y
                          if (minlon lt min(lon(0,*))) then overvector_x=-overvector_x
                        endif
        lon =  minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1)
        lat =  minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1)
  endif else begin
        lon = reform(lon(*,0))
        lat = reform(lat(0,*))
        ;npoints=n_elements(lon(*,0)) + n_elements(lon(0,*))
        ;what_I_plot = REBIN( what_I_plot, npoints, npoints )
        ;overvector_x = REBIN( overvector_x, npoints, npoints )
        ;overvector_y = REBIN( overvector_y, npoints, npoints )
        ;overcontour = REBIN( overcontour, npoints, npoints )
        ;lon =  minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1)
        ;lat =  minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1)
  endelse
endelse
help, what_I_plot, lon, lat
print, min(what_I_plot)
print, max(what_I_plot)
print, min(overcontour)
print, max(overcontour)
;
;
;
map_latlon, $
        what_I_plot, $                          ; 2D field
        lon, $                                  ; 1D latitude
        lat, $                                  ; 1D longitude
        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
        overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
        overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
        overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
        ct=pal, $                               ; color table (33-rainbow is default)
        colors=colors, $                        ; number of colors/levels (32 is default)
        title=title_user, $                     ; title of the plot ('' is default)
        format=format                           ; format of colorbar annotations ('(F6.2)' is default)
;
;
;
if ( coord2d eq 'polar' ) then begin    ;;; pourrait aller dans map_latlon
  loadct, 0
  MAP_GRID, $
          CHARSIZE = 1., $
          COLOR    = 0,   $
          ;lats=-60, $
          LABEL    = 1,   $   ;; /LABEL or LABEL=2 (one label any 2 grid lines)
          LATDEL   = intervaly, $   ;;5 10
          LONDEL   = intervalx, $   ;;15
          LONLAB   = latmin + intervaly/2., $ ;5. + (latmin+latmax)/2., $ ;0.
          LATLAB   = -0.001, $
          GLINESTYLE = 2, $
          GLINETHICK = 0.3
          ;LONALIGN = 0., $
          ;LATALIGN = 1.
endif
;
;
;
PS_End, /PNG
                endfor
end