#!/bin/csh

#-----------------------------------------------------------------------
#   This is commented, self executing f77 example program for reading
#   the OpenGGCM 3d data.  Make this file executable and check/modify the
#   shell commands (compiler instructions, file names) at the end.  Then
#   execute this file to see the example working.
#   Jimmy Raeder, IGPP/UCLA.  Adapted for the WEB July 22, 1998
#   OpenGGCM update  2006-11  JR
#-----------------------------------------------------------------------


set echo

cat >! tmp.f <<'EOF'
c
c.. this is an example program to read the MHD simulation data
c   and interpolate the to an arbitrary point within the simulation
c   box.
c      
c.. What you really need to know:
c
c  0.) The coordinate system is 'inverted GSE', i.e., -X_gse,-Y_gse,Z_gse.
c  1.) data are stored in a `compressed ascii' format, and thus are basically
c      machine independent (except for vintage mainframes). The 8'th bit
c      is used, thus `ftp' transfers should be done in `binary' mode since
c      some ftp clients strip the 8'th bit. However, there are no control
c      characters ( ascii 0 -- 31 ) in the files, except for line feed.
c  2.) open files as 'formatted'.
c  3.) read a 1d field with:
c      call getf11(iunit,a,n,l1,l2,it)
c      iunit (integer, input):            fortran file unit
c      a (real, output):                  field
c      n (integer, output):               size of a
c      l1 (character*80, input/output):   field identifier
c      l2 (character*80, output):         character tag
c      it (integer, output):              integer tag
c  4.) note:
c      a) you dont need to know the filed size, it is returned in `n',
c         however, `a' must be dimensioned sufficiently large.
c      b) `l1' and `l2' MUST be declared `character*80' in the calling 
c         program.
c      c) if `l1' is set to 'any', then the next field will be read.
c      d) if `l1' is set to a known identifier (like 'gridx'), the routine
c         will skip if necessary to that field and read it.
c      e) on `EOF' or any read error `n' returns less than zero.
c      f) `l2' usually contains some additional information, like the UT time.
c      g) `it' usually contains a value of seconds since the start of the run.
c  5.) read 2d and 3d fields with:
c      call getf21(iunit,a,nx,ny,l1,l2,it) and
c      call getf31(iunit,a,nx,ny,nz,l1,l2,it) , respectively.
c      Everything works the same way as with 1d fields, except that the multiple
c      field dimensions (nx,ny) or (nx,ny,nz) are returned.
c  6.) The returned values in `a' are of limited precision, about 3.e-4 in the
c      worst case.  That is good enough for virtually any purpose, because 
c      interpolation errors are generally larger.  `a' also contains no true
c      zeroes.  Zero values are returned as the smallest possible positive
c      value.
c  7.) The routines have been used on Intel-Linux boxes (pgf77 or fort77), SunOS,
c      Solaris, IBM RS6000 (AIX), CRAY-C90, DEC/Ultrix, CRAY-T3E, INTEL-Paragon,
c      and should run without problems on any UNIX workstation.
c      And, yes, GNU g77 works too.
c      I dont know about Microsoft, Apple, and what else there might exist.
c  8.) Corresponding write routines (datf11, datf21, datf31) are included.  Their
c      use should be straightforward.
c      
c.. What you may want to know:
c
c 1.)  Each field has a 5 line header that stores some pertinent
c      information (magic number, dimensions, etc.).
c 2.)  The field of values `a' is cropped to a dynamic range (ratio of
c      maximum to minimum abs value) of no more than 1.0e16, then is chopped 
c      into chunks of 64.  For each chunk the log of the abs value is taken and
c      linearly discretized in about 4400 steps (that number, with the dynamic range
c      determines the precision).  The resulting integer is factored into 2
c      integers and combined with the sign, resulting in 2 integers in the range
c      0-94. These are converted into ascii characters (33-126), which are then
c      run length encoded and written to the file (2 lines per 64 chunk).
c      The 8'th bit is used to indicate run length counts, thus all stored
c      characters lie in the range (33-126,172-254). In addition, a checksum is 
c      stored for each row and compared upun reading.  Thus, damaged files are 
c      easily detected.  Compression is usually
c      pretty good, `gzipping', i.e., nearly optimal compression of the files,
c      usually yields only 25% further compression.  In the worst case, 2 bytes
c      per value are used, typical is 1-1.5 bytes per value, and constant or nearly
c      constant fields compress to 0.1 bytes per value or better.
c 3.)  Using the software with C may be possible in at least 3 ways:
c      a) Compiling the routines with a f77 compiler, compiling the C routines
c      (including main()) with a C compiler and linking with the C compiler.
c      In this case, the proper calling sequence for the fortran routines must be
c      used, i.e., parameters passed by pointers, and character strings pass 2
c      pointers, possibly appending an underscore character to the fortran routine
c      name.  Also, the fortran libraries (something like libU77.a and libF77.a
c      must be provided in the link step.
c      b) Using f2c (from ftp.netlib.org) to convert the f77 code to C code and
c      compiling and linking everything with a C compiler.  In this case the
c      proper f2c libraries must be used because the f77 code uses IO and math
c      functions.
c      c) Rewrite the f77 routines in C.
c      I have not tried either and I would appreciate to hear about the results
c      if anyone tries out the above steps.
c 4.)  Things to do:
c      a) write a native C version
c      b) use different encoding (non-log) for fields with limited range, like
c      B-components and V-components.  That should improve compression and allow
c      for better precision.
c
c--------------------------------------------------------
      program mhdread
c--------------------------------------------------------
c
c   example program to read OpenGGCM 3d fields
c
      common /a1/bx(1500000),by(1500000),bz(1500000),
     *vx(1500000),vy(1500000),vz(1500000),rr(1500000),
     *pp(1500000)
      common /a2/gx(200),gy(200),gz(200)
      character*80 l1,l2,fgrid,f3d
c
      call getarg(1,fgrid)
      open(10,file=fgrid,status='old')
c
c... read grid parameters
c    gridpoint (ix,iy,iz) is at (gx(ix),gy(iy),gz(iz)) in RE
c
      l1='gridx'
      call getf11(10,gx,nnx,l1,l2,it)
      nx=nnx
      write(0,*)' got grid x ',l1(1:jlen(l1,80)),' ',nx
      write(0,*)' from ',gx(1),' to ',gx(nx)
      l1='gridy'
      call getf11(10,gy,nny,l1,l2,it)
      ny=nny
      write(0,*)' got grid y ',l1(1:jlen(l1,80)),' ',ny
      write(0,*)' from ',gy(1),' to ',gy(ny)
      l1='gridz'
      call getf11(10,gz,nnz,l1,l2,it)
      nz=nnz
      write(0,*)' got grid z ',l1(1:jlen(l1,80)),' ',nz
      write(0,*)' from ',gz(1),' to ',gz(nz)
      close(10)
c
c..... read fields
c      Bx,By,Bz in nT, rr (density) in cm**-3
c
      call getarg(2,f3d)
      open(10,file=f3d,status='old')
c
      l1='bx'
      call getf31(10,bx,nnx,nny,nnz,l1,l2,it)
      write(0,*)' got ',l1(1:jlen(l1,80)),' ',nnx,nny,nnz
c
      l1='by'
      call getf31(10,by,nnx,nny,nnz,l1,l2,it)
      write(0,*)' got ',l1(1:jlen(l1,80)),' ',nnx,nny,nnz
c
      l1='bz'
      call getf31(10,bz,nnx,nny,nnz,l1,l2,it)
      write(0,*)' got ',l1(1:jlen(l1,80)),' ',nnx,nny,nnz
c
      l1='rr'
      call getf31(10,rr,nnx,nny,nnz,l1,l2,it)
      write(0,*)' got ',l1(1:jlen(l1,80)),' ',nnx,nny,nnz
c
c... read next field, which should be 'pp_sm', i.e., the pressure
c
      l1='any'
      call getf31(10,rr,nnx,nny,nnz,l1,l2,it)
      write(0,*)' got ',l1(1:jlen(l1,80)),' ',nnx,nny,nnz
      close(10)
c
c..... ok, everything is in memory, do one interpolation for example
c      "iout" is returned 0 if (x,y,z) is in the box, 1 otherwise
c      do not use anything with
c      x**2+y**2+z**2 < 5.0**2, it may be trash
c
      y=6.2
      z=7.1
      do 100 x=-15.0,-5.0,0.1
      call ipol3a(bx,nx,ny,nz,bx_xyz,gx,gy,gz,x,y,z,iout)
      call ipol3a(by,nx,ny,nz,by_xyz,gx,gy,gz,x,y,z,iout)
      call ipol3a(bz,nx,ny,nz,bz_xyz,gx,gy,gz,x,y,z,iout)
      call ipol3a(rr,nx,ny,nz,rr_xyz,gx,gy,gz,x,y,z,iout)
      write(0,'(5(1x,g12.5))')x,rr_xyz,bx_xyz,by_xyz,bz_xyz
100   continue
      stop
      end
c.--------------------------------------------------------------
      subroutine ipol3a(a,nx,ny,nz,b,gx,gy,gz,x,y,z,iout)
c.--------------------------------------------------------------
c
c     routine for linear interpolation of a scalar value on a
c     MHD grid.  Somewhat optimized to avoid index searches between
c     calls with similar x,y,z input.
c     Jimmy Raeder, IGPP/UCLA.  Modified for the WEB, July 22, 1997.
c
c     a (real, input):                    grid function
c     nx,ny,nz (integer, input):          dimensions of a
c     b (real, output):                   interpolated value
c     gx,gy,gz (real 1d arrays, input):   grid coordinates
c     x,y,z (real, input):                coordinates to which to interpolate
c     iout (integer, output):             error flag, ne zero if x,y,z outside grid
c
c
c
      real a(nx,ny,nz),gx(nx),gy(ny),gz(nz)
      data ixl,iyl,izl /1,1,1/ 
c
c.... optimize search: most likely index is unchanged or in neighbor cell
c
      iout=0
      ix=ixl
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94979 
      ix=min0(ixl+1,nx-1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94979 
      ix=max0(ixl-1,1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94979 
      do 94971 i=1,nx-1
      if(x.ge.gx(i).and.x.le.gx(i+1)) then
      ix=i
      goto 94979 
      endif
94971 continue
      iout=1
      return
94979 continue
      ixl=ix
      iy=iyl
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94950 
      iy=min0(iyl+1,ny-1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94950 
      iy=max0(iyl-1,1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94950 
      do 94942 i=1,ny-1
      if(y.ge.gy(i).and.y.le.gy(i+1)) then
      iy=i
      goto 94950 
      endif
94942 continue
      iout=1
      return
94950 continue
      iyl=iy
      iz=izl
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94921 
      iz=min0(izl+1,nz-1)
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94921 
      iz=max0(izl-1,1)
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94921 
      do 94913 i=1,nz-1
      if(z.ge.gz(i).and.z.le.gz(i+1)) then
      iz=i
      goto 94921 
      endif
94913 continue
      iout=1
      return
94921 continue
      izl=iz
      dxi=(x-gx(ix))/(gx(ix+1)-gx(ix))
      dyi=(y-gy(iy))/(gy(iy+1)-gy(iy))
      dzi=(z-gz(iz))/(gz(iz+1)-gz(iz))
      x00=a(ix,iy,iz)+dxi*(a(ix+1,iy,iz)-a(ix,iy,iz))
      x10=a(ix,iy+1,iz)+dxi*(a(ix+1,iy+1,iz)-a(ix,iy+1,iz))
      x01=a(ix,iy,iz+1)+dxi*(a(ix+1,iy,iz+1)-a(ix,iy,iz+1))
      x11=a(ix,iy+1,iz+1)+dxi*(a(ix+1,iy+1,iz+1)-a(ix,iy+1,iz+1))
      y0=x00+dyi*(x10-x00)
      y1=x01+dyi*(x11-x01)
      b=y0+dzi*(y1-y0)
      return
      end
c-----------------------------------------------------------
      subroutine skipii(iu,istep,isdim)
c-----------------------------------------------------------
      character*80 l1,l2
      character*10 f1,f2
      if(isdim.eq.1)f2='FIELD-1D-1'
      if(isdim.eq.2)f2='FIELD-2D-1'
      if(isdim.eq.3)f2='FIELD-3D-1'
100   continue
      read(iu,'(a)',end=190)f1
      if(f1.ne.f2) goto 100
      read(iu,'(a)',end=190)l1
      read(iu,'(a)',end=190)l2
      if(isdim.eq.1)read(iu,*)it,nx
      if(isdim.eq.2)read(iu,*)it,nx,ny
      if(isdim.eq.3)read(iu,*)it,nx,ny,nz
      if(it.lt.istep) goto 100
190   continue
      backspace(iu)
      backspace(iu)
      backspace(iu)
      backspace(iu)
      backspace(iu)
      return
      end
c-----------------------------------------------------------
      subroutine getf11(iu,a1,nx,l1,l2,it)
c-----------------------------------------------------------
      real a1(*)
      character*80 rec
      character*80 l1,l2
      m=len(l1)
      call end0(l1,m)
      isany=0
      if(l1(1:m).eq.'any') isany=1
100   continue
      read(iu,1000,end=190) rec
      if(rec(1:10).eq.'FIELD-1D-1') then
      read(iu,1000,end=190) rec
      if((isany.eq.0).and.l1(1:m).ne.rec(1:m)) goto 100
      l1=rec
      read(iu,1000,end=190) l2
      read(iu,*,end=190)it,nx
      call rdn2(iu,a1,nx,rec,it,rid)
      return
      endif
      goto 100
190   continue
      nx=-1
1000  format(a)
      return
      end
c-----------------------------------------------------------
      subroutine getf21(iu,a1,nx,ny,l1,l2,it)
c-----------------------------------------------------------
      real a1(*)
      character*80 rec
      character*80 l1,l2
      m=len(l1)
      call end0(l1,m)
      isany=0
      if(l1(1:m).eq.'any') isany=1
100   continue
      read(iu,1000,end=190) rec
      if(rec(1:10).eq.'FIELD-2D-1') then
      read(iu,1000,end=190) rec
      if((isany.eq.0).and.l1(1:m).ne.rec(1:m)) goto 100
      l1=rec
      read(iu,1000,end=190) l2
      read(iu,*,end=190)it,nx,ny
      call rdn2(iu,a1,nn,rec,it,rid)
      if(nn.lt.0)nx=nn
      return
      endif
      goto 100
190   continue
      nx=-1
1000  format(a)
      return
      end
c-----------------------------------------------------------
      subroutine end0(r,m)
c-----------------------------------------------------------
      character*80 r
      n=min0(m,len(r))
      do 100 i=1,n
      m=i-1
      if(r(i:i).eq.' ') return
100   continue
      m=n
      return
      end
c-----------------------------------------------------------
      subroutine datf11(iu,a1,nx,l1,l2,it)
c-----------------------------------------------------------
      character*80 l1,l2
      real a1(*)
      write(iu,'(a)')'FIELD-1D-1'
      write(iu,'(a)') l1
      write(iu,'(a)') l2
      write(iu,*) it,nx
      call wrn2(iu,a1,nx,'FUNC-1-1',it,float(it))
      return
      end
c-----------------------------------------------------------
      subroutine datf21(iu,a1,nx,ny,l1,l2,it)
c-----------------------------------------------------------
      character*80 l1,l2
      real a1(*)
      write(iu,'(a)')'FIELD-2D-1'
      write(iu,'(a)') l1
      write(iu,'(a)') l2
      write(iu,*)it,nx,ny
      call wrn2(iu,a1,nx*ny,'FUNC-2-1',it,float(ny))
      return
      end
c-----------------------------------------------------------
      subroutine datf31(iu,a1,nx,ny,nz,l1,l2,it)
c-----------------------------------------------------------
      character*80 l1,l2
      real a1(*)
      write(iu,'(a)')'FIELD-3D-1'
      write(iu,'(a)') l1
      write(iu,'(a)') l2
      write(iu,*)it,nx,ny,nz
      call wrn2(iu,a1,nx*ny*nz,'FUNC-3-1',it,float(ny))
      return
      end
c-----------------------------------------------------------
      subroutine getf31(iu,a1,nx,ny,nz,l1,l2,it)
c-----------------------------------------------------------
      real a1(*)
      character*80 rec
      character*80 l1,l2
      m=len(l1)
      call end0(l1,m)
      isany=0
      if(l1(1:m).eq.'any') isany=1
100   continue
      read(iu,1000,end=190) rec
      if(rec(1:10).eq.'FIELD-3D-1') then
      read(iu,1000,end=190) rec
      if((isany.eq.0).and.l1(1:m).ne.rec(1:m)) goto 100
      l1=rec
      read(iu,1000,end=190) l2
      read(iu,*,end=190) it,nx,ny,nz
      call rdn2(iu,a1,nn,rec,it,rid)
      if(nn.lt.0)nx=nn
      return
      endif
      goto 100
190   continue
      nx=-1
1000  format(a)
      return
      end
c.---------------------------------------
      subroutine wrn2(iu,a,n,cid,it,rid)
c.---------------------------------------
      real a(n)
      character*8 cid
      real a1(0:63)
      integer i1(0:63),i2(0:63),i3(0:63)
      zmin=1.e33
      zmax=1.e-33
      do 100 i=1,n
      b=abs(a(i))
      zmin=amin1(zmin,b)
      zmax=amax1(zmax,b)
100   continue
1000  format(a,i8,3e14.7,i8,a)
      zmin=amax1(1.e-33,zmin)
      zmax=amin1( 1.e33,zmax)
      z2=alog(zmax)
      z1=-76.
      if(zmin.ne.0.)z1=alog(zmin)
      z1=amax1(z1,z2-37.)
      if(abs(z2-z1).le.1.e-5) then
      z1=z1-1.0
      z2=z2+1.0
      endif
      z0=exp(z1)
      dz=float(4410)/(z2-z1)
      write(iu,1000)'WRN2',n,z1,z2,rid,it,cid
      do 200 k=1,n,64
      nk=min0(63,n-k)
      do 210 i=0,nk
      a1(i)=amin1(1.e33,abs(a(i+k)))
      a1(i)=amax1(z0,a1(i))
      a1(i)=dz*(alog(a1(i))-z1)+0.5
      i3(i)=a1(i)
      i3(i)=max0(0,i3(i))
      i3(i)=min0(4414,i3(i))
      i1(i)=i3(i)/94
      i2(i)=i3(i)-94*i1(i)
      if(a(i+k).lt.0.)i1(i)=i1(i)+47
      i1(i)=i1(i)+33
      i2(i)=i2(i)+33
210   continue
      call wrnenc(iu,i1,nk)
      call wrnenc(iu,i2,nk)
200   continue
      return
      end
c.---------------------------------------
      subroutine wrnenc(iu,i1,n)
c.---------------------------------------
      integer i1(0:63)
      integer i2(-1:63)
      character c*72
      ick=i1(0)
      ir=i1(0)
      ic=1
      k=-1
      do 100 i=1,n
      ick=ick+i1(i)
      if(i1(i).eq.ir) then
      ic=ic+1
      else
      if(ic.eq.1) then
      k=k+1
      i2(k)=ir
      else
      k=k+1
      i2(k)=ic+170
      k=k+1
      i2(k)=ir
      endif
      ic=1
      ir=i1(i)
      endif
100   continue
      if(ic.eq.1) then
      k=k+1
      i2(k)=ir
      else
      k=k+1
      i2(k)=ic+170
      k=k+1
      i2(k)=ir
      endif
      i2(-1)=33+mod(ick,92)
      j=0
      do 200 i=-1,k
      j=j+1
      c(j:j)=char(i2(i))
200   continue
      write(iu,'(a)')c(1:j)
      return
      end
c.---------------------------------------
      subroutine wrndec(iu,i1,n)
c.---------------------------------------
      integer i1(0:63)
      common /wrdqq9/i2(-2:67)
      character c*72
      c=' '
      read(iu,'(a)',end=900,err=900) c
      i=-2
      j=0
100   continue
      j=j+1
      if(j.gt.67) then
      write(0,*)'wrndec: cant find end of encoded record'
      n=-5
      return
      endif
      ir=ichar(c(j:j))
      if(ir.eq.32) goto 190
      if(ir.gt.127) then
      ic=ir-170
      j=j+1
      ir=ichar(c(j:j))
      do 110 l=1,ic
      i=i+1
      i2(i)=ir
110   continue
      else
      i=i+1
      i2(i)=ir
      endif
      goto 100
190   continue
      ick=0
      n=i
      if(n.gt.63) then
      write(0,*)'wrndec: n gt 63, n=',n
      n=-5
      write(0,'(a)')'rec:'
      write(0,'(a)')c
      return
      endif
      do 200 i=0,n
      i1(i)=i2(i)
      ick=ick+i1(i)
200   continue
      if(i2(-1).ne.33+mod(ick,92)) then
      write(0,*)'wrndec: checksum error '
      write(0,'(a,a)')'rec:',c
      write(0,*)i2(-1),33+mod(ick,92)
      n=-5
      return
      endif
      return
900   continue
      write(0,*)' wrndec eof/err '
      n=-5
      return
      end
c.---------------------------------------
      subroutine rdn2(iu,a,n,cid,it,rid)
c.---------------------------------------
c character decoding
      real a(*)
      character*8 cid
      character*4 did
      real a1(0:63)
      integer i1(0:63),i2(0:63),i3(0:63)
100   continue
      read(iu,1000,err=100,end=900)did,n,zmin,zmax,rid,it,cid
1000  format(a,i8,3e14.7,i8,a)
      if(did.ne.'WRN2') goto 100
      if(zmin.eq.zmax) then
      do 200 i=1,n
      a(i)=zmin
200   continue
      return
      endif
      l=0
      do 300 k=1,n,64
      nk=min0(63,n-k)
      call wrndec(iu,i1,nn)
      if(nn.ne.nk) then
      write(0,*)'rdn2: nn .ne. nk ',nn,nk,n,k
      n=-2
      return
      endif
      call wrndec(iu,i2,nn)
      if(nn.ne.nk) then
      write(0,*)'rdn2: nn .ne. nk ',nn,nk
      n=-2
      return
      endif
      dzi=(zmax-zmin)/float(4410)
      do 330 i=0,nk
      i1(i)=i1(i)-33
      i2(i)=i2(i)-33
      sig=1.
      if(i1(i).ge.47) then
      sig=-1.
      i1(i)=i1(i)-47
      endif
      i3(i)=i2(i)+94*i1(i)
      a1(i)=i3(i)
      a1(i)=dzi*a1(i)+zmin
      a1(i)=sig*exp(a1(i))
      l=l+1
      a(l)=a1(i)
330   continue
300   continue
      return
900   continue
      n=-1
      return
      end
c-------------------------------------
      integer function jlen(c,m)
c-------------------------------------
      character*80 c
      k=1
      do 100 i=1,m
      if(c(i:i).ne.' ')k=i
100   continue
      jlen=k
      return
      end
'EOF'

g77 tmp.f