pro read_f90dat
;
; Read and plot epot or fac from ascii file written by test_w05sc.F90
; write(lu,"(2i4)") nlon,nlat
; write(lu,"(6f8.2)") mlt
; write(lu,"(6f18.13)") mlat
; write(lu,"(6e12.4)") [epot or fac]
;
; file = 'wei05sc_epot_f90.dat'
  file = 'wei05sc_fac_f90.dat'
  openr,lu,file,/get_lun
  readf,lu,nmlt,nlat
  mlt = fltarr(nmlt)
  mlat = fltarr(nlat)
  readf,lu,mlt
  print,format="('nmlt=',i4,' mlt=',/,(10f8.3))",nmlt,mlt
  readf,lu,mlat
  print,format="('nlat=',i4,' mlat=',/,(10f8.3))",nlat,mlat
  data = fltarr(nmlt,nlat)
  readf,lu,data
  spval = 1.e36
  data(where(data eq spval)) = 0.
  fmin = min(data) & fmax = max(data)
  close,lu
  free_lun,lu
  print,'Completed read of file ',file

  lon = fltarr(nmlt)
  for i=0,nmlt-1 do lon[i] = mlt[i]*15. ; for plotting

  if (keyword_set(ps)) then begin
    set_plot,'PS'
  endif else begin
    set_plot,'X'
  endelse
;
; Center of projection:
  projcen = [90.,0.] ; lat,lon of proj center
  map_stperim = 50.  ; perimeter lat
;
; Position in normalized coords:
;
  xll = .15 & yll = .15 & xur = .85 & yur = .85
  mappos = [xll,yll,xur,yur]
; print,'mappos (xll,yll,xur,yur) =',mappos

  map_set,projcen[0],projcen[1],/stereo,position=mappos,/noborder,$
    limit=[map_stperim,-180.,90.,180.]

  fmin = min(data) & fmax = max(data)
  print,'data fmin,max=',fmin,',',fmax
;
; Hardwire contour levels and interval:
;
; cmin = -60. & cmax = 60. & cint = 10. ; hardwire min,max,int
; cmin = fix(fmin) & cmax = fix(fmax)   ; use fmin,fmax
  cmin = fmin & cmax = fmax
  cint = (cmax-cmin)/10.                ; make 10 levels
  nlevels = fix((cmax-cmin)/cint)
  print,'cmin=',cmin,' cmax=',cmax,' cint=',cint,' nlevels=',nlevels
  clevels = fltarr(nlevels)
  for i=0,nlevels-1 do clevels[i] = cmin+i*cint
;
; Negative contours are dashed, positive solid,
; and 0 is long dashes or dotted.
;
  clinestyle = intarr(nlevels)
  clinestyle[*] = 0 ; default solid
  for i=0,nlevels-1 do begin
    if clevels[i] lt 0. then clinestyle[i] = 2 ; dashed
    if clevels[i] eq 0. then clinestyle[i] = 1 ; dotted
;   if clevels[i] eq 0. then clinestyle[i] = 5 ; long dashes
  endfor

  contour,data,lon,mlat,/follow,xstyle=1,ystyle=1,$
    position=mappos,/noerase,/overplot,$
    levels=clevels,c_linestyle=clinestyle

  continents = 0
  grid = 0
  if continents eq 1 then map_continents
  if grid eq 1 then map_grid,/label

  labpol_chsize = 1.0
  mag = 1
  hem = 'N'
; fixpolar = 'LON'
  fixpolar = 'SLT'
  labpol,mappos,hem,labpol_chsize,fixpolar,mag=mag

  chsize = 2.0
  if (strpos(file,'epot') ne -1) then $
    title = 'Weimer 2005 Electric Potential (kV)' $
  else $
    title = 'Weimer 2005 Field-Aligned Current'

  xyouts,xll+.5*(xur-xll),yur+.10,title,align=.5,/norm,$
    charsize=chsize

  chsize = 1.2
; title = 'BZ='+string(format="(f8.3)",bz[isolar])+$
;        ' BY='+string(format="(f8.3)",by[isolar])+$
;        ' SWVEL='+string(format="(f8.3)",swvel)+$
;        ' SWDEN='+string(format="(f8.3)",swden)
; xyouts,xll+.5*(xur-xll),yur+.05,title,align=.5,/norm,$
;   charsize=chsize

  chsize = 1.2
  minmaxlab = 'MIN,MAX = '+string(format="(g10.4)",fmin)+','+$
                           string(format="(g10.4)",fmax)+'  '+$
             'INTERVAL = '+string(format="(g10.4)",cint)
  xyouts,xll+.5*(xur-xll),yll-.12,minmaxlab,align=.5,/norm,$
    charsize=chsize

  if (keyword_set(ps)) then device,/close_file

end