;-------------------------------------------------------------------------
pro addmapgrid,map,field
;
; Add continental outlines and/or map grid if not magnetic grid.
;
if (map.continents eq 1 or map.grid eq 1) and field.grid_type ne 'magnetic' $
  then begin
  if map.continents eq 1 then map_continents,color=map.continents_color
  if map.grid eq 1 then map_grid,/label,color=map.grid_color
endif else if (map.continents eq 1 or map.grid eq 1) and $
  field.grid_type eq 'magnetic' then begin
  print,'WARNING: cannot plot continents or map grid on magnetic grid.'
  map.continents = 0
  map.grid = 0
endif
end
;-----------------------------------------------------------------------
;
pro pltce,info,animate=animate
;
file = *info.pfile
map = *info.pmap
fields = *info.fields
field = map.field
data = *map.data
lons = *map.lons & nlon = n_elements(lons)
lats = *map.lats & nlat = n_elements(lats)

;
;print,' '
;print,'pltce:'
;print,'field.name = ',field.name
;print,'map.projection = ',map.projection
;print,'map.zvert=',map.zvert
;print,'map.mtime=',map.mtime
;print,'map.missing_value=',map.missing_value
;print,'map.stperim=',map.stperim
;print,'nlon=',nlon,' lons=' & print,lons
;print,'nlat=',nlat,' lats=' & print,lats
;print,'data min,max = ',min(data),',',max(data)
;
; Write binary file for debug:
;file='for_rsi.dat'
;openw,lu,file,/get_lun
;writeu,lu,nlon
;writeu,lu,lons
;writeu,lu,nlat
;writeu,lu,lats
;writeu,lu,data
;close,lu
;free_lun,lu
;print,'Wrote binary file ',file,' lu=',lu
;
; Check for all missing value:
;
nmissing = 0
text = 1 & if keyword_set(animate) then text = 0
nodata = check_missing(data,float(map.missing_value),field.name,text=text,$
  nmissing=nmissing)
if nodata eq 1 then return
missing_value = float(map.missing_value)
if nmissing gt 0 then print,'Global Maps: Field ',field.name,': Found ',$
  nmissing,' missing data points (out of ',n_elements(data),')'
;
; Set map position (normalized coords):
;
xll = .15 & yll = .48
xur = .95 & yur = yll+.5*(xur-xll)
mappos = [xll,yll,xur,yur]
dlat = lats[1]-lats[0]
;
; Geographic field lats start at +/- 87.5, whereas magnetic fields go to +/- 90.
;
if (lats[0]-dlat/2. eq -90.) and (lats[nlat-1]+dlat/2. eq +90.) then begin
  poledeg = (dlat/2.)/180.*(yur-yll) ; distance between top lat and pole
  conpos = [xll,yll+poledeg,xur,yur-poledeg]
; conpos = [xll,yll+poledeg,xur,yur-poledeg*2.5]
endif else begin
  conpos = [xll,yll,xur,yur]
endelse
xmid = .5*(xll+xur)
;
; Set up map projection:
;
; If censlt >= 0, then it is user-specified local time to fix at center
; of projection.
;
censlt = map.censlt ; local time (hrs) of center of lon axis of projection
;
plat=0 & plon=map.cenlon & rot = 0.

if(map.fixpolar eq 'LON') then censlt=-1 

if censlt ge 0. then begin
  uthrs = float(map.mtime[1])+float(map.mtime[2])/60.
  plon = fslt(censlt,uthrs,plon,'getlon')
; rot = censlt*15.
; print,'pltce: censlt=',censlt,' uthrs=',uthrs,' plon=',plon,' rot=',rot
endif
map_set,plat,plon,rot,/cylindrical,position=mappos,/noborder
;
; Titles, labels, tic marks, etc:
;
if field.grid_type eq 'geographic' then begin
  xtitle = 'LONGITUDE (DEGREES GEOGRAPHIC)'
  if censlt ge 0. then xtitle = 'SOLAR LOCAL TIME (HOURS)' 
  ytitle = 'LATITUDE (DEGREES GEOGRAPHIC)'
endif else begin
  xtitle = 'LONGITUDE (DEGREES MAGNETIC)'
  if censlt ge 0. then xtitle = 'MAG LOCAL TIME (HOURS)' 
  ytitle = 'LATITUDE (DEGREES MAGNETIC)'
endelse
xticlabs = ['-180','-120','-60','0','60','120','180']
if censlt ge 0. then begin
  sltlabs = fltarr(7)
  sltlabs[3] = censlt
  for i=2,0,-1 do begin
    sltlabs[i] = sltlabs[i+1]-4.
    if(sltlabs[i] ge 24.) then sltlabs[i]=sltlabs[i]-24
    if(sltlabs[i] lt 0.) then sltlabs[i]=sltlabs[i]+24
  endfor
  for i=4,6    do begin
    sltlabs[i] = sltlabs[i-1]+4.
    if(sltlabs[i] ge 24.) then sltlabs[i]=sltlabs[i]-24
    if(sltlabs[i] lt 0.) then sltlabs[i]=sltlabs[i]+24
  endfor
  xticlabs = strarr(7)
  for i=0,6 do begin
	xticlabs[i] = string(format="(f5.1)",sltlabs[i])
	xticlabs[i] = strcompress(xticlabs[i],/remove_all)
  endfor
; print,'pltce: censlt=',censlt,' xticlabs=' & print,xticlabs
endif
if(map.fixpolar eq 'LON') then begin
  lonlabs=fltarr(7)
  lonlabs[3]=map.cenlon
  for i=2,0,-1 do begin 
    lonlabs[i] = lonlabs[i+1]-60.
    if(lonlabs[i] gt 180.) then lonlabs[i]=lonlabs[i]-360
    if(lonlabs[i] lt -180.) then lonlabs[i]=lonlabs[i]+360
  endfor
  for i=4,6    do begin
    lonlabs[i] = lonlabs[i-1]+60.
    if(lonlabs[i] gt 180.) then lonlabs[i]=lonlabs[i]-360
    if(lonlabs[i] lt -180.) then lonlabs[i]=lonlabs[i]+360
  endfor
  xticlabs = strarr(7)
  for i=0,6 do begin
	xticlabs[i] = string(format="(f6.1)",lonlabs[i])
	xticlabs[i] = strcompress(xticlabs[i],/remove_all)
  endfor
  xtitle = 'LONGITUDE (DEGREES GEOGRAPHIC)
endif

yticlabs = ['-90','-60','-30','0','30','60','90']
;
fminmax,data,framemin,framemax,missing_value
if map.fmin ge map.fmax then begin
  fminmax,data,fmin,fmax,missing_value
  if map.fmin eq map.fmax then begin
    map.fmin=fmin & map.fmax=fmax
    if not keyword_set(animate) then $
      print,'Global CE maps reset map.fmin,max=',map.fmin,map.fmax
  endif else $
    if not keyword_set(animate) then $
      print,'Global CE maps did NOT reset min,max: map.fmin,max=',$
        map.fmin,map.fmax
endif else begin
  fmin = map.fmin & fmax = map.fmax
  if not keyword_set(animate) then $
    print,'Global CE maps using preset map.fmin,max=',map.fmin,map.fmax
endelse
;
; Byte-scale to an image. Make any missing data black:
;
if (map.plottype eq 'image_only') or (map.plottype eq 'image+contours') then begin
  missing_indices = where(data eq map.missing_value,nmissing)
  if nmissing gt 0 then begin
    print,'Field ',field.name,' CE map: Found ',nmissing,' missing data points.'
    ave = total(data(where(data ne map.missing_value)))/n_elements(data)
    data(missing_indices) = ave
    image = bytscl(data,top=!d.table_size-3,min=fmin,max=fmax)+1
    image(missing_indices) = 0
  endif else $
    image = bytscl(data,top=!d.table_size-3,min=fmin,max=fmax)+1
;
; Scale image up and warp image to map projection:
;
; As of 3/27/03: congrid + map_image results in offset contours and image
;   (contours do not match colors in the image, esp high north latitudes).
;   (This is apprently because of interpolating twice before contour)
;   Map_patch eliminates this problem, but results in a vertical dashed
;   black line at longitude=0. This is apparently due to a bug in map_patch.
;   A work-around is to reduce nlon by one with congrid before map_patch.
;   I will use this until the bugfix for map_patch is available. For
;   correspondence with RSI, see ~foster/tgcmproc/consult/image+contour.
; It is also necessary to use map_patch for irregular latitude grid with
;   fields on mag grid.
;
; imscale = 8
; image = congrid(image,nlon*imscale,nlat*imscale,/interp)
; image = map_image(image,startx,starty,xsize,ysize,$
;   latmin=lats[0],latmax=lats[nlat-1],$
;   lonmin=lons[0],lonmax=lons[nlon-1],compress=1)

; image = map_patch(image,lons,lats,xstart=startx,ystart=starty,$
;   xsize=xsize,ysize=ysize) ; dashed black line bug

  image_patch = congrid(image,nlon-1,nlat) ; bug fix hack
  image = map_patch(image_patch,dummy,lats,$
    xstart=startx,ystart=starty,xsize=xsize,ysize=ysize)
;
; Set up color bar (pro clrbar):
;
  barw = .90*(xur-xll) & barh = .03
  xbar = xmid-.5*barw & ybar = yll-barh-.1
  cblabels = [strcompress(string(fmin,format="(g10.4)"),/remove_all),$
    strcompress(string(0.5*(fmin+fmax),format="(g10.4)"),/remove_all),$
    strcompress(string(fmax,format="(g10.4)"),/remove_all)]
endif ; plottype is image_only or image+contours
;
; Field label:
;
if strlen(field.long_name) gt 0 then fieldlab = field.long_name $
  else fieldlab = field.name
if file.diffs le 0 then begin
  if strlen(field.units) gt 0 then $
    fieldlab = fieldlab + ' (' + strcompress(field.units,/remove_all) + ')'
endif else fieldlab = fieldlab + ' (' + field.difftype + ' ' + $
  strcompress(field.units,/remove_all) + ')'
;
; Top label (under field label):
;
uthrs = float(map.mtime[1])+float(map.mtime[2])/60.
toplab = string(format="('DAY = ',i3,'  UT = ',f5.2)",map.mtime[0],uthrs)
if field.nlev gt 0 then $
  toplab = toplab + string(format="('  ',a,' = ',f6.2)",$
    field.levshortname,map.zvert)
;
; Min,max label:
;
minmaxlab = 'MIN,MAX = '+string(format="(g10.4)",framemin)+','+$
                         string(format="(g10.4)",framemax)
if map.plottype ne 'image_only' then $
  minmaxlab = minmaxlab + string(format="(' INTERVAL = ',g10.4)",map.cint)
;
; Bottom label:
;
 botlab1 = strcompress(field.ncfile,/remove_all)
;botlab1 = strcompress(field.ncfile,/remove_all)+ $
; string(format="(' (',i3,',',i2,',',i2,')')",map.mtime[*])
;
; Contour levels were made by mkclevels.pro.
;
clevels = *map.clevels
nlevels = n_elements(clevels)
clinestyle = intarr(nlevels)
clinestyle[*] = 0 ; solid
for i=0,nlevels-1 do if clevels[i] lt 0. then clinestyle[i] = 2 ; dashed
c_charsize = 1.1
chsize = 1.2
;
; Display image and/or contour.
; Use position=mappos when drawing axes and labels, and use
;   position=conpos when drawing contours.
;
;
;Rotate data in longitude if necessary (i.e., fixed slt has been requested):
;
if(map.fixpolar eq 'SLT') then begin
  ix0 = ixfind_nearest(lons,0.)
  ixcen = ixfind_nearest(lons,plon)
  if ixcen lt ix0 then ishift = ix0-ixcen else ishift = -(ixcen-ix0)
; print,'censlt=',censlt,' ix0=',ix0,' ixcen=',ixcen,' shift data by ',ishift
  data = shift(data,ishift,0)
endif
if(map.fixpolar eq 'LON') then begin
  ix0 = ixfind_nearest(lons,0.)
  ixcen = ixfind_nearest(lons,plon)
  if ixcen lt ix0 then ishift = ix0-ixcen else ishift = -(ixcen-ix0)
; print,'censlt=',censlt,' ix0=',ix0,' ixcen=',ixcen,' shift data by ',ishift
  data = shift(data,ishift,0)
endif
case map.plottype of
;
; Display image only. 
; Use contour,/nodata to draw axes, ticmarks, labels:
;
  'image_only': begin
    tv,image,startx,starty,xsize=xsize,ysize=ysize
    addmapgrid,map,field
    contour,data,lons,lats,/noerase,/nodata,position=mappos,$
      xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
      xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
      yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
    clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.21,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.24,botlab1  ,/norm,align=0.5,charsize=chsize
  end
;
; Display image with contour overlay:
;
  'image+contours': begin
    tv,image,startx,starty,xsize=xsize,ysize=ysize
    addmapgrid,map,field
; axes:
    contour,data,lons,lats,/noerase,position=mappos,/nodata,$
      xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
      xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
      yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
; contours:
    contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
      /follow,c_charsize=c_charsize,levels=clevels,c_linestyle=clinestyle,$
      c_colors=[map.clineclr]
    clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.21,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.25,botlab1,/norm,align=0.5,charsize=chsize
  end
;
; Contour only (B&W):
;
  'monochrome_contours': begin
    addmapgrid,map,field
    contour,data,lons,lats,/noerase,/nodata,position=mappos,$
      xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
      xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
      yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
    contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
      /follow,c_charsize=c_charsize,levels=clevels,c_linestyle=clinestyle,$
      c_colors=[map.clineclr]
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.11,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.15,botlab1,/norm,align=0.5,charsize=chsize
  end
;
; Color fill contour:
;
  'colorfill_contours': begin
    contour,data,lons,lats,/noerase,position=mappos,/nodata,$
      xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
      xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
      yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
    contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
      c_charsize=c_charsize,levels=clevels,/follow,/fill
    contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
      c_charsize=c_charsize,levels=clevels,/follow,c_linestyle=clinestyle,$
      c_colors=[map.clineclr]
    map_set,/cylindrical,position=mappos,/noborder,/noerase
    addmapgrid,map,field
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.11,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.15,botlab1,/norm,align=0.5,charsize=chsize
  end
  else: print,'>>> Global CE maps: unknown plot type ',map.plottype
endcase
;
; map.fmin,max and map.continents may have been changed, 
; so update state:
;
*info.pmap = map
end