pro rdgpi,ps=ps
;
; Read file made by fgrep "Step" tgcm23_task0.out, when gpi==1, e.g.:
;  0:Step     1 of  2400 mtime  80, 0, 3, 0 f107d=128.15 power= 22.05 ctpoten= 58.28
; Make line plots of f107d, power, and ctpoten
;
;openr,lu,'gpi.d80-85.step.out',/get_lun
;openr,lu,'gpi.d85-90.step.out',/get_lun
;openr,lu,'gpi.d80-90.step.out',/get_lun
;openr,lu,'gpi.step.out',/get_lun
;openr,lu,'gpi_ncep.step.out',/get_lun
;openr,lu,'gpi.step.out',/get_lun
;openr,lu,'roble.step.out',/get_lun
 openr,lu,'gpi.step.cmq.out',/get_lun
mxgpi = 1000
f107d = fltarr(mxgpi)
hpower = fltarr(mxgpi)
ctpoten = fltarr(mxgpi)
mtime = intarr(3,mxgpi)
dmtime = fltarr(mxgpi) ; decimal day (from mtime)
iday=-1
ihr=-1
imin=-1
i = 0
;0 Step    10 of  2400 mtime 100, 0,30, 0 f107d=172.64 power= 24.90 ctpoten= 61.63 wcsecs/step=  1.89
while not eof(lu) do begin
  readf,lu,$
;   format="(31x,i3,1x,i2,1x,i2,10x, f6.2,7x,f6.2,9x,f6.2)",$
;     iday,ihr,imin,f107d_rd, hpower_rd, ctpoten_rd
    format="(28x,i4,1x,i2,1x,i2,10x, f6.2,7x,f6.2,9x,f6.2)",$
      iday,ihr,imin,f107d_rd, hpower_rd, ctpoten_rd
  mtime(0,i) = iday
  mtime(1,i) = ihr
  mtime(2,i) = imin
  dmtime(i) = float(iday)+(float(ihr)/24.)+(float(imin)/(24.*60.))
  f107d(i) = f107d_rd 
  hpower(i) = hpower_rd 
  ctpoten(i) = ctpoten_rd 
  i = i+1
end
ngpi = i
free_lun,lu
;
if (keyword_set(ps)) then begin
  set_plot,'ps'
  device,file='rdgpi.ps',xoffset=1.,yoffset=3.
endif else begin
  set_plot,'x'
endelse
;
;print,'f107d=',f107d(0:ngpi-1)
;print,'hpower=',hpower(0:ngpi-1)
;print,'ctpoten=',ctpoten(0:ngpi-1)
;
print,format="('ngpi=',i5,' mtime=',3i4,' to ',3i4)",$
  ngpi,mtime(0:2,0),mtime(0:2,ngpi-1)
title = string(format="(' mtime ',3i4,' to ',3i4)",$
  mtime(0:2,0),mtime(0:2,ngpi-1))
xtitle = string(format="('Model day (ngpi = ',i5,')')",ngpi)
;if (keyword_set(ps)) then device,font_size=10.,xsize=19.,ysize=20.
!p.multi = [0,2,2]
;
; f107d:
rmin = min(f107d(0:ngpi-1))
rmax = max(f107d(0:ngpi-1))
ytitle = 'F107D (' + string(format="(f5.1)",rmin) + ' to ' + $
                     string(format="(f5.1)",rmax) + ')'
plot,dmtime(0:ngpi-1),f107d(0:ngpi-1),ytitle=ytitle,$
  xtitle=xtitle,title=title,$
; yrange=[100.,150.]
  yrange=[130.,175.]
; yrange=[65.,85.]
;xcurs = 0. & ycurs = 0.
;if (not keyword_set(ps)) then cursor,xcurs,ycurs,1
;
; hpower:
rmin = min(hpower(0:ngpi-1))
rmax = max(hpower(0:ngpi-1))
ytitle = 'HPOWER (' + string(format="(f5.1)",rmin) + ' to ' + $
                      string(format="(f5.1)",rmax) + ')'
plot,dmtime(0:ngpi-1),hpower(0:ngpi-1),ytitle=ytitle,$
  xtitle=xtitle,title=title
;if (not keyword_set(ps)) then cursor,xcurs,ycurs,1
;
; ctpoten:
rmin = min(ctpoten(0:ngpi-1))
rmax = max(ctpoten(0:ngpi-1))
ytitle = 'CTPOTEN (' + string(format="(f5.1)",rmin) + ' to ' + $
                       string(format="(f5.1)",rmax) + ')'
plot,dmtime(0:ngpi-1),ctpoten(0:ngpi-1),ytitle=ytitle,$
  xtitle=xtitle,title=title
if (keyword_set(ps)) then device,/close
end