$Id: open.pro,v 1.4 2008/01/21 20:48:14 bob Exp $

FUNCTION findit,file,ext,debug=debug
  ;print,file,systime(1)
  if exists(file) then return,file
  file1=file
  f=findfile(file1,count=nf)
  if nf gt 0 then file1=f(0)
  ;print,file,nf,f,file1
  if keyword_set(debug) then print,'findit:',file1
  if fexists(file1) then return,file1

  file1=fileroot(file)+ext
  if keyword_set(debug) then print,'findit:',file1
  if fexists(file1) then return,file1

  file1='data/'+fileroot(file)+ext
  if keyword_set(debug) then print,'findit:',file1
  if fexists(file1) then return,file1

  file1='data/'+file
  if keyword_set(debug) then print,'findit:',file1
  if fexists(file1) then return,file1

  file1=file
  if keyword_set(debug) then print,'findit:',file1
  return,file1
END

PRO openp,p,file,_extra=extra
  COMMON cpart,npart
  default, file, 'particles.dat'
  open,p,file,dim=npart,nv=3,_extra=extra
END

PRO open,a,file,unit=unit,xy=xy,quiet=quiet,mhdsnaps=mhdsnaps,packed=packed,$
  update=update,boundaries=boundaries,debug=debug,ntime=ntime,times=times, $
  xmesh=xmesh,ymesh=ymesh,zmesh=zmesh,dim=dim,ds=ds,nvar=nvar1, $
  verbose=verbose,_extra=extra
@common_cdata
@common_cmesh
  COMMON cvert, cxvertical
  COMMON ctime, ctimes
  COMMON cpart, npart
  COMMON cmpi, mpi_size, mpi_ny ,mpi_nz, mpi_y, mpi_z

  default,file,'snapshot.dat'

  if fexists(file) then begin
    file1=file
  end else begin
    file1=findit(file,'.fe')
    if not fexists(file1) then file1=findit(file,'.dat')
    if not fexists(file1) then begin
      print,'file not found:',file
      return
    endif
  end

  if n_elements(ds) eq 0 then begin
    dx=1.
    dy=1.
    dz=1.
  end

  if n_elements(dim) eq 0 then begin

    fdim=findit(fileroot(file1,head=head)+'.dx','.dx',debug=debug)
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) and strpos(file1,'snapshot-') eq 0 then fdim='grid'+strmid(file1,8,6)+'.txt'
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) and strpos(file1,'scratch-') eq 0 then fdim='grid'+strmid(file1,7,6)+'.txt'
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) then fdim='grid.txt'
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) then fdim=head+'grid.txt'
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) then fdim='open.dx'
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) then fdim='*.dx'
    if keyword_set(debug) then print,fdim
    if not fexists(fdim) then fdim=head+'*.dx'
    if keyword_set(debug) then print,fdim

    if fexists(fdim) then begin
      dxfile=1

    end else begin
      dxfile=0
      fdim=findit(fileroot(file1,head=head)+'.dim','.dim',debug=debug)
      if keyword_set(debug) then print,fdim
      if not fexists(fdim) then fdim=head+'open.dim'
      if keyword_set(debug) then print,fdim
      if not exists(fdim) then fdim=head+'*.dim'
      if not fexists(fdim) then fdim='open.dim'
      if keyword_set(debug) then print,fdim
      if not fexists(fdim) then fdim='*.dim'
      if keyword_set(debug) then print,fdim
      if not fexists(fdim) then begin
    print,'file not found:',fdim
    return
      endif
    end

    default,quiet,0
    if quiet lt 2 then print,file1,':',fdim
    openr,l2,/get,fdim,/swap_if_big_endian
    tmp=lonarr(3)
    readf,l2,tmp
    nx=tmp(0)
    ny=tmp(1)
    nz=tmp(2)
    sz=3
    tmp=fltarr(3)
    readf,l2,tmp
    if n_elements(boundaries) ne 0 then dx=tmp(0)/(nx-boundaries-1) else dx=tmp(0)
    dy=tmp(1)
    dz=tmp(2)
    if not dxfile then begin
      dx=dx/nx
      dy=dy/ny
      dz=dz/nz
    end
    x=dx*(findgen(nx)-3)
    if n_elements(boundaries) ne 0 then x=dx*(findgen(nx)-3) else x=dx*findgen(nx)
    y=dy*findgen(ny)
    z=dz*findgen(nz)
    gamma=0.
    if eof(l2) eq 0 then readf,l2,gamma,cxvertical
    if keyword_set(verbose) then print,'gamma=',gamma,cxvertical,eof(l2)
    nvar=5L
    if eof(l2) eq 0 then readf,l2,nvar
    if keyword_set(verbose) then print,'nvar=',nvar,eof(l2)
    nvar1=nvar
    npart=0L
    if eof(l2) eq 0 then readf,l2,npart
    if keyword_set(verbose) then print,'npar=:',npart,eof(l2)
    mpi_size=0L & mpi_ny=0L & mpi_nz=0L & mpi_y=0L & mpi_z=0L
    if eof(l2) eq 0 then readf,l2,mpi_size,mpi_ny,mpi_nz,mpi_y,mpi_z
    close,l2
    free_lun,l2
    dim=[nx,ny,nz]

  end else begin
    sz=n_elements(dim)
    nx=long(dim(0))
    if sz ge 2 then ny  =dim(1) else ny=1L
    if sz ge 3 then nz  =dim(2) else nz=1L
    if sz ge 4 then nvar=dim(3) else nvar=1
    if n_elements(nvar1) ne 0 then nvar=nvar1
    gamma=5./3.
    cxvertical=1.
  end

  if n_elements(a) gt 0 then begin
    fs=fstat(a)
    if fs.open then free_lun,fs.unit
  end

  if n_elements(unit) gt 0 then begin
    fs=fstat(unit)
    if fs.open then free_lun,fs.unit
  end

  if keyword_set(update) then $
    openu,unit,/get,file1,error=error,/swap_if_big_endian,_extra=extra $
  else $
    openr,unit,/get,file1,error=error,/swap_if_big_endian,_extra=extra
  if (error ne 0) then begin
    print,'no permission'
    return
  end

  fs=fstat(unit)
  mhdsnaps=fs.size/(4L*nvar)/(nx*ny*nz)
  if not keyword_set(quiet) then begin
    xz=['z','x']
    print,'    FILE : ',file1
    print,'nx ny nz :',nx,ny,nz,nvar,format='(a,4i5)'
    print,'dx dy dz :',dx,dy,dz,format='(a,3f8.5)'
    print,'   gamma :',gamma,format='(a,f8.5)'
    print,'vertical : ',xz(cxvertical)

    if n_elements(nvar) gt 0 then begin
      if ( nvar eq 5 or nvar eq 6 ) then $
        print,'HD  snaps:',fs.size/(nvar*4L)/(nx*ny*nz) $
      else $
        print,'MHD snaps:',mhdsnaps
    end else begin
      print,'HD  snaps:',fs.size/20L/(nx*ny*nz)
      print,'MHD snaps:',mhdsnaps
    end
  end
  ntime=(mhdsnaps>1)

  if keyword_set(packed) then begin
    p=assoc(unit,lonarr(4))
    q=p(0)
    ndim=q(0) & nx=q(1) & ny=q(2) & nz=q(3)
    print,ndim,nx,ny,nz
    a={lun:unit, size:size(lonarr(nx,ny,nz,/nozero)), $
       recl:4L*nx*ny*nz, offset:4L*(4+4)}
  end else begin
    if n_elements(xy) eq 0 then begin
      if      sz eq 1 then a=assoc(unit,fltarr(nx,/nozero)) $
      else if sz eq 2 then a=assoc(unit,fltarr(nx,ny,/nozero)) $
      else if sz eq 3 then a=assoc(unit,fltarr(nx,ny,nz,/nozero)) $
      else if sz ge 4 then a=assoc(unit,fltarr(nx,ny,nz,nvar,/nozero))
    end else begin
      a=assoc(unit,fltarr(nx,ny,/nozero))
    end
  end

  timfile = fileroot(file)+'.tim'
  if not fexists(timfile) then timfile='time.dat'
  if fexists(timfile) then begin
    openr,/get,u,timfile,/swap_if_big_endian,_extra=extra
    fs=fstat(u)
    print,'unit=',u
    if fs.size ge ntime*3*4 then begin
      if fs.size ge ntime*4*4 then d=assoc(u,fltarr(4)) else d=assoc(u,fltarr(3))
      times=fltarr(ntime)
      for i=0,ntime-1 do times[i]=d[0,i]
      ctimes=times
      if not keyword_set(quiet) then print,'time['+str(ntime-1)+'] = '+str(times[ntime-1])
    end
    free_lun,u
  end

  meshfile = fileroot(file1)+'.mesh'
  if fexists(meshfile) then read_mesh, meshfile, $
     xmesh=xmesh, ymesh=ymesh, zmesh=zmesh, _extra=extra
  meshfile = fileroot(file1)+'.mshtxt'
  if fexists(meshfile) then begin
    openr,/get,munit,meshfile,/swap_if_big_endian
    mx=0L & my=0L & mz=0L
    readf,munit,mx & xmesh=fltarr(mx) & readf,munit,xmesh
    readf,munit,my & ymesh=fltarr(my) & readf,munit,ymesh
    readf,munit,mz & zmesh=fltarr(mz) & readf,munit,zmesh
    free_lun,munit
  end
  meshfile = 'mesh.dat'
  if not fexists(meshfile) then meshfile = fileroot(file1)+'.msh'
  if fexists(meshfile) then begin
    readmesh,meshfile,xmesh=xmesh,ymesh=ymesh,zmesh=zmesh,_extra=extra
    x=xmesh & y=ymesh & z=zmesh
  end
END