;
pro plt_integ,ps=ps
;
; Read and plot height-integrated global means calculated by tgcmproc.
; Read ascii text file made from grep on "glbmean vertical integ" in 
;   tgcmproc_f90 stdout prints from mkutvert.F.
;
; Sample first day from glbmdifj.dat (obtained by grep of glbmdifj.out):
;
;Field TN       mtime=  10   0   0 glbmean vertical integ (pert) =   0.6387E+11 total=  0.3258E+30
;Field TN       mtime=  10   0   0 glbmean vertical integ (cntr) =   0.5390E+11 total=  0.2749E+30
;Field H2O      mtime=  10   0   0 glbmean vertical integ (pert) =   0.8886E+21 total=  0.4533E+40
;Field H2O      mtime=  10   0   0 glbmean vertical integ (cntr) =   0.1311E+19 total=  0.6688E+37
;Field NOZ      mtime=  10   0   0 glbmean vertical integ (pert) =   0.1949E+19 total=  0.9940E+37
;Field NOZ      mtime=  10   0   0 glbmean vertical integ (cntr) =   0.1239E+17 total=  0.6322E+35
;Field O3       mtime=  10   0   0 glbmean vertical integ (pert) =   0.1331E+19 total=  0.6787E+37
;Field O3       mtime=  10   0   0 glbmean vertical integ (cntr) =   0.1297E+19 total=  0.6615E+37
;Field NO       mtime=  10   0   0 glbmean vertical integ (pert) =   0.1257E+17 total=  0.6411E+35
;Field NO       mtime=  10   0   0 glbmean vertical integ (cntr) =   0.9105E+16 total=  0.4644E+35
;Field NO2      mtime=  10   0   0 glbmean vertical integ (pert) =   0.3749E+16 total=  0.1912E+35
;Field NO2      mtime=  10   0   0 glbmean vertical integ (cntr) =   0.3437E+16 total=  0.1753E+35
;
; file = 'glbmdifj.dat'
; file = 'glbmdifj_d10-360.dat'
; file = 'glbmdifq.dat'
  file = 'glbmdifr.dat'
;
  openr,lu,file,/get_lun
  print,'Opened file ',file,' lu=',lu
  fields = ['TN','H2O','NOZ','O3','NO','NO2']
  nflds = n_elements(fields)
  print,'Assuming ',nflds,' fields: ',fields
;
; Read 2 lines (pert,cntr) per field per day.
; Determine number of days on the file: 
;
  line = ''
  ndays = 0                         ; number of days on the file
  mtime = intarr(3)
  while (not eof(lu)) do begin
    for ifld = 0,nflds-1 do begin
      for i=0,1 do begin            ; pert,cntr
        readf,lu,line
      endfor
    endfor
    ndays = ndays+1
  end
  print,'ndays=',ndays
  point_lun,lu,0
;
; Read data:
;
  finteg_pert = dblarr(ndays,nflds)  ; perturbed integrations
  finteg_cntr = dblarr(ndays,nflds)  ; control integrations
  finteg_difs = dblarr(ndays,nflds)  ; difference fields
  finteg = double(0.)
  days = intarr(ndays)
  mtimes = intarr(3,ndays)
  for iday=0,ndays-1 do begin
    for ifld = 0,nflds-1 do begin
      for i=0,1 do begin            ; pert,cntr
        readf,lu,format="(21x,3i4,52x,e12.4)",mtime,finteg
        if (i eq 0) then begin
          finteg_pert[iday,ifld] = finteg 
        endif else begin
          finteg_cntr[iday,ifld] = finteg 
        endelse
        days[iday] = mtime[0]
        mtimes[*,iday] = mtime[*]
      endfor ; i=0,1
    endfor ; ifld=0,nflds-1
;   print,'mtime ',mtimes[*,iday]
;   for ifld=0,nflds-1 do $
;     print,'  ',fields[ifld],' pert=',finteg_pert[iday,ifld],' cntr=',finteg_cntr[iday,ifld]
  endfor ; iday=0,ndays-1
  free_lun,lu
;
; Set up for plots:
;
  if (keyword_set(ps)) then begin
    set_plot,'PS'
    psfile = file+'.ps'
    device,file=psfile
    print,'Set plot to PS: will make plot file ',psfile
  endif else begin
    set_plot,'X'
    print,'Set plot to X'
  endelse
  filelab = 'File = ' + file
  case file of 
   'glbmdifj.dat': begin
      filelab = filelab + ' (100% case)'
   end
   'glbmdifj_d10-360.dat': begin
      filelab = filelab + ' (100% case)'
   end
   'glbmdifq.dat': begin
      filelab = filelab + ' (10% case)'
   end
   'glbmdifr.dat': begin
      filelab = filelab + ' (1% case)'
   end
   else: print,'>>> WARNING unknown file name ',file,' for case designation.'
  endcase 
  chsize = 1.2
  re = 6371.e5
  pi = 3.14159
  earth_area = 4.*pi*re^2
  print,'earth_area = 4.*pi*re^2 = ',earth_area
;
; Field loop:
;
  for ifld = 0,nflds-1 do begin
;
; Plot perturbed case:
;
    title_pert = 'Height-integrated Global Mean ' + fields[ifld] + ' (perturbed case)'
    xtitle='Model Day'
    ytitle='Height-integrated ' + fields[ifld] + ' (perturbed)'
    f = finteg_pert[*,ifld]
    if (fields[ifld] eq 'TN') then f = f / earth_area
    fmin = min(f) & fmax = max(f)
    plot,mtimes[0,*],f,title=title_pert,xtitle=xtitle,$
      ytitle=ytitle,yrange=[fmin,fmax]
    xyouts,0.55,-.05,filelab,/norm,align=0.5,charsize=chsize
    minmaxlab = 'Min,max = ' + string(fmin) + ', ' + string(fmax)
    xyouts,0.55,-.10,minmaxlab,/norm,align=0.5,charsize=chsize
;
; Plot control case:
;
    title_cntr = 'Height-integrated Global Mean ' + fields[ifld] + ' (control case)'
    xtitle='Model Day'
    ytitle='Height-integrated ' + fields[ifld] + ' (control)'
    f = finteg_cntr[*,ifld]
    if (fields[ifld] eq 'TN') then f = f / earth_area
    fmin = min(f) & fmax = max(f)
    plot,mtimes[0,*],f,title=title_cntr,xtitle=xtitle,ytitle=ytitle,$
      yrange=[fmin,fmax]
    xyouts,0.55,-.05,filelab,/norm,align=0.5,charsize=chsize
    minmaxlab = 'Min,max = ' + string(fmin) + ', ' + string(fmax)
    xyouts,0.55,-.10,minmaxlab,/norm,align=0.5,charsize=chsize
;
; Plot diffs:
;
    title_difs = 'Height-integrated Global Mean ' + fields[ifld]
    title_difs = title_difs + ' (difference fields)'
    xtitle='Model Day'
    ytitle='Height-integrated ' + fields[ifld] + ' (differences)'
    finteg_difs = finteg_pert[*,ifld] - finteg_cntr[*,ifld] ; whole array diffs
    if (fields[ifld] eq 'TN') then $
      finteg_difs = (finteg_pert[*,ifld] / earth_area) - $
                    (finteg_cntr[*,ifld] / earth_area)
    fmin = min(finteg_difs) & fmax = max(finteg_difs)
    plot,mtimes[0,*],finteg_difs,title=title_difs,$
      xtitle=xtitle,ytitle=ytitle,yrange=[fmin,fmax]
    xyouts,0.55,-.05,filelab,/norm,align=0.5,charsize=chsize
    minmaxlab = 'Min,max = ' + string(fmin) + ', ' + string(fmax)
    xyouts,0.55,-.10,minmaxlab,/norm,align=0.5,charsize=chsize
    plots,[days[0],days[ndays-1]],[0.,0.],linestyle=2
;
; End field loop:
  endfor ; ifld=0,nflds-1
;
  if (keyword_set(ps)) then begin
    device,/close
  endif
end