;------------------------------------------------------------------------------
;+
; NAME:
;       PCONTOUR
; PURPOSE:
;       Display or plot a color contoured image with axis and a color bar, 
;       using one of several projections (default is cylindrical equidistant).
; LAST UPDATE: 20 January 2000 
; CATEGORY:
; CALLING SEQUENCE:
;       pcontour    -> Displays help information for pcontour.
;       pcontour, a
; COMMON BLOCKS:
; MODIFICATION HISTORY:
;       W. Berg - April 3, 1995
;       D. Anderson, Feb 1999, Added _extra structure to pass add'l args
;       D. Anderson, Mar 1999, Added projection options
;       D. Anderson, June 1999, Added GIF options
;       D. Anderson, July 1999, Added region and smooth options
;       D. Anderson, Aug 1999, Device available colors check 
;       D. Anderson, Aug 1999, Change label defaults and imgsize handling
;       D. Anderson, Jan 2000, Correct reversal of 'cyl' + /order GIFs 
;       G. Wick, Apr 2000, Added /islands flag for hi-resolution map and coasts
;
; VARIABLE DESCRIPTIONS:
;   dcnt      Count of valid data values in passed array
;   dndex     Array index to valid data values in passed array
;   ls, ls2   Character size for labels.
;   px        Lower-left x coordinate for plotted image, in full page coords
;   py        Lower-left y coordinate for plotted image, in full page coords
;   pxc       Lower-left x coordinate for plotted colorbar, in full page coords
;   pyc       Lower-left y coordinate for plotted colorbar, in full page coords
;   scrnflag  0 if drawn to screen; 1 for PS, EPS, or GIF output 
;   ts, ts2   Character size for title.
;   xs        Width of plotted image, in full page coordiates
;   ys        Height of plotted image, in full page coordiates
;   xsc       Width of plotted colorbar, in full page coordiates
;   ysc       Height of plotted colorbar, in full page coordiates
;
;---------------------------------------------------------------------------
pro pcontour, a, title=title, subtitle=sttl, units=units, project=prj, $
              xlabel=xlab, ylabel=ylab, vmin=vmn, vmax=vmx, region=region, $ 
              minlat=minlat, maxlat=maxlat, minlon=minlon, maxlon=maxlon,  $
              nxminor=mxt, nyminor=myt, mis_val=misval, msk_val=mskval, $
              mis_col=miscol, msk_col=mskcol, cthick=cth, nofill=nofill, $
              lines=lines, multi=multi, nlevels=nlev, gifout=gifout, $
              ps_port=psp, ps_land=psl, eps_port=epsp, $
              eps_land=epsl, tick_off=tick_off, coast=coast, smooth=smooth, $
              ccolor=ccol, left_off=left_off, right=right, exp=exp, $
              top=top, bottom_off=bottom_off, cbar_off=cbar_off, $
              label=label, order=order, back_col=bcol, axis_col=acol, $
              lsize=ls, tsize=ts, position=pos, col_file=cfile, $
              ps_file=psfile, aspect=aspect, xtickname=xtn, ytickname=ytn, $
              xtickv=xtv, ytickv=ytv, vertical=vbar, horizontal=hbar, $
              ctickname=ctn, imgpos=imgpos, islands=islands, pngout=pngout, $
              outfile=outfile, _extra=xtra 

if (n_params(0) EQ 0) then begin
  print,' PCONTOUR  -> Displays help information for pcontour.'
  print,' PCONTOUR,a -> Displays a color image with axis and a color bar.'
  print,' Keywords:'
  print,'   ASPECT=aspect  Aspect ratio for the plot (def=(maxlon-minlon)/(maxlat-minlat).'
  print,'   SCALE=scale    Scale factor for the data.'
  print,'   OFFSET=offset  Offset factor for the data.'
  print,'   TITLE=title    Plot title'
  print,'   SUBTITLE=sttl  Secondary plot title'
  print,'   UNITS=units    Units label for the color bar'
  print,'   PROJECT=prj    Map projection (default: cylindrical equidistant)'
  print,'                  cyl  Cylindrical Equidistant'
  print,'                  mer  Mercator'
  print,'                  mol  Mollweide'
  print,'                  sat  Satellite'
  print,'                  polN Polar stereographic, north pole'
  print,'                  polS Polar stereographic, south pole'
  print,'                  xy   XY plot (no map)'
  print,'   XLABEL=xlab    X axis label.'
  print,'   YLABEL=ylab    Y axis label.'
  print,'   VMIN=vmn       Minimum data value.'
  print,'   VMAX=vmx       Maximum data value.'
  print,'   CMIN=cmn       Color value corresponding to the minimum data value.'
  print,'   CMAX=cmx       Color value corresponding to the maximum data value.'
  print,'   MINLAT=minlat  Minimum latitude value.'
  print,'   MAXLAT=maxlat  Maximum latitude value.'
  print,'   MINLON=minlon  Minimum longitude value.'
  print,'   MAXLON=maxlon  Maximum longitude value.'
  print,'   REGION=region  [minlat,maxlat,minlon,maxlon] array.'
  print,'   LEVEL=nl       Number of the level to display (1..n).'
  print,'   TIME=nt        Number of the time band to display (1..n).'
  print,'   VAR=nv         Number of the variable to display (1..n).'
  print,'   NLEVELS=nlev   Number of contour levels.'
  print,'   NXTICK=nxt     Number of major tick mark intervals along x'
  print,'   NYTICK=nyt     Number of major tick mark intervals along y.'
  print,'   NCTICK=nct     Number of tick mark intervals for the color bar.'
  print,'   NXMINOR=mxt    Minor tick intervals per major tick along x.'
  print,'   NYMINOR=myt    Minor tick intervals per major tick along y.'
  print,'   XTICKNAME=xtn  Array of tick names for the x axis.'
  print,'   YTICKNAME=ytn  Array of tick names for the y axis.'
  print,'   CTICKNAME=ctn  Array of tick names for the colorbar axis.'
  print,'   TSIZE=ts       Character size for title.'
  print,'   LSIZE=ls       Character size for labels.'
  print,'   BACK_COL=bcol  Color for the background.'
  print,'   AXIS_COL=acol  Color for axis lines and labels.'
  print,'   MIS_VAL=misval The missing data value. (def=-999)'
  print,'   MSK_VAL=mskval The data value to mask out. (def=-998)'
  print,'   MIS_COL=miscol The color value for missing data. (def=0)'
  print,'   MSK_COL=mskcol The color value for mask data. (def=2)'
  print,'   CCOLOR=ccol    Color for coastal outlines.'
  print,'   GCOLOR=gcol    Color for grid lines.'
  print,'   CTHICK=cth     Thickness of coastal outlines.'
  print,'   GTHICK=gth     Thickness of grid lines.'
  print,'   ATHICK=ath     Thickness of axis lines.'
  print,'   CTYPE=ctyp     Line style for coastal outlines.'
  print,'   GTYPE=gtyp     Line style for grid lines.'
  print,'   POSITION=pos   Position of plot within window'
  print,'                  (def=[0,0,1,1]->whole window).'
  print,'   IMGPOS=imgpos  Position of the image within the plot window.'
  print,'                  (image x1,y1,wid,ht & colorbar x1,y1,wid,ht)'
  print,'   COL_FILE=cfile Name of the ascii color table file.'
  print,'   PS_FILE=psfile Name of the output PostScript file.'
  print,'   OUTFILE=outfileName of the output GIF or PNG file.'
  print,'   SMOOTH=bcar    Boxcar smooth data before contour (i.e. 1,3,5 etc.)'
  print,'   /NOFILL        Turn off the contour fill.'
  print,'   /LINES         Turn on the contour lines and labels.'
  print,'   /MULTI         Set flag for doing multiple plots (do not set for last plot).'
  print,'   /PS_PORT       Write output to PostScript file in portrait mode.'
  print,'   /EPS_PORT      Write output to Encapsulated PS file, portrait mode.'
  print,'   /PS_LAND       Write output to PostScript file in landscape mode.'
  print,'   /EPS_LAND      Write output to Encapsulated PS file, landscape mode.'
  print,'   /GIFOUT        Write output as GIF image.'
  print,'   /COAST         Draw coastlines over the image.'
  print,'   /ISLANDS       Use high-resolution map and draw islands.'
  print,'   /GRID          Draw a grid over the map.'
  print,'   /BOTTOM_OFF    Turn off axis on the bottom of plot.'
  print,'   /LEFT_OFF      Turn off axis on the left side of plot.'
  print,'   /TOP           Add and label axis on top of plot.'
  print,'   /RIGHT         Add and label axis on right side of plot.'
  print,'   /TICK_OFF      Turn off the tick marks on the axis.'
  print,'   /CBAR_OFF      Turn off the color bar.'
  print,'   /VERTICAL      Force the color bar to be vertically oriented.'
  print,'   /HORIZONTAL    Force the color bar to be horizontally oriented.'
  print,'   /LABEL         Label axis E/W and N/S (default is  +/-)'
  print,'   /ORDER         Reverse order of image (top to bottom)'
  print,'   /VERBOSE       Print status information'
  return
endif

; Default variable values

 cmn=5
 cmx=255
 nl=1
 nt=1
 nv=1
 nxt=6
 nyt=6
 offset=0.0
 sublat=0.0
 sublng=-135.0
 scale=1.0
 scrnflag=0
 xws=800
 yws=600

; Check for above variables passed as "extra" (non-explicit) arguments

 varcnt = n_tags(xtra)
 if varcnt gt 0 then begin
   varstr=tag_names(xtra)
   for i = 1, varcnt do begin
     if strpos(varstr(i-1), 'LEVEL') gt -1 then nl=xtra.level
     if strpos(varstr(i-1), 'TIME') gt -1 then nt=xtra.time
     if strpos(varstr(i-1), 'VAR') gt -1 then nv=xtra.var
     if strpos(varstr(i-1), 'NCTICK') gt -1 then nct=xtra.nctick
     if strpos(varstr(i-1), 'NXTICK') gt -1 then nxt=xtra.nxtick
     if strpos(varstr(i-1), 'NYTICK') gt -1 then nyt=xtra.nytick
     if strpos(varstr(i-1), 'SUBLAT') gt -1 then sublat=xtra.sublat
     if strpos(varstr(i-1), 'SUBLNG') gt -1 then sublng=xtra.sublng
     if strpos(varstr(i-1), 'SCALE') gt -1 then scale=xtra.scale
     if strpos(varstr(i-1), 'OFFSET') gt -1 then offset=xtra.offset
     if strpos(varstr(i-1), 'CMIN') gt -1 then cmn=xtra.cmin
     if strpos(varstr(i-1), 'CMAX') gt -1 then cmx=xtra.cmax
     if strpos(varstr(i-1), 'XWS') gt -1 then xws=xtra.xws
     if strpos(varstr(i-1), 'YWS') gt -1 then yws=xtra.yws
     if strpos(varstr(i-1), 'NPREC') gt -1 then nprec=xtra.nprec
     if strpos(varstr(i-1), 'USA') gt -1 then usa=xtra.usa
     if strpos(varstr(i-1), 'XSIZE') gt -1 then xsz=xtra.xsize
     if strpos(varstr(i-1), 'YSIZE') gt -1 then ysz=xtra.ysize
     if strpos(varstr(i-1), 'LCOLOR') gt -1 then lcol=xtra.lcolor
     if strpos(varstr(i-1), 'GTHICK') gt -1 then gth=xtra.gthick
     if strpos(varstr(i-1), 'ATHICK') gt -1 then ath=xtra.athick
     if strpos(varstr(i-1), 'VERBOSE') gt -1 then verbose=xtra.verbose
   end
 endif

; Read in default header values from common block

if n_elements(nl) eq 0 then nl=1
if n_elements(nt) eq 0 then nt=1
if n_elements(nv) eq 0 then nv=1
help, a  
a2=hdr_intf(a,title,sttl,units,minlat,maxlat,minlon,maxlon,misval,mskval,$
            order,nl,nt,nv,verbose=verbose)

;-------  Set defaults  ---------------

if n_elements(scale) eq 0 then scale=1.0
if n_elements(offset) eq 0 then offset=0.0
if n_elements(title) eq 0 then title=''
if n_elements(sttl) eq 0 then sttl=''
if n_elements(prj) eq 0 then prj='cyl'
if n_elements(units) eq 0 then units=''
if n_elements(xlab) eq 0 then xlab=''
if n_elements(ylab) eq 0 then ylab=''
if n_elements(region) eq 0 then begin
  if n_elements(minlat) eq 0 then minlat=-90
  if n_elements(maxlat) eq 0 then maxlat=90
  if n_elements(minlon) eq 0 then minlon=0
  if n_elements(maxlon) eq 0 then maxlon=360
endif
if n_elements(nprec) eq 0 then nprec=0
if n_elements(cmn) eq 0 then cmn=5
if n_elements(cmx) eq 0 then cmx=255
if n_elements(nxt) eq 0 then nxt=6
if n_elements(nyt) eq 0 then nyt=6
if n_elements(mxt) eq 0 then mxt=0
if n_elements(myt) eq 0 then myt=0
if n_elements(misval) eq 0 then misval=-999
if n_elements(mskval) eq 0 then mskval=-998
if n_elements(miscol) eq 0 then miscol=0
if n_elements(mskcol) eq 0 then mskcol=2
if n_elements(cth) eq 0 then cth=1
if n_elements(gth) eq 0 then gth=1
if n_elements(ath) eq 0 then ath=1
if n_elements(ctyp) eq 0 then ctyp=0
if n_elements(gtyp) eq 0 then gtyp=1
if n_elements(ccol) eq 0 then ccol=0
if n_elements(gcol) eq 0 then gcol=0
if n_elements(bcol) eq 0 then bcol=1
if n_elements(acol) eq 0 then acol=0
if n_elements(lcol) eq 0 then lcol=0
if n_elements(ts) eq 0 then ts=1.0
if n_elements(ls) eq 0 then ls=1.0
if n_elements(pos) lt 4 then pos=[0.0,0.0,1.0,1.0]
if n_elements(cfile) eq 0 then cfile=''
if n_elements(outfile) eq 0 then begin
  if keyword_set(gifout) then begin
    outfile='idl.gif'
  endif else begin
    outfile='idl.png'
  endelse
endif
if n_elements(psfile) eq 0 then begin
  if keyword_set(psp) or keyword_set(psl) then psfile='idl.ps'
  if keyword_set(epsp) or keyword_set(epsl) then psfile='idl.epsi'
endif
if n_elements(xsz) eq 0 then begin
  if keyword_set(psp) or keyword_set(epsp) then xsz=6.5
  if keyword_set(psl) or keyword_set(epsl) then xsz=9.0
endif
if n_elements(ysz) eq 0 then begin
  if keyword_set(psp) or keyword_set(epsp) then ysz=9.0
  if keyword_set(psl) or keyword_set(epsl) then ysz=6.5
endif
if n_elements(nlev) eq 0 then nlev=10
if n_elements(nct) eq 0 then nct=nlev
if keyword_set(nofill) then lines=1
if prj ne 'cyl' AND prj ne 'mer' AND prj ne 'xy' then begin
  left_off=1                      ; No need for rectangular x-y axis 
  bottom_off=1                    ; No need for rectangular x-y axis 
  endif
if n_elements(region) eq 4 then begin
  minlat = region(0)
  maxlat = region(1)
  minlon = region(2)
  maxlon = region(3)
  if (minlat gt maxlat) OR (minlon gt maxlon) then begin
    print, 'Error: problem with specified values for region!'
    stop
  endif
endif
if n_elements(aspect) eq 0 then begin
  case prj of 
    'xy':  aspect=float(maxlon-minlon)/float(maxlat-minlat)
    'cyl': aspect=float(maxlon-minlon)/float(maxlat-minlat)
    'mer': aspect=float(maxlon-minlon)/float(maxlat-minlat)
    'mol': aspect=float(maxlon-minlon)/float(maxlat-minlat)
    'sat': aspect=1.0
    'polN': aspect=1.0
    'polS': aspect=1.0
  endcase
endif
if n_elements(satpos) eq 0 then satpos=[6.61,0,0] ;GOES elev in globe radii

; Open the window or PostScript file or Z buffer for plotting

 ts2=ts
 ls2=ls
 if keyword_set(psp) OR keyword_set(psl) OR $
    keyword_set(epsp) OR keyword_set (epsl) then begin
   if (!p.font LE 0) then !p.font=0
   set_plot, 'ps'
   load_ctbl,cfile,ncolors,cmin=cmn,cmax=cmx
   if keyword_set(psp) then begin
     if (!d.unit EQ 0) then device,filename=psfile,/COLOR,BITS_PER_PIXEL=8, $
         /PORTRAIT,xsize=xsz,ysize=ysz,xoffset=1.0,yoffset=1.0,/inches, $
         /Helvetica,font_index=6
   endif 
   if keyword_set(psl) then begin
     if (!d.unit EQ 0) then device,filename=psfile,/COLOR,BITS_PER_PIXEL=8, $
         /LANDSCAPE,xsize=xsz,ysize=ysz,xoffset=1.0,yoffset=10.0,/inches, $
         /Helvetica,font_index=6
   endif
 
   if keyword_set(epsp) then begin
     if (!d.unit EQ 0) then device, filename=psfile, /COLOR, BITS_PER_PIXEL=8, $
        /PORTRAIT, xsize=xsz, ysize=ysz, xoffset=0.0, yoffset=0.0, /inches, $
        /encapsul, /Helvetica,font_index=6
   endif
   if keyword_set(epsl) then begin
     if (!d.unit EQ 0) then device, filename=psfile, /COLOR, BITS_PER_PIXEL=8, $
       /LANDSCAPE, xsize=xsz, ysize=ysz, xoffset=0.0, yoffset=9.0, /inches, $
       /encapsul, /Helvetica,font_index=6
   endif
   if (!p.multi(0) EQ 0) then begin
     tvlct,r,g,b,/get
     if (r(bcol) NE 255 OR g(bcol) NE 255 OR b(bcol) NE 255) then begin
       back=bytarr(10,10)+bcol
       TV,back,0.0,0.0,XSIZE=1.0,YSIZE=1.05,/norm
     endif
   endif
   xws=!d.x_size
   yws=!d.y_size
   scrnflag=1
 endif 

 if keyword_set(gifout) OR keyword_set(pngout) then begin
   if (!p.multi(0) LE 0) then begin
     set_plot, 'z'
     device, set_resolution=[xws,yws], z_buffering=0
   endif else begin
     xws=!d.x_size
     yws=!d.y_size
   endelse
   load_ctbl, cfile, ncolors, cmin=cmn, cmax=cmx
;  if !d.n_colors gt 256 then device, decomposed=0, /install_color
   if !d.n_colors gt 256 then device, decomposed=0
   if (!p.multi(0) EQ 0) then erase, bcol
   scrnflag=1
   ts2=ts2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
   ls2=ls2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
 endif

 if ( scrnflag EQ 0 ) then begin
   set_plot,'x'
   if (!d.window LT 0) then begin
     window,xsize=xws,ysize=yws
   endif else begin
     xws=!d.x_size
     yws=!d.y_size
   endelse
   load_ctbl,cfile,ncolors,cmin=cmn,cmax=cmx
   if !d.n_colors gt 256 then device, decomposed=0, /install_color
   if (!p.multi(0) EQ 0) then erase,bcol
   ts2=ts2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
   ls2=ls2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
 endif

 if keyword_set(multi) then begin
   !p.multi(0)=!p.multi(0)+1
 endif else begin
   !p.multi(0)=0
 endelse
 pmulti=!p.multi(0)
 if keyword_set(order) then !order=1 else !order=0

; Apply scale and offset to input data; index the missing and masked pixels

 n=size(a2)
 a3=fltarr(n(1),n(2))
 a3=a2*scale-offset
 dndex=where(a2 NE misval AND a2 NE mskval,dcnt)
 mndex=where(a2 EQ misval OR a2 EQ mskval,mcnt)
 if (dcnt GT 0) then begin
   if n_elements(vmn) eq 0 then vmn=min(a3(dndex))
   if n_elements(vmx) eq 0 then vmx=max(a3(dndex))
 endif else begin
   if n_elements(vmn) eq 0 then vmn=min(a3)
   if n_elements(vmx) eq 0 then vmx=max(a3)
 endelse

; Smooth the input data if so requested

 if ( keyword_set(smooth) ) then begin
   if ( smooth GT 1 ) then begin
     if mcnt GT 0 then a3(mndex) = !VALUES.F_NAN
     a3s = smooth ( a3, smooth, /NAN, /edge_truncate )
     if mcnt GT 0 then a3s(mndex) = misval
     a3 = a3s
   endif
 endif

; Scale the input file to byte values

 x=fltarr(n(1))
 y=fltarr(n(2))
 x(0)=minlon
 x(n(1)-1)=maxlon
 y(0)=minlat
 y(n(2)-1)=maxlat
 for i=1,n(1)-2 do $
   x(i)=float(minlon)+(float(i)+0.5)*float(maxlon-minlon)/float(n(1))
 for j=1,n(2)-2 do $
   y(j)=float(minlat)+(float(j)+0.5)*float(maxlat-minlat)/float(n(2))
 index=where(a3 GT vmx,cnt)
 if (cnt GT 0) then a3(index)=vmx
 index=where(a3 LT vmn,cnt)
 if (cnt GT 0) then a3(index)=vmn
 if (mcnt GT 0) then a3(mndex)=vmx+2

; Calculate the position and size of the output image

 imgsize, pos, xws, yws, xs, ys, px, py, xsc, ysc, pxc, pyc, tpos1, tpos2, $
          orient, vmn, vmx, minlon, maxlon, minlat, maxlat, aspect, ts2, ls2, $
          xtl, ytl, ctl, xlab, ylab, title, sttl, units, cbar_off, $
          tick_off, left_off, bottom_off, right, top, ytn, vbar, hbar , $
          verbose=verbose

 if ( n_elements(imgpos) EQ 0 ) then begin
   xscl=pos[2]-pos[0]
   yscl=pos[3]-pos[1]
   imgpos=fltarr(8)
   imgpos[0]=(px-pos[0])/xscl
   imgpos[1]=(py-pos[1])/yscl
   imgpos[2]=xs/xscl
   imgpos[3]=ys/yscl
   imgpos[4]=(pxc-pos[0])/xscl
   imgpos[5]=(pyc-pos[1])/yscl
   imgpos[6]=xsc/xscl
   imgpos[7]=ysc/yscl
   if keyword_set(verbose) then print,'imgpos = ',imgpos
 endif else begin
   tpos1=tpos1-py-ys
   tpos2=tpos2-py-ys
   xscl=pos[2]-pos[0]
   yscl=pos[3]-pos[1]
   px=imgpos[0]*xscl+pos[0]
   py=imgpos[1]*yscl+pos[1]
   xs=imgpos[2]*xscl
   ys=imgpos[3]*yscl
   pxc=imgpos[4]*xscl+pos[0]
   pyc=imgpos[5]*yscl+pos[1]
   xsc=imgpos[6]*xscl
   ysc=imgpos[7]*yscl
   tpos1=tpos1+py+ys
   tpos2=tpos2+py+ys
 endelse

; Set up map coordinates

 !P.COLOR=acol

 case prj of
  'cyl': begin
           map_set, 0.0, (minlon+maxlon)/2.0, 0, $
                   /cylindrical, limit=[minlat,minlon,maxlat,maxlon], $
                   /noerase, xmargin=[0,0], ymargin=[0,0], $
                   pos=[px,py,px+xs,py+ys], /noborder
           ; For some reason gifouts are reverse-mirrored w/o the following!
           if keyword_set(order) and (keyword_set(gifout) or $
              keyword_set(pngout)) then begin
             !order=0
           end
         end
  'mer': begin
           map_set, 0.0, (minlon+maxlon)/2.0, 0, $
           /mercator, limit=[minlat,minlon,maxlat,maxlon], $
           /noerase, xmargin=[0,0], ymargin=[0,0], $
           pos=[px,py,px+xs,py+ys], /noborder
           ; Projection requires N-S ordering of rows
           if keyword_set(order) then begin
             a3new=reverse(a3,2)
             a3=a3new
             !order=0
           end
         end
  'mol': begin
           map_set, 0.0, (minlon+maxlon)/2.0, 0, $
           /mollweide, limit=[minlat,minlon,maxlat,maxlon], $
           /noerase, xmargin=[0,0], ymargin=[0,0], $
           pos=[px,py,px+xs,py+ys], /noborder, /horizon
           ; Projection requires N-S ordering of rows
           if keyword_set(order) then begin
             a3new=reverse(a3,2)
             a3=a3new
             !order=0
           end
         end
  'sat': begin
           map_set, sublat, sublng, /satellite, $
           limit=[minlat,minlon,maxlat,maxlon], $
           /noerase, xmargin=[0,0], ymargin=[0,0], $
           pos=[px,py,px+xs,py+ys], sat_p=satpos, /noborder
           ; Projection requires N-S ordering of rows
           if keyword_set(order) then begin
             a3new=reverse(a3,2)
             a3=a3new
             !order=0
           end
         end
  'polN': begin
            map_set, 90.0, (minlon+maxlon)/2.0, 0, $
            /stereographic, limit=[0.0, minlon, 90.0, maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], $
            pos=[px,py,px+xs,py+ys], /noborder
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              a3new=reverse(a3,2)
              a3=a3new
              !order=0
            end
          end
  'polS': begin
            map_set, -90.0, (minlon+maxlon)/2.0, 0, $
            /stereographic, limit=[-90.0, minlon, 0.0, maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], $
            pos=[px,py,px+xs,py+ys], /noborder
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              a3new=reverse(a3,2)
              a3=a3new
              !order=0
            end
          end
  'xy':  begin
           plot,[minlon,maxlon],[minlat,maxlat],position=[px,py,px+xs,py+ys], $
                xstyle=5,ystyle=5,/nodata,/noerase,/normal
         end
   else: begin
           print, 'Requested projection not recognized'
           stop
        end
 endcase

 xyouts, px+xs/2.0, tpos1, string(title,format='("!3",a)'), charsize=ts2*1.5, $
         /normal, alignment=0.5
 xyouts, px+xs/2.0, tpos2, string(sttl,format='("!3",a)'), charsize=ts2*1.0, $
         /normal,alignment=0.5
 !P.CHARSIZE=ls2

; Resize and plot contours

 concol=fltarr(nlev+1)
 concl2=fltarr(nlev+1)
 conlev=fltarr(nlev+1)
 conlab=intarr(nlev+1)
 consty=intarr(nlev+1)
 for i=0,nlev do begin
   if (keyword_set(nofill)) then begin
     concl2(i)=cmn+i*(ncolors-cmn-1)/float(nlev-1)
     if (i MOD 2 EQ 0) then conlab(i)=1 else conlab(i)=0
   endif else begin
     concl2(i)=lcol
     conlab(i)=0
   endelse
   concol(i)=cmn+i*(ncolors-cmn-1)/float(nlev-1)
   conlev(i)=vmn+i*float(vmx-vmn)/float(nlev)
   if (conlev(i) LT 0.0) then consty[i]=1 else consty[i]=0
 endfor
 if keyword_set(exp) then begin
   scl=vmx/2.0^nct
   conlev[0]=1.0*scl
   for i=1,nlev do begin
     xt=float(i*nct)/float(nlev)
     conlev[i]=2.0^xt*scl
   endfor
 endif
 if (keyword_set(nofill)) then begin
   concl2(nlev)=ncolors-1
 endif else begin
   concl2(nlev)=lcol
 endelse
 concol(nlev)=ncolors-1
 conlbl=fltarr(nct+1)
 if keyword_set(exp) then begin
   scl=vmx/float(2^nct)
   ctn=strarr(nct+1)
   for i=0,nct do begin
     xt=2^i
     xt=xt*scl
     ctn[i] = fixfmt(xt)
   endfor
 endif
 for i=0,nct do begin
   conlbl[i]=float(vmn)+float(i)*float(vmx-vmn)/float(nct)
 endfor

 if (dcnt GT 0) then begin
   if (mcnt GT 0 AND NOT keyword_set(nofill)) then begin
     polyfill,[px,px+xs,px+xs,px,px],[py,py,py+ys,py+ys,py],color=miscol,/norm
   endif
   if keyword_set(order) AND prj eq 'cyl' then begin
     if (NOT keyword_set(nofill)) then CONTOUR,rotate(a3,7),x,y,/overplot, $
       position=[px,py,px+xs,py+ys],/cell_fill,nlevels=nlev,c_colors=concol, $
       levels=conlev,max_value=vmx+1
     if (keyword_set(lines)) then CONTOUR,rotate(a3,7),x,y,/overplot, $
       position=[px,py,px+xs,py+ys],/follow,nlevels=nlev,c_colors=concl2, $
       levels=conlev,max_value=vmx+1,c_labels=conlab,c_linestyle=consty
   endif else begin
     if (NOT keyword_set(nofill)) then CONTOUR,a3,x,y,/overplot, $
       position=[px,py,px+xs,py+ys],/cell_fill,nlevels=nlev,c_colors=concol, $
       levels=conlev,max_value=vmx+1
     if (keyword_set(lines)) then CONTOUR,a3,x,y,/overplot, $
       position=[px,py,px+xs,py+ys],nlevels=nlev,c_colors=concl2, $
       levels=conlev,max_value=vmx+1,c_labels=conlab,c_linestyle=consty
   endelse
 endif else begin
     XYOUTS,(minlon+maxlon)/2.0,(minlat+maxlat)/2.0,'Data unavailable',$
            ALIGNMENT=0.5,SIZE=1.4
 endelse

; Draw the axis

 case prj of
  'xy':  xt_axis,nxt,nyt,minlon,maxlon,minlat,maxlat,xtl,ytl,xtn,ytn,xtv,ytv, $
                 xlab,ylab,mxt,myt,left_off,bottom_off,label,right,top, $
                 tlabels,order,thick=ath
  'cyl': map_axis, nxt, nyt, minlon, maxlon, minlat, maxlat, xtl, $
                   ytl, xtn, ytn, xtv, ytv, xlab, ylab, mxt, myt, left_off, $
                   bottom_off, label, right, top, nprec=nprec, thick=ath
  'mer': if keyword_set(label) then begin 
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
         endif else begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
         endelse
  'mol': if keyword_set(label) then begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
         endif else begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
         endelse
  'sat': if keyword_set(label) then begin 
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
         endif else begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
         endelse
  'polN': if keyword_set(label) then begin 
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
          endif else begin
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
          endelse
  'polS': if keyword_set(label) then begin 
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
          endif else begin
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
          endelse
 endcase


; Draw the continental outlines and map grid

 if (keyword_set(coast) AND (prj NE 'xy')) then begin
   if keyword_set(usa) then begin
     if keyword_set(islands) then begin
       map_continents, mlinethick=cth, mlinestyle=ctyp, /coasts, /hires, $
                       /usa, color=ccol
     endif else begin
       map_continents, mlinethick=cth, mlinestyle=ctyp, /usa, /continents, $
                       color=ccol
     endelse
   endif else begin
     if keyword_set(islands) then begin
       map_continents, mlinethick=cth, mlinestyle=ctyp, /coasts, /hires, $
                       color=ccol
     endif else begin
       if (coast EQ 2) then begin
         map_continents, mlinethick=cth, color=ccol, mlinestyle=ctyp, $
                         /fill_continents
       endif else begin
         map_continents, mlinethick=cth, color=ccol, mlinestyle=ctyp
       endelse
     endelse
   endelse
 endif


; Draw outline around plot

 case prj of
   'xy':  plot,[minlon,maxlon],[minlat,maxlat],/nodata,/xstyl,/ystyl,xticks=1, $
          xtickn=[' ',' '],yticks=1,ytickn=[' ',' '], $
          position=[px,py,px+xs,py+ys],/noerase,xthick=ath,ythick=ath
   'cyl': plot,[minlon,maxlon],[minlat,maxlat],/nodata,/xstyl,/ystyl,xticks=1, $
          xtickn=[' ',' '],yticks=1,ytickn=[' ',' '], $
          position=[px,py,px+xs,py+ys],/noerase,xthick=ath,ythick=ath
   'mer': plot,[minlon,maxlon],[minlat,maxlat],/nodata,/xstyl,/ystyl,xticks=1, $
          xtickn=[' ',' '],yticks=1,ytickn=[' ',' '], $
          position=[px,py,px+xs,py+ys],/noerase,xthick=ath,ythick=ath
   'mol': map_grid, color=acol,glinestyle=0,glinethick=gth, $
          lats=[minlat,maxlat], $
          lons=[minlon,maxlon],xthick=ath,ythick=ath
   'sat': map_horizon, thick=ath, noclip=1
   'polN': map_horizon, thick=ath, noclip=1
   'polS': map_horizon, thick=ath, noclip=1
    else: print,'Plot outline not drawn!'
 endcase

 !p.multi(0)=pmulti

; Draw the colorbar

 if (cmx GT ncolors) then cmx=cmx-(256-ncolors)
 if (NOT keyword_set(cbar_off)) then begin
   if (orient EQ 1) then begin
     if (n_elements(ctn) EQ 0) then begin
       cont_cbar,pxc,pxc+xsc,pyc,pyc+ysc,conlbl,concol,/norm,unit=units, $
         charsize=ls2,ticklen=ctl,thick=ath
     endif else begin
       cont_cbar,pxc,pxc+xsc,pyc,pyc+ysc,conlbl,concol,/norm,unit=units, $
         ctickname=ctn,charsize=ls2,ticklen=ctl,thick=ath
     endelse
   endif else begin
     if (n_elements(ctn) EQ 0) then begin
       cont_cbar,pxc,pxc+xsc,pyc,pyc+ysc,conlbl,concol,/norm,unit=units, $
         charsize=ls2,ticklen=ctl,thick=ath
     endif else begin
       cont_cbar,pxc,pxc+xsc,pyc,pyc+ysc,conlbl,concol,/norm,unit=units, $
         ctickname=ctn,charsize=ls2,ticklen=ctl,thick=ath
     endelse
   endelse
 endif

; Close the output PostScript file

 if (keyword_set(psp) OR keyword_set(psl)) AND (!p.multi(0) LE 0) then begin
   device,/close_file
   set_plot,'x'
 endif

 if keyword_set(gifout) AND !p.multi(0) LE 0 then begin
  image = tvrd()
  tvlct, /get, r, g, b
  print, outfile
  write_gif, outfile, image, r, g, b
 endif

 if keyword_set(pngout) AND !p.multi(0) LE 0 then begin
  image = tvrd()
  tvlct, /get, r, g, b
  print, outfile
  write_png, outfile, image, r, g, b
 endif

; Set the POSVEC system variable, which makes the page position of the last
; image drawn available to the next called routine (e.g., overplot.pro)

 pos=[px,py,px+xs,py+ys]
 defsysv, '!POSVEC', exists=i
 if i eq 1 then begin !POSVEC = pos
 endif else defsysv, '!POSVEC', pos

; Set the PROJ system variable, which makes the projection of the last image
; drawn available to the next called routine (e.g., overplot.pro)

 defsysv, '!PROJ', exists=i
 if i eq 1 then begin !PROJ = prj
 endif else defsysv, '!PROJ', prj

; Set the SUBLAT and SUBLNG system variables

 defsysv, '!SUBLAT', exists=i
 if i eq 1 then begin !SUBLAT = sublat
 endif else defsysv, '!SUBLAT', sublat

 defsysv, '!SUBLNG', exists=i
 if i eq 1 then begin !SUBLNG = sublng
 endif else defsysv, '!SUBLNG', sublng

 defsysv, '!CMIN', exists=i
 if i eq 1 then begin !CMIN = cmn
 endif else defsysv, '!CMIN', cmn

 defsysv, '!CMAX', exists=i
 if i eq 1 then begin !CMAX = cmx
 endif else defsysv, '!CMAX', cmx

; Record the bounding lat/lons in the LIMIT system variable

 defsysv, '!LIMIT', exists=i
 if i eq 1 then begin !LIMIT = [minlat,minlon,maxlat,maxlon]
 endif else defsysv, '!LIMIT', [minlat,minlon,maxlat,maxlon]

end