function param
  COMPILE_OPT IDL2,HIDDEN
; Dummy to keep IDL from complaining. The real param() routine will be
; compiled below
end
;
function param2
COMPILE_OPT HIDDEN
; Dummy to keep IDL from complaining. The real param2() routine will be
; compiled below
end
;
spawn,'pwd'
default, quiet, 0
default, proc, 0
default, datatopdir, 'data'
default, varfile, 'var.dat'
default, dimfile, 'dim.dat'
;datadir = datatopdir+'/proc'+str(proc)
datadir = datatopdir
;
mx=0L & my=0L & mz=0L & nvar=0L & naux=0L
prec=''
nghostx=0L & nghosty=0L & nghostz=0L
quiet=3
;
close,1
openr,1,datadir+'/'+dimfile
readf,1,mx,my,mz,nvar,naux
readf,1,prec
readf,1,nghostx,nghosty,nghostz
close,1
;
prec = (strtrim(prec,2))        ; drop leading zeros
prec = strmid(prec,0,1)
if ((prec eq 'S') or (prec eq 's')) then begin ; single precision
  one = 1.e0
  nl2idl_d_opt = ''
endif else if ((prec eq 'D') or (prec eq 'd')) then begin ; double precision
  one = 1.D0
  nl2idl_d_opt = '-d'
endif else begin
  if (quiet le 4) then print, "prec = `", prec, "' makes no sense to me"
  STOP
endelse
;
  lhydro=1
  ldensity=1
  lentropy=0
  lmagnetic=1
  lionization_fixed=0
  leos_fixed_ionization=0
;
@readstartpars
;
nintervals=5
default,tit,''
default,t1,10
default,t2,1e9
default,kf,5.
default,yrSt,[0,15]
@param
quiet=0
te=max(tt) < t2
tlast=te-.25*(te-t1)
good=where(tt ge t1 and tt le t2)
good1=where(tt ge tlast and tt le t2)
urms=urms>1e-9
;St=-uxbm/(jxbxbm+uxDxuxbm)*urms*kf
alp=uxbm/urms
oumn=oum/(kf*urms^2)
jbmn=jbm/(kf*urms^2)
u2u13mn=u2u13m/(kf*urms^2)
b2b13mn=b2b13m/(kf*urms^2)
if quiet > 0 then begin
  ;print,'alp/urms,St=',mean(alp(good)),mean(St(good)),fo='(a,f8.5,f10.3)'
  print,'<o.u>, <j.b>=',mean(oum(good)),mean(jbm(good)),fo='(a,f8.4,f8.4)'
  print,'<u2u13m>, <b2b13m>=',mean(u2u13m(good)),mean(b2b13m(good)),fo='(a,f8.5,f8.5)'
  print,'urms=',mean(urms(good)^2)^.5,fo='(a,f7.3)'
  if i_uxjm ne 0 then print,mean(uxjm(good))
endif
;
;  plot 6 panels
;
!p.multi=[0,2,3]
!p.charsize=2
;
;  negative alpha value
;
plot,tt,alp,ytit='!6-!7a!6 / !8u!6!drms!n',yr=minmax(alp(good))
oplot,tt(good),alp(good),col=122
oplot,tt(good1),alp(good1),col=188
meanerror,alp(good),alpm,alpdn,alpup,nintervals,err=alperr
oplot,tt(good),tt(good)*0.+alpm,col=188
oplot,tt(good),tt(good)*0.+alpdn
oplot,tt(good),tt(good)*0.+alpup
;
;  urms
;
plot,tt,urms,ytit='urms',yr=minmax(urms(good))
oplot,tt(good),urms(good),col=122
oplot,tt(good1),urms(good1),col=188
meanerror,urms(good),urmsm,urmsdn,urmsup,nintervals,err=urmserr
oplot,tt(good),tt(good)*0.+urmsm,col=188
oplot,tt(good),tt(good)*0.+urmsdn
oplot,tt(good),tt(good)*0.+urmsup
;
;  oum
;
plot,tt,oumn,ytit='oumn',yr=minmax(oumn(good))
oplot,tt(good),oumn(good),col=122
oplot,tt(good1),oumn(good1),col=188
meanerror,oumn(good),oumnm,oumndn,oumnup,nintervals,err=oumnerr
oplot,tt(good),tt(good)*0.+oumnm,col=188
oplot,tt(good),tt(good)*0.+oumndn
oplot,tt(good),tt(good)*0.+oumnup
;
;  jbm
;
plot,tt,jbmn,ytit='jbm',yr=minmax(jbmn(good))
oplot,tt(good),jbmn(good),col=122
oplot,tt(good1),jbmn(good1),col=188
meanerror,jbmn(good),jbmnm,jbmndn,jbmnup,nintervals,err=jbmnerr
oplot,tt(good),tt(good)*0.+jbmnm,col=188
oplot,tt(good),tt(good)*0.+jbmndn
oplot,tt(good),tt(good)*0.+jbmnup
;
;  u2u13m
;
plot,tt,u2u13mn,ytit='u2u13m',yr=minmax(u2u13mn(good))
oplot,tt(good),u2u13mn(good),col=122
oplot,tt(good1),u2u13mn(good1),col=188
meanerror,u2u13mn(good),u2u13mnm,u2u13mndn,u2u13mnup,nintervals,err=u2u13mnerr
oplot,tt(good),tt(good)*0.+u2u13mnm,col=188
oplot,tt(good),tt(good)*0.+u2u13mndn
oplot,tt(good),tt(good)*0.+u2u13mnup
;
;  b2b13m
;
plot,tt,b2b13mn,ytit='b2b13m',yr=minmax(b2b13mn(good))
oplot,tt(good),b2b13mn(good),col=122
oplot,tt(good1),b2b13mn(good1),col=188
meanerror,b2b13mn(good),b2b13mnm,b2b13mndn,b2b13mnup,nintervals,err=b2b13mnerr
oplot,tt(good),tt(good)*0.+b2b13mnm,col=188
oplot,tt(good),tt(good)*0.+b2b13mndn
oplot,tt(good),tt(good)*0.+b2b13mnup
;
fo='(f7.4,2f8.5)'
if quiet > 0 then begin
  print,mean(alp(good1)),mean(u2u13m(good1)),mean(b2b13m(good1)),fo=fo
  print,mean(alp(good1)),mean(oum(good1)),mean(jbm(good1)),fo=fo
  print,'max(tt)-t1=',max(tt)-t1
  print,'max(tt)=',max(tt)
endif
!p.multi=0
;
fo='(f4.2,e8.1,f4.1,2f6.3,2f8.4,8f7.3,2i5,1x,a)'
B0=par2.B_ext(2)
eta=par2.eta
tend=long(max(tt))
delt=long(tend-t1)
openw,9,'pstrouhal.tmp'
print,   B0,eta,kf,urmsm,urmserr,alpm,alperr,oumnm,oumnerr,jbmnm,jbmnerr,u2u13mnm,u2u13mnerr,b2b13mnm,b2b13mnerr,tend,delt,tit,fo=fo
printf,9,B0,eta,kf,urmsm,urmserr,alpm,alperr,oumnm,oumnerr,jbmnm,jbmnerr,u2u13mnm,u2u13mnerr,b2b13mnm,b2b13mnerr,tend,delt,tit,fo=fo
close,9
;
spawn,'cat pstrouhal.tmp >> ../pstrouhal.dat'
END