pro pmap_lev_mls_100hpa,data,longrid,latgrid,dlon,dlat,levels=levels,title=title, $
 region=region,figname=figname,colorindex=colorindex

; Make a postscript file of a map variable
; Input:  
;  data = an array of data (longrid,latgrid) elements
;  latgrid = one dimensional array of latitudes
;  longrid = one dimensional array of longitudes
;  levels= the contour levels to be plotted (maxium 26 elements is appropriate)
;  labels= labels on the color scale
;  title = title on figure
;  region = plot region, e.g. glb, asa, ind, usa, pac, nor
;  figname = postscript figure name (without .ps)
;
; Output:
;  figname.ps
;-

if n_elements(levels) eq 0 and n_elements(labels) eq 0 then $
 levels=findgen(12)*(max(data)-min(data))/12+min(data)
if n_elements(title) eq 0 then title = figname
if n_elements(region) eq 0 then region = 'glb'

if region eq 'glb' then begin
 xr1=-180 & xr2=180
 yr1= -90 & yr2= 90
 xx=indgen(7)*60 -180
 yy=indgen(7)*30 - 90
endif

if region eq 'ind' then begin
 xr1= 40 & xr2=100
 yr1=  0 & yr2= 40
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'asa' then begin
 xr1= 60 & xr2=160
 yr1= -10 & yr2= 50
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'sea4rc' then begin
 xr1= 60 & xr2=160
 yr1= -10 & yr2= 50
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'sas' then begin
 xr1= 100 & xr2=180
 yr1= 0 & yr2= 60
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'usa' then begin
 xr1=-130 & xr2=-60
 yr1=  24 & yr2= 54
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'usa1' then begin
 xr1=-170 & xr2=-40
 yr1=  10 & yr2= 80
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'cal' then begin
 xr1=-125 & xr2=-110
 yr1=  25 & yr2= 45
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'eur' then begin
 xr1= -20 & xr2= 60
 yr1=  35 & yr2= 75
 xx=indgen( (xr2-xr1)/10 +1 )*10 +xr1
 yy=indgen( (yr2-yr1)/10 +1 )*10 +yr1
endif

if region eq 'nor' then begin
 xr1=-180 & xr2=180 
 yr1=   0 & yr2= 90
 xx=indgen(7)*60 -180
 yy=indgen( (yr2-yr1)/15 +1 )*15 +yr1
endif

xlonnm = strcompress(abs(xx),/remove_all)
xlatnm = strcompress(abs(yy),/remove_all)
for i=0,n_elements(xx)-1 do begin
 if xx(i) lt 0 and xx(i) ne -180 then xlonnm(i) = xlonnm(i)+'W'
 if xx(i) gt 0 and xx(i) ne  180 then xlonnm(i) = xlonnm(i)+'E'
endfor
for j=0,n_elements(yy)-1 do begin
 if yy(j) lt 0 then xlatnm(j) = xlatnm(j)+'S'
 if yy(j) gt 0 then xlatnm(j) = xlatnm(j)+'N'
endfor
none  = replicate(' ',13)

nq=n_elements(levels)
q =fltarr(nq,nq)
labels=strarr(nq)

if n_elements(colorindex) eq 0 then begin
 if nq eq 12 then color_ind=[30,70,86,104,144,180,194,202,212,250,10,10]
;if nq eq 11 then color_ind=[30,70,86,104,144,180,194,202,212,250,10]
 if nq eq 11 then color_ind=[30,70,86,104,144,180,194,    212,250,10,10]
 if nq eq 10 then color_ind=[30,70,86,104,144,180,194,    212,250,10,10]
 if nq eq 9 then color_ind=[30,70,104,144,180,194,    212,250,10,10]
endif else begin
 color_ind=colorindex
endelse

set_plot,'ps'
ict=39
loadct,ict
if region ne 'nor' then begin
 x1=0.07 & x2=0.94 & y1=0.13 & y2=0.95
 xc1=x1 & xc2=x2 & yc1=0.04 & yc2=0.06
endif
if region eq 'nor' then begin
 x1=0.07 & x2=0.94 & y1=0.74 & y2=0.95
 xc1=x1 & xc2=x2 & yc1=y1-0.07 & yc2=yc1+0.02
endif
!psym=0
!p.font=0
device,/symbol,font_index=4
;device,/portrait,bits=8
;device,/portrait,yoffset=1.,ysize=9.,/inches
device,/landscape,bits=8
device,/landscape,xoffset=0.5, yoffset=10.5, xsize=10, ysize=7.,/inches
device,filename=figname+'.ps',/color

if max(data(*,*,6)) gt max(levels) then levels(nq-1)=max(data(*,*,6))
for l=0,nq-1 do begin
 if levels(l) lt 0.0001 and levels(l) gt 0 then $
 labels(l)=string(levels(l),format='(e10.1)')
 if levels(l) ge 0.0001 and levels(l) lt 0.001 then $
 labels(l)=string(levels(l),format='(f6.4)')
 if levels(l) ge 0.001 and levels(l) lt 0.01 then $
 labels(l)=string(levels(l),format='(f5.3)')
 if levels(l) ge 0.01 and levels(l) lt 0.1 then $
 labels(l)=string(levels(l),format='(f4.2)')
 if levels(l) ge 0.1 and levels(l) lt 1 then $
 labels(l)=string(levels(l),format='(f4.2)')
 if levels(l) ge 1. and levels(l) lt 10. then $
 labels(l)=string(levels(l),format='(f3.1)')
 if levels(l) ge 10 and levels(l) lt 100. then $
 labels(l)=string(levels(l),format='(i2)')
 if levels(l) ge 100 and levels(l) lt 1000. then $
 labels(l)=string(levels(l),format='(i3)')
 if levels(l) ge 1000 and levels(l) lt 10000. then $
 labels(l)=string(levels(l),format='(i4)')
 if levels(l) ge 10000 and levels(l) lt 100000. then $
 labels(l)=string(levels(l),format='(i5)')
endfor

;set_viewport,x1,x2,y1,y2
;contour,data,longrid,latgrid,/cell_fill,levels=levels, $

;set_viewport,x1,x2,y1,y2
set_viewport
!x.omargin=[3,3]
!y.omargin=[2,2]
!p.charsize=1.2
!p.charthick=5
!p.multi=[0,1,1,1,0]

;gocart_level
;slev=['992','970','929','867','787','697','600','510','433','369',$
;      '313','266','226','192','163','139','118','100','85','67',  $
;	  '48','34','24','14','6.7','2.9','1.1','0.41','0.14','0.04']
slev=['1000','681','464','316','215','146','100','68','46','31',   $
      '21','14','10','6.8','4.64','3.16','2.15','1.46','1.00','0.681', $
	  '0.464','0','0','0','0','0','0','0','0','0','0','0','0','0','0','0','0']
nlev=n_elements(slev)

;for ilev=0,nlev-1 do begin
for ilev=6,6 do begin
map_set,/continent,/cyl,/noerase,/advance, $
title='MLS CO (ppbv) '+slev(ilev)+' hpa '+title , $
limit=[yr1,xr1,yr2,xr2],mlinethick=2,color=0
contour,data(*,*,ilev),longrid,latgrid,/overplot,/nodata,levels=levels, $
c_colors=color_ind,/fol, $ 
xrange=[xr1,xr2], yrange=[yr1,yr2], xstyle=1,ystyle=1, $
xticklen= -0.02, yticklen= -0.01, $
xtickname=xlonnm, ytickname=xlatnm, $
xticks=n_elements(xlonnm)-1, xminor=6, $
yticks=n_elements(xlatnm)-1, yminor=6, charsize=1.2,color=0
;map_set,/continent,/cyl,/noerase,$;/noborder, $
;limit=[yr1,xr1,yr2,xr2],mlinethick=2,color=0

 nlon=n_elements(longrid)
 nlat=n_elements(latgrid)
 ;dlon=2.5
 ;dlat=1.25
 icolor2d=intarr(nlon,nlat)
 for ilon=0,nlon-1 do begin
  ;print, 'ilon',ilon
 for ilat=0,nlat-1 do begin
   ;print, 'ilat',ilat
    if(data(ilon,ilat,ilev) lt levels(0)) then icolor2d(ilon,ilat)=255
     ;print, 'here'
   for l=0,nq-2 do begin
     ;print, 'here2'
    if(data(ilon,ilat,ilev) ge levels(l) and data(ilon,ilat,ilev) lt levels(l+1)) $
      then icolor2d(ilon,ilat)=color_ind(l)
   endfor
    ;print, 'here3'
    if(data(ilon,ilat,ilev) ge levels(nq-1)) then icolor2d(ilon,ilat)=0
   x0=longrid(ilon)-dlon/2.
    if(x0 lt -180.) then x0=x0+360.
    if(x0 gt  180.) then x0=x0-360.
   x1=longrid(ilon)+dlon/2.
    if(x1 lt -180.) then x1=x1+360.
    if(x1 gt  180.) then x1=x1-360.
   y0=latgrid(ilat)-dlat/2.
    if(y0 lt -90.) then y0=-90.
    if(y0 gt  90.) then y0= 90.
   y1=latgrid(ilat)+dlat/2.
    if(y1 lt -90.) then y1=-90.
    if(y1 gt  90.) then y1= 90.
    ;print, x0,x1,y0,y1
   polyfill,[x0,x0,x1,x1], [y0,y1,y1,y0],color=icolor2d(ilon,ilat)
 endfor
 endfor
map_continents,/continents,/coasts,/countries,mlinethick=2,color=0
 ;map_set,/continent,/cyl,/noerase, $;/noborder, $
 ;limit=[yr1,xr1,yr2,xr2],color=0
 ;if region eq 'usa' or region eq 'cal' then $
 ;map_set,/continent,/usa,/cyl,/noerase, $;/noborder, $
 ;limit=[yr1,xr1,yr2,xr2],color=0
endfor

for l = 0,nq-1 do q(l,*)=levels(l)
set_viewport,xc1,xc2,yc1,yc2
contour,q,/fil,levels=levels, $
c_colors=color_ind, $
/noerase, $
xstyle=1,ystyle=1, $
xticklen=0., $
yticklen=0., $
xtickname=labels, ytickname=replicate(' ',5), $
xticks=n_elements(levels)-1, xminor=1, $
yticks=1, yminor=1, charsize=1.2,color=0
contour,q,/ov,levels=levels,color=0

;xyouts, 0.5, 0.97, title,/normal
device,/close

end