;***************************************************************************
;***************************************************************************
;Read SGI files on DEC
;***************************************************************************
;***************************************************************************

;procedure file,var,fpln,cdataz,cdatas
;are converted such that they byteswap the SAV- data before they are stored
;this for using machines that operates on byteswapped order relative to the stored data

;***************************************************************************
;
;  Write variables in a plane, unpacked, ascii or binary packed
;
;***************************************************************************
pro wpln,f,j,i,pos
	common cmov,cvar,cpos,cfirst,clast,cstep,csize,cfile,ctmp
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
	common cplane,nc

	if n_params() eq 0 then begin
		print,'wpln(f,itime,ivar,n) [0 indexed]'
		return
	end
	if n_params() lt 3 then pos=cpos
	cpos=pos
	ntot=nx*ny*nz
	nxy=nx*ny

	; unpacked scr file
	if (incr gt 10) then begin   
		a=assoc(1,fltarr(nx,ny),irec(j*incr+i)+pos*nxy*4L)
		a(0)=f
		return

	;  ascii packed sav file, ipack=1
	endif else if (ipack eq 1) then begin
		nc=((ntot+3)/4)*2
		a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L)
		b=a(0)
		c1=string(b(0:15))
		c2=string(b(16:31))
		vmin=float(c1)
		vscl=float(c2)
		cscl=1./vscl
		print,b,min(f),vmin,vscl
		print,'wr.offset=',irec(j*incr+i)+pos*nxy*4L,nx*ny
		a=assoc(1,intarr(nx,ny),irec(j*incr+i)+pos*nxy*4L)
		m15=2L^15
  		a(0)=fix((f-vmin)*cscl+0.5)-m15
  		return
	endif else if (ipack eq 2) then begin
		ntot=nx*ny*nz
		nc=(ntot+1)/2
		a=assoc(1,fltarr(2),irec(j*incr+i)+nc*4L)
		b=a(0)
		vmin=b(0)
		vscl=b(1)
		vscl=1./vscl
		print,b,vmin,vscl
		a=assoc(1,intarr(nx,ny),irec(j*incr+i)+pos*nxy*2L)
		m15=2L^15
		a(0)=fix((f-vmin)*vscl+0.5)-m15
		return
	endif else if (ipack eq 3) then begin
		ntot=nx*ny
		nc=((ntot+3)/4)*2
		a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L+pos*(nc*4L+32))
		b=a(0)
		c1=string(b(0:15))
		c2=string(b(16:31))
		;print,c1,c2
		vmin=min(f)
		m15=2L^15
		vscl=(max(f)-vmin)/(2.^16-1.)
		if vscl eq 0. then vscl=1.
		vscl=1./vscl
		c1=string(vmin,format='(e16.8)')
		c2=string(vscl,format='(e16.8)')
		b(0:15)=byte(c1)
		b(16:31)=byte(c2)
		a(0)=b
		;print,string(b)
		a=assoc(1,intarr(nx,ny),irec(j*incr+i)+pos*(nc*4L+32))
		tmp=fix(long((f-vmin)*vscl+0.5)-m15)
		;print,min(tmp),max(tmp)
		a(0)=tmp
		return
	end
end
;***************************************************************************
;
;  Extract variables in a plane, unpacked, ascii or binary packed
;
;***************************************************************************
function fpln,j,i,pos
	common cmov,cvar,cpos,cfirst,clast,cstep,csize,cfile,ctmp
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
	common cplane,nc

	if n_params() eq 0 then begin
		print,'f=fpln(itime,ivar,n) [0 indexed]'
		return,0
	end
	if n_params() lt 3 then pos=cpos
	cpos=pos
	ntot=nx*ny*nz
	nxy=nx*ny

	; unpacked scr file
	if (incr gt 10) then begin   
		a=assoc(1,fltarr(nx,ny),irec(j*incr+i)+pos*nxy*4L)
		return,a(0)

	;  ascii packed sav file, ipack=1
	endif else if (ipack eq 1) then begin
		nc=((ntot+3)/4)*2
		a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L)

	 	b=a(0)
; 	        byteorder,b,/ntohl ;@

		c1=string(b(0:15))
		c2=string(b(16:31))
		vmin=float(c1)
		vscl=float(c2)
		;print,b,vmin,vscl
		vscl=1./vscl
		m15=2L^15
		; print,irec(j*incr+i)+pos*nxy*4L

;  		tmp=vmin+vscl*(m15+a(0))



		if cpos ge 0 then begin
			;print,'rd.offset=',irec(j*incr+i)+pos*nxy*4L,nx*ny
			a=assoc(1,intarr(nx,ny),irec(j*incr+i)+pos*nxy*2L)

 	 	b=a(0) ;@
	        byteorder,b,/ntohs ;@
  		

		;	return,vmin+vscl*(m15+a(0))
 	
		return,vmin+vscl*(m15+b) ;@

		end else begin
			tmp1=fltarr(nx,nz)
			for n=0,nz-1 do begin
			  a=assoc(1,intarr(nx),irec(j*incr+i)+(n*nxy-cpos*nx)*2L)

 	 	b=a(0) ;@
	        byteorder,b,/ntohs ;@


;			  tmp1(*,n)=vmin+vscl*(m15+a(0))
			  tmp1(*,n)=vmin+vscl*(m15+b)

			end
			return,tmp1		
		end
	endif else if (ipack eq 3) then begin
		ntot=nx*ny
		nc=((ntot+3)/4)*2
		m15=2L^15
		if cpos ge 0 then begin
		  a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L+pos*(nc*4L+32))
		  b=a(0)
		  c1=string(b(0:15))
		  c2=string(b(16:31))
		  vmin=float(c1)
		  vscl=float(c2)
		  vscl=1./vscl
		; print,irec(j*incr+i)+pos*nxy*4L
			a=assoc(1,intarr(nx,ny),irec(j*incr+i)+pos*(nc*4L+32))
	b=a(0) ;@
        byteorder,b,/ntohs ;@

;			return,vmin+vscl*(m15+a(0))
		return,vmin+vscl*(m15+b) ;@
	
        
 	end else begin
			tmp1=fltarr(nx,nz)
			for n=0,nz-1 do begin
			  a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L+n*(nc*4L+32))
			  b=a(0)
			  c1=string(b(0:15))
			  c2=string(b(16:31))
			  vmin=float(c1)
			  vscl=float(c2)
			  vscl=1./vscl
			  a=assoc(1,intarr(nx),irec(j*incr+i)+n*(nc*4L+32)-cpos*nx*2L)
			  tmp1(*,n)=vmin+vscl*(m15+a(0))
			end
			return,tmp1		
		end
	endif else if (ipack eq 2) then begin
		ntot=nx*ny*nz
		nc=(ntot+1)/2
		a=assoc(1,fltarr(2),irec(j*incr+i)+nc*4L)
		b=a(0)
		vmin=b(0)
		vscl=b(1)
		vscl=1./vscl
		;print,b,vmin,vscl
		m15=2L^15
		a=assoc(1,intarr(nx,ny),irec(j*incr+i)+pos*nxy*2L)
		return,vmin+vscl*(m15+a(0))
	end
end
;**************************************************************************
;
;  Float array record
;
;**************************************************************************
function frec,j
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	off=irec(j)
	len=(irec(j+1)-off)/4
	a=assoc(1,fltarr(len),off)
	return,a(0)
end
;***************************************************************************
;
;  Extract 3-d variables, unpacked, ascii or binary packed
;
;***************************************************************************
function var,j,i
;
;  extract var i, time j from file
;
common coff,mpar,mav,ipack,incr,ioff,time,voff,cfile,irec
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0
;
	if n_params() lt 2 then begin
		print,'tmp=var(itime,ivar)'
		return,0
	end
	ntot=nx*ny*nz
	if (incr eq 12) then begin   ;  unpacked hd scratch file
  		a=assoc(1,fltarr(nx,ny,nz),irec(j*incr+i+5))
  		tmp=a(0)
      endif else if (incr eq 18) then begin ; unpacked mhd scratch file
            a=assoc(1,fltarr(nx,ny,nz),irec(j*incr+i+8))
            tmp=a(0)
	endif else if ipack eq 0 then begin    ;  unpacked save file
  		a=assoc(1,fltarr(nx,ny,nz),irec(j*incr+i))
  		tmp=a(0)
	endif else if ipack eq 1 then begin    ;  ascii packed scr file
 	 	nc=((ntot+3)/4)*2
 	 	a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L)

 	 	b=a(0)
; 	        byteorder,b,/ntohl ;@

 	 	c1=string(b(0:15))
  		c2=string(b(16:31))
  		vmin=float(c1)
  		vscl=float(c2)
  		;print,b,vmin,vscl
  		vscl=1./vscl
  		m15=2L^15
  		a=assoc(1,intarr(nx,ny,nz),irec(j*incr+i))

;  		tmp=vmin+vscl*(m15+a(0))

 	 	b=a(0) ;@
	        byteorder,b,/ntohs ;@
  		tmp=vmin+vscl*(m15+b) ;@



	endif else if ipack eq 2 then begin    ;  binary packed scr file
		nc=(ntot+1)/2
		a=assoc(1,fltarr(2),irec(j*incr+i)+nc*4L)
		b=a(0)
		vmin=b(0)
		vscl=b(1)
		vscl=1./vscl
  		;print,b,vmin,vscl
		m15=2L^15
		a=assoc(1,intarr(nx,ny,nz),irec(j*incr+i))
		tmp=vmin+vscl*(m15+a(0))
	endif else if (ipack eq 3) then begin  
		tmp=fltarr(nx,ny,nz)
		for n=0,nz-1 do tmp(*,*,n)=fpln(j,i,n)
	end
	return,reform(tmp)
end
;***************************************************************************
pro wvar,f,j,i
;
;  write var i, time j to file
;
common coff,mpar,mav,ipack,incr,ioff,time,voff,cfile,irec
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0
;
	if n_params() eq 0 then begin
		print,'wvar,f,itime,ivar'
		return
	end
	close,1
	openu,1,cfile
	ntot=nx*ny*nz
	if (incr eq 12) then begin   ;  unpacked hd scratch file
  		a=assoc(1,fltarr(nx,ny,nz),irec(j*incr+i+5))
		a(0)=f
  		return
      endif else if (incr eq 18) then begin ; unpacked mhd scratch file
  		a=assoc(1,fltarr(nx,ny,nz),irec(j*incr+i+8))
		a(0)=f
  		return
	endif else if ipack eq 1 then begin    ;  ascii packed sav file
 	 	nc=((ntot+3)/4)*2
		vmin=min(f)
		vscl=max(f)-vmin
  		m15=2L^15
		cscl=(m15+m15-1)/vscl
 	 	a=assoc(1,bytarr(32),irec(j*incr+i)+nc*4L)
		c1=string(vmin,format='(e16.8)')
		c2=string(cscl,format='(e16.8)')
		b=bytarr(32)
		b(0:15)=byte(c1)
		b(16:31)=byte(c2)
  		print,b,vmin,cscl
		a(0)=b
  		a=assoc(1,intarr(nx,ny,nz),irec(j*incr+i))
  		a(0)=fix((f-vmin)*cscl+0.5)-m15
  		return
	endif else if ipack eq 2 then begin    ;  binary packed sav file
		b=fltarr(2)
		vmin=min(f)
		vscl=max(f)-vmin
		m15=2L^15
		cscl=(m15+m15-1)/vscl
		b(0)=vmin
		b(1)=cscl
		print,vmin,cscl,b
		nc=(ntot+1)/2
		a=assoc(1,fltarr(2),irec(j*incr+i)+nc*4L)
		a(0)=b
		a=assoc(1,intarr(nx,ny,nz),irec(j*incr+i))
		a(0)=fix((f-vmin)*cscl+0.5)-m15
		return
      endif else if ipack eq 3 then begin    ;  binary packed sav file
            s=size(f)
            for n=0,s(3)-1 do wpln,f(*,*,n),j,i,n
	end
end

function fint,j
	common coff,mpar,mav,ipack,incr,ioff,time,voff,cfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
	a=assoc(1,fltarr(nx,ny),irec(j*incr+incr-1)+4L*(mpar+mav*nz))
 	b=a(0) ;@
	byteorder,b,/ntohl ;@
	return,b
end
function flnrho,i
	return,var(i,0)
end
function fux,i
	return,var(i,1)
end
function fuy,i
	return,var(i,2)
end
function fuz,i
	return,var(i,3)
end
function fee,i
	return,var(i,4)
end
function fax,i
	return,var(i,5)
end
function fay,i
	return,var(i,6)
end
function faz,i
	return,var(i,7)
end
function fbx,i
        bx=yder(faz(i))
        bx=temporary(bx)-zder(fay(i))
        return,bx
end
function fby,i
        by=zder(fax(i))
        by=temporary(by)-xder(faz(i))
        return,by
end
function fbz,i
        bz=xder(fay(i))
        bz=temporary(bz)-yder(fax(i))
        return,bz
end

function ftt,i
       common coff,mpar,mav,ipack,incr,ioff,time,voff,cfile,irec
       if (incr le 10) then return,var(i,incr-2)    ; sav file
       if (incr eq 12) then return,var(i,incr-5-2)  ; scr file hd
       if (incr eq 18) then return,var(i,incr-8-2)  ; scr file mhd
end

;
;  Movie 
;
pro movie,var=var,pos=pos,first=first,last=last,step=step,size=size,vmin=vmin,vmax=vmax,sample=sample,file=file

	common cmov,cvar,cpos,cfirst,clast,cstep,csize,cfile,ctmp
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
	common cplane,nc

	flg=1
	print,!d.name

	if n_elements(var)    ne 0 then begin & cvar=var     & flg=0 & end
	if n_elements(pos)    ne 0 then begin & cpos=pos     & flg=0 & end
	if n_elements(first)  ne 0 then begin & cfirst=first & flg=0 & end
	if n_elements(last)   ne 0 then begin & clast=last   & flg=0 & end
	if n_elements(step)   ne 0 then begin & cstep=step   & flg=0 & end
	if n_elements(size)   ne 0 then begin & csize=size   & flg=0 & end
	if n_elements(file)   ne 0 then begin & cfile=file   & flg=0 & end

	if n_elements(sample) eq 0 then begin & sample=0             & end

	if n_elements(file)   eq 0 then begin
		if n_elements(cfile)  eq 0 then begin
			cfile="cmpsav"
			print,"file ?"
			read,cfile
			savefile,cfile
			flg=0
		end
	end else begin
		cfile=file
		savefile,cfile
		flg=0
	end
	if n_elements(cvar)   eq 0 then begin
		print,"variabel no (-1=int 0=lnrho 1=ux 2=uy 3=uz 4=ee 5=tt 8=u 9=m) ?"
		read,cvar
		cvar=long(cvar)
	end
	if n_elements(cfirst) eq 0 then begin
		print,"depth index ?"
		read,cpos
		cpos=long(cpos)
	end
	if n_elements(cfirst) eq 0 then begin
		print,"first snapshot ?"
		read,cfirst
		cfirst=long(cfirst)
	end
	if n_elements(clast) eq 0 then begin
		print,"last snapshot ?"
		read,clast
		clast=long(clast)
	end
	if n_elements(cstep) eq 0 then begin
		print,"step snapshot ?"
		read,cstep
		cstep=long(cstep)
	end
	if n_elements(csize) eq 0 then begin
		print,"size of frame (pixels) ?"
		read,csize
		csize=long(csize)
	end

	blow=long(csize/nx)
	sx=blow*nx
	sy=blow*ny
	nm=(clast-cfirst)/cstep+1
	tmax=nm-1

	if flg eq 0 then begin
		ntot=nx*ny*nz
		nc=(ntot+1)/2
		res=fltarr(nx,ny,nm>2)

		t=-1
		print,"reading ..."
		print,cfirst,clast,cstep
		for j=cfirst,clast,cstep do begin
			t=t+1
			if cvar eq 8 then begin
				res(*,*,t)=fpln(j,1)^2+fpln(j,2)^2+fpln(j,3)^2
				res(*,*,t)=sqrt(res(*,*,t))
			end else if cvar eq 9 then begin
				res(*,*,t)=fpln(j,1)^2+fpln(j,2)^2+fpln(j,3)^2
				b=1.67*1.38e-16*fpln(j,5)/(1.2*1.66e-24)*1e-12
				res(*,*,t)=sqrt(res(*,*,t)/b)
			end else if cvar ge 0 then begin
				b=fpln(j,cvar)
				res(*,*,t)=b
			end else begin
				intoff=4L*(mpar+mav*nz)
				a=assoc(1,fltarr(nx,ny),irec(j*incr+incr-1)+intoff)
				res(*,*,t)=a(0)
			end
			print,t,j
		endfor

		if n_elements(vmin) ne 0 then rmin=vmin else begin
			rmin=min(res)
			print,"vmin=",rmin
		end
		if n_elements(vmax) ne 0 then rmax=vmax else begin
			rmax=max(res)
			print,"vmax=",rmax
		end

		if sample eq 3 then begin
			tmpx=res
			res=fltarr(2*nx,2*ny,nm)
			for t=0,nm-1 do begin
				for j=0,ny-1 do res(*,j,t)=[tmpx(*,j,t),tmpx(*,j,t)]
				res(*,ny:ny+ny-1,t)=res(*,0:ny-1,t)
			endfor
			tmpx=0
			sx=sx*2
			sy=sy*2
		endif

		print,"rebinning ..."
		if !d.name eq "TEK" then begin
			if rmin eq rmax then begin
				ctmp=bytarr(sx,sy,nm)
				for i=0,nm-1 do $
					ctmp(*,*,i)=bytscl(res(*,*,i),top=255-64)+64
			end else $
				ctmp=bytscl(res,min=rmin,max=rmax,top=255-64)+64
		end else begin
			if rmin eq rmax then begin
				ctmp=bytarr(sx,sy,nm)
				for i=0,nm-1 do $
					ctmp(*,*,i)= $
					rebin(bytscl(res(*,*,i),top=250),sx,sy)
			end else $
				ctmp= $
				rebin(bytscl(res,min=rmin,max=rmax,top=250),sx,sy,nm)
		end
	end

	repeat begin
		for t=0,tmax do begin
			if !d.name eq "TEK" then begin
				am,ctmp(*,*,t)
			end else if sample eq 0 then begin
				tv,ctmp(*,*,t),0,sy
				tv,ctmp(*,*,t),sx,sy
				tv,ctmp(*,*,t)
				tv,ctmp(*,*,t),sy,0
			end else if sample eq 1 then begin
				tmpx=bytarr(2*sx,2*sy)
				tmpx(0:sx-1,0:sy-1)=ctmp(*,*,t)
				tmpx(sx:2*sx-1,0:sy-1)=ctmp(*,*,t)
				tmpx(0:sx-1,sy:2*sy-1)=ctmp(*,*,t)
				tmpx(sx:2*sx-1,sy:2*sy-1)=ctmp(*,*,t)
				tv,rebin(tmpx,4*sx,4*sy,/sam)
			end else if sample eq 2 then begin
				tmpx=bytarr(2*sx,2*sy)
				tmpx(0:sx-1,0:sy-1)=ctmp(*,*,t)
				tmpx(sx:2*sx-1,0:sy-1)=ctmp(*,*,t)
				tmpx(0:sx-1,sy:2*sy-1)=ctmp(*,*,t)
				tmpx(sx:2*sx-1,sy:2*sy-1)=ctmp(*,*,t)
				tv,tmpx
			end else if sample eq 3 then begin
				tv,ctmp(*,*,t)
			end
			xyouts,2*sx-60,2*sy-10,string(t),/device
			if !d.name eq "SUN" then empty
		endfor
		print,"again ?"
	end until get_kbrd(1) ne "y"
	print,cfile,cvar,cpos,format="(' file=',a,' var=',i1,' pos=',i2)"
	print,cfirst,clast,cstep,csize,format="(' first=',i2,' last=',i2,' step=',i2,'   size=',i3)"
end

function intens,j
;  intensity
	common cmov,cvar,cpos,cfirst,clast,cstep,csize,cfile,ctmp
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 

	ntot=nx*ny
	nc=((ntot+3)/4)*2
	a=assoc(1,fltarr(nx,ny),irec(j*incr)+4L*(mpar+mav*mz))
	return,a(0)
	end


function intpck,j
;  ascii packed intensity, ipack=3
	common cmov,cvar,cpos,cfirst,clast,cstep,csize,cfile,ctmp
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 

	ntot=nx*ny
	nc=((ntot+3)/4)*2
	a=assoc(1,bytarr(32),irec(j*incr)+nc*4L)
	b=a(0)
	;print,b
	c1=string(b(0:15))
	c2=string(b(16:31))
	vmin=float(c1)
	vscl=float(c2)
	;print,b,vmin,vscl
	vscl=1./vscl
	m15=2L^15
	a=assoc(1,intarr(nx,ny),irec(j*incr))
	return,vmin+vscl*(m15+a(0))
	end

pro file,name,writeaccess=writeaccess,debug=debug,quiet=quiet
	; load 3D sav or scr file, extract pointers

	common coff,mpar,mav,ipack,incr,ioff,time,voff,cfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0

	if n_elements(debug) eq 0 then on_error,1

; open file (automatically reads fixed record length)
	close,1
	if n_elements(writeaccess) ne 0 then openu,1,name else openr,1,name
	cfile=name

; get indices
	i=assoc(1,lonarr(17))
	inext=i(0,0)
	byteorder,inext,/ntohl

	lcfile=i(1,0)
        byteorder,lcfile,/ntohl

	ndata=i(2,0)
        byteorder,ndata,/ntohl

	nw=i(4,0)
        byteorder,nw,/ntohl

	;lrec=min([nw,4096])
	idata=i(5,0)
        byteorder,idata,/ntohl

	incr=i(7,0)
        byteorder,incr,/ntohl

	j=assoc(1,lonarr(10),4L*(1+lcfile+1000))
	lrec=j(0,0)
        byteorder,lrec,/ntohl

	if lrec eq 0 then lrec=4096
	;number of time steps in file
	ntmax=(idata-1)/incr+1
	if n_elements(quiet) eq 0 then $
	print,'incr  =',incr,'	ntmax =',ntmax,'	lrec =',lrec
	if n_elements(quiet) eq 0 then $
	print,'idata =',idata

	; get file pointers
	if inext eq 1 then inext = inext + 2*incr	; for 3d.scr wrap

	irecp=assoc(1,lonarr(inext),4L*(lcfile+1))
        tmp=irecp(0)
        byteorder,tmp,/ntohl
        
;	irec=4L*lrec*irecp(0)
	irec=4L*lrec*tmp

	; get data record, values ; cdata_88
	mpar=55			; number of scalar parameters in /cdata/
	mav=17			; number of (mz) arrays in /cdata/
	timeoffset=18		; offset of time ("t") param in /cdata/	
	compoffset=51		; offset of ipack

	data=assoc(1,lonarr(mpar),irec(idata-1))
	nx=data(0,0)
        byteorder,nx,/ntohl	

	ny=data(1,0)
        byteorder,ny,/ntohl	

	nz=data(2,0)
        byteorder,nz,/ntohl	

	if ndata gt mpar+nz*mav+nx*ny then begin
		mpar=58		; number of scalar parameters in /cdata/
		mav=17		; number of (mz) arrays in /cdata/
		timeoffset=18	; offset of time ("t") param in /cdata/	
		compoffset=51	; offset of ipack
		data=assoc(1,lonarr(mpar),irec(idata-1))
	endif
	if ndata gt mpar+nz*mav+nx*ny or incr eq 2 then begin
	;if ndata gt mpar+nz*mav+nx*ny then begin ; cdata_hd_891117
		mpar=73		; number of scalar parameters in /cdata/
		mav=22		; number of (mz) arrays in /cdata/
		timeoffset=17	; offset of time ("t") param in /cdata/	
		compoffset=31	; offset of ipack
		data=assoc(1,lonarr(mpar),irec(idata-1))
	endif
	if ndata gt mpar+nz*mav+nx*ny then begin
		mpar=75		; number of scalar parameters in /cdata/
		mav=22		; number of (mz) arrays in /cdata/
		timeoffset=17	; offset of time ("t") param in /cdata/	
		compoffset=31	; offset of ipack
		data=assoc(1,lonarr(mpar),irec(idata-1))
	endif
	if ndata gt mpar+nz*mav+nx*ny then begin
		mpar=75		; number of scalar parameters in /cdata/
		mav=35		; number of (mz) arrays in /cdata/
		timeoffset=17	; offset of time ("t") param in /cdata/	
		compoffset=31	; offset of ipack
		data=assoc(1,lonarr(mpar),irec(idata-1))
	endif
	if ndata gt mpar+nz*mav+nx*ny then begin ; cdata_hd_9105
		mpar=76		; number of scalar parameters in /cdata/
		mav=35		; number of (mz) arrays in /cdata/
		timeoffset=17	; offset of time ("t") param in /cdata/	
		compoffset=31	; offset of ipack
		data=assoc(1,lonarr(mpar),irec(idata-1))
	endif
	if ndata gt mpar+nz*mav+nx*ny then begin ; cdata_{9003,9007}
		mpar=76		; number of scalar parameters in /cdata/
		mav=39		; number of (mz) arrays in /cdata/
		timeoffset=17	; offset of time ("t") param in /cdata/	
		compoffset=31	; offset of ipack
		data=assoc(1,lonarr(mpar),irec(idata-1))
	endif
        if ndata gt mpar+nz*mav+nx*ny then begin ; cdata_mhd_9109
                mpar=81         ; number of scalar parameters in /cdata/
                mav=43          ; number of (mz) arrays in /cdata/
                timeoffset=17   ; offset of time ("t") param in /cdata/
                compoffset=31   ; offset of ipack
                data=assoc(1,lonarr(mpar),irec(idata-1))
        endif
        if ndata gt mpar+nz*mav+nx*ny then begin ; cdata_mhd
                mpar=83         ; number of scalar parameters in /cdata/
                mav=48          ; number of (mz) arrays in /cdata/
                timeoffset=17   ; offset of time ("t") param in /cdata/
                compoffset=31   ; offset of ipack
                data=assoc(1,lonarr(mpar),irec(idata-1))
        endif

	if compoffset gt 0 then begin
	  ipack=data(compoffset,0) 
          byteorder,ipack,/ntohl	
        end

	if compoffset le 0 then ipack=1
	
	if n_elements(quiet) eq 0 then $
	print,'nx    =',nx,'	ny   =',ny,  '	nz =',nz
	if n_elements(quiet) eq 0 then $
	print,'ndata =',ndata,'	mpar =',mpar,'		ipack =',ipack
	r1=assoc(1,fltarr(ndata),irec(idata-1-(ntmax-1)*incr))
	tbegin=r1(timeoffset,0)
        byteorder,tbegin,/ntohl	

	; vertical and horizontal scales
	rr=assoc(1,fltarr(ndata),irec(idata-1))
	z=rr(mpar:mpar-1+nz,0)
        byteorder,z,/ntohl	

	dx=rr(6,0)
        byteorder,dx,/ntohl	

	dy=rr(7,0)
        byteorder,dy,/ntohl	

	x=indgen(nx)*dx
	y=indgen(ny)*dy

	tend=rr(timeoffset,0)
        byteorder,tend,/ntohl	

	if n_elements(quiet) eq 0 then $
	print,'tbegin =',tbegin,'   tend =',tend

	;
	if ipack eq 2 then nw=nw/2  ; compressed records
	nbuf=(nw-1)/lrec+1
	nw=nbuf*lrec

	;
	if inext gt 1 then voff=irec(1)-irec(0) else voff=0
	ioff=lonarr(ntmax)
	time=fltarr(ntmax)
	for i=0,ntmax-1 do begin
		ioff(i)=irec(i*incr)
		data=assoc(1,fltarr(1),4L*timeoffset+irec(idata-1-(ntmax-1-i)*incr))
                tmp=data(0,0)

        byteorder,tmp,/ntohl	

;		time(i)=data(0,0)
		time(i)=tmp
	end

	;print,'times =',time
	;print,'ioff =',ioff

	ilnr=0
	ivx=1
	ivy=2
	ivz=3
	if incr eq 7 then begin
		ibx=-1
		iby=-1
		ibz=-1
		iee=4
		itt=5
	end else begin
		ibx=4
		iby=5
		ibz=6
		iee=7
		itt=8
	end

return
end


;************************************************************************
;  Extract variables of cdata
;
;************************************************************************

function cindex,iv
;
	common cdata_names, array

	s=size(iv)
	if s(0) eq 0 and s(1) ne 7 then return,iv

	if n_elements(array) le 1 then spawn,'cdata.help -1',array

	index=where(array eq iv, n)
	if n le 0 or index(0) le 0 then begin
		print,iv+' not in cdata:'
		print,array
		return,-1
	end
	return,index(0)
end

pro wdata,f,it,iv ;Write cdata at time it variable iv
;
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 

	if n_params() eq 0 then begin
		print,'wdata,f,isnap,iarray'
		return
	end
;
	iv=cindex(iv)
	a=assoc(1,fltarr(nz),ioff(it)+(incr-1)*voff+4L*(mpar+iv*nz))
	a(0)=f
	return
end


function cdataz,it,iv,help=help
;
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
	common cdata_names, array

	if n_params() eq 0 then print,'cdata(isnap,iarray)'
	if n_elements(help) gt 0 then begin
		if n_elements(array) le 1 then spawn,'cdata.help -1',array
		for i=0,n_elements(array)-1 do $
			print,i,array(i),format='(i3,2x,a)'
		if n_params() lt 2 then return,array
	end
	if n_params() eq 0 then return,0

	iv=cindex(iv)
	a=assoc(1,fltarr(nz),ioff(it)+(incr-1)*voff+4L*(mpar+iv*nz))
        b=a(0)
        byteorder,b,/ntohl

;	return,a(0)
	return,b
end

function cdatat,iv       ;cdataz variable iv  as fiction of time
;
	common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 

	if n_params() eq 0 then begin
		print,'tmp=cdatat(iarray)'
		return,0
	end
;
	a=fltarr(nz,ntmax,/nozero)
	iv=cindex(iv)
	for t=0,ntmax-1 do a(*,t)=cdataz(t,iv)
	return,a
end


;Subject: scatas.pro

;************************************************************************
;  Extract scalar all variables from cdata
;
;************************************************************************
 function cdatas,it
;function cdatas,it,iv
;
;  extract var iv, time it from file
;
common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
;
i=assoc(1,lonarr(10))
idata=i(5,0)
        byteorder,idata,/ntohl

;a=assoc(1,fltarr(75),irec(idata-1-(ntmax-1-it)*incr))
;a=assoc(1,fltarr(1),ioff(it)+(incr-1)*voff+4L*iv)

a=assoc(1,lonarr(mpar),irec(idata-1-(ntmax-1-it)*incr))
a=a(0)
byteorder,a,/ntohl
b=assoc(1,fltarr(mpar),irec(idata-1-(ntmax-1-it)*incr))
b=b(0)
byteorder,b,/ntohl

c=[a(0:2),b(3:4),a(5),b(6:9),a(10),b(11:21),a(22:25),b(26:30),$
  a(31:41),a(42:45),b(46:50),a(51:52),b(53:75)]

return,c

end
;
; nx(0),ny(1),nz(2),height(3),width(4),n0(5),dx(6),dy(7),dz(8),bzone(9),npot(10)
; teff(11),grav(12),abund(13),fnom(14),tmass0(15),tee0(16)
; t(17),dt(18),tstop(19),dtold(20),fdtc(21),ntime(22),nst1(23),nst3(24),nst5(25)
; dtsave(26),dtscr(27),dtprnt(28),dtsmpl(29),dtsmth(30),ipack(31),nsmpl(32)
; ip(33),idl(34),idm(35),idn(36),ldl(37),ldm(38),ldn(39),ledbug(40),idebug(41)
; dato(42),time(44)
; uur(46),uul(47),uut(48),uub(49),mu0(50)
; nmu(51),nph(52),taumax(53),dabmax(54)
; cnu1(55),cnu2(56),cnu3(57),cnu4(58),cnu5(59)
; eetop(60),eebot(61),rhotop(62),rhobot(63),pptop(64),ppbot(65),
; period(66),nwavel(67),amp0(68),bfm(69),bfke(70),bfe(71),
; spare1(72),spare2(73),spare3(74)

;************************************************************************
 function wdatas,f,it
;
;  write parameters, f, in cdata
;
common coff,mpar,mav,ipack,incr,ioff,time,voff,xfile,irec
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
;
i=assoc(1,lonarr(10))
idata=i(5,0)

a=assoc(1,lonarr(mpar),irec(idata-1-(ntmax-1-it)*incr))
;b=assoc(1,fltarr(mpar),irec(idata-1-(ntmax-1-it)*incr))

;c=[a(0:2),b(3:4),a(5),b(6:9),a(10),b(11:21),a(22:25),b(26:30),$
;  a(31:41),a(42:45),b(46:50),a(51:52),b(53:75)]


end
;
; nx(0),ny(1),nz(2),height(3),width(4),n0(5),dx(6),dy(7),dz(8),bzone(9),npot(10)
; teff(11),grav(12),abund(13),fnom(14),tmass0(15),tee0(16)
; t(17),dt(18),tstop(19),dtold(20),fdtc(21),ntime(22),nst1(23),nst3(24),nst5(25)
; dtsave(26),dtscr(27),dtprnt(28),dtsmpl(29),dtsmth(30),ipack(31),nsmpl(32)
; ip(33),idl(34),idm(35),idn(36),ldl(37),ldm(38),ldn(39),ledbug(40),idebug(41)
; dato(42),time(44)
; uur(46),uul(47),uut(48),uub(49),mu0(50)
; nmu(51),nph(52),taumax(53),dabmax(54)
; cnu1(55),cnu2(56),cnu3(57),cnu4(58),cnu5(59)