#!/bin/csh

#
#..... this is a self executing shell script
#      that compiles and executes a f77 program
#      to access global MHD ionosphere data
#
#      To run make sure 'ioread.com' is executable (chmod +x ioread.com)
#      and edit at the end of this file for your favorite
#      f77 compiler and for a correct input file
#      the type 'ioread.com'
#
#      Joachim Raeder, IGPP/UCLA  7/14/98
#
#  

set echo

cat >! tmp.f <<'eeeeee'
c----------------------------------------------
      program main
c----------------------------------------------
c
c.... no of gridpoints in azimuth
      parameter (nnp=121)
c
c.... no of gridpoints in co-latitude
      parameter (nnt=361)
c
      character*80 infile,l1,l2
      common /aaa1/a(361*121)
c
c..... get input file from command line and open file
c
      call xparse('-fin',1,1,'required',infile,idum,rdum)
      open(10,file=infile,status='old')
c
c..... what do we want?
c..... choices:
c                prec_e_e0_1:  mean particle energy of diffuse electron
c                              precipitation (eV)
c                prec_e_fe_1:  energy flux of diffuse electron
c                              precipitation (W/m**2)
c                rrio:         mapped particle density (m**-3)
c                prec_e_e0_2:  mean particle energy of discrete electron
c                              precipitation (eV)
c                prec_e_fe_2:  energy flux of discrete electron
c                              precipitation (W/m**2)
c                pacurr:       field aligned current density, A/m**2
c                              downward positive
c                pot:          potential (V)
c                delphi:       downward accelerating potential drop (V)
c                ppio:         mapped MHD pressure (Pa)
c                sigp:         Pedersen conducatance (S)
c                sigh:         Hall conducatance (S)
c                delbt:        ground magnetic perturbation, 
c                              nortward component (nT)
c                delbr:        ditto, vertical component
c                delbp:        ditto, east-west component
c                xjh:          Joule heating, not unnormalized yet
c                any:          getf21 will read the next record and
c                              return in 'l1' what it got
c
c     note that the sequnce in which things are read may be important.
c     you may need a rewind statement.
c
c
c...... getf21 will read 'pacurr' from unit 10 and return:
c       a:      the array
c       np,nt:  the array dimensions, or something <=0 on eof or error
c       l2:     80 char string with some info like UT time
c       istep:  which is in seconds the time since the start of the simulation
c               should be the same number as encoded in the filename
c
      l1='any'   
      call getf21(10,a,np,nt,l1,l2,istep)
      write(0,*)'got: ',np,nt,istep,' ',l1(1:jlen(l1,80))
      write(0,*)'got: ',l2(1:jlen(l2,80))
      if(np.le.0) then
      write(0,*)'reading error or end-of-file ..... stopping.'
      stop
      endif
c
      l1='any'   
      call getf21(10,a,np,nt,l1,l2,istep)
      write(0,*)'got: ',np,nt,istep,' ',l1(1:jlen(l1,80))
      write(0,*)'got: ',l2(1:jlen(l2,80))
      if(np.le.0) then
      write(0,*)'reading error or end-of-file ..... stopping.'
      stop
      endif
c
c.... the way things are stored in a(np,nt):
c
c     1. coordinate is longitude in SM coordinates from 
c     -180 degrees (midnight, 24MLT) going eastward (increasing MLT)
c      to 180 degrees
c
c     2. coordinate is magnetic co-latitude starting from the
c     magnetic north pole
c
      phi1=-180.0
      phi2=180.0
      the1=0.0
      the2=180.0
      dphi=(phi2-phi1)/float(np-1)
      dthe=(the2-the1)/float(nt-1)
c
c..... interpolate to p,t coordinates
c
      do 100 p=-180.0,180.0,15.0
      do 100 t=25.0,25.0,0.5
      call inp2d1(a,np,nt,phi1,dphi,the1,dthe,p,t,w,wx,wy)
c
c....  inp2d1 returns the interpolated value in 'w'.
c      wx and wy are the dreivatives
c
      xmlt=(p+180.0)/15.0
      write(0,'(a,a,a,f7.2,1x,f7.2,a,g11.5)')
     *'got interpolated value of ',
     * l1(1:jlen(l1,80)),' at ',xmlt,t,' ==> ',w
c
100   continue
c
      stop
      end
c
c...... this is all library stuff
c
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 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
c      write(0,*)' zmin,zmax ',zmin,zmax
      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
cdir$ shortloop
      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---------------------------------------
      subroutine wrndec(iu,i1,n)
c---------------------------------------
      integer i1(0:63)
      common /wrdqq9/i2(-2:67)
      character c*72
      nsa=n
      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,*)'error:wrndec: cant find end of encoded record'
      goto 910
      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,*)'error:wrndec: n gt 63, n=',n
c      n=-5
c      write(0,'(a)')'rec:'
c      write(0,'(a)')c
      goto 910
      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,*)'error:wrndec: checksum error '
c      write(0,'(a,a)')'rec:',c
c      write(0,*)i2(-1),33+mod(ick,92)
c      n=-5
      goto 910
      endif
      return
910   continue
      do 911 i=0,63
      i1(i)=34
911   continue
      n=nsa
      return
900   continue
      write(0,*)' wrndec eof/err '
      n=-5
      return
      end
c-------------------------------------
      integer function jlen(c,m)
c-------------------------------------
      character*80 c
      k=1
      do 100 i=1,m
      if(c(i:i).ne.' '.and.c(i:i).ne.char(0))k=i
100   continue
      jlen=k
      return
      end
c-------------------------------------------------------------
      subroutine inp2d1(w,nxp1,nyp1,xa,hx,ya,hy,x,y,ww,wx,wy)
c-------------------------------------------------------------
c
c  2d bilinear interpolation within a regular array
c
      real w(*)
c
      xx=(x-xa)/hx
      i1=xx
      i1=max0(0,min0(nxp1-2,i1))
      dx=xx-float(i1)
      yy=(y-ya)/hy
      j1=yy
      j1=max0(0,min0(nyp1-2,j1))
      dy=yy-float(j1)
      k1=1+i1+j1*nxp1
      k2=k1+1
      ww=1.e30
      if(w(k1).ge.1.e28.or.w(k2).ge.1.e28) return
      k3=k2+nxp1
      k4=k1+nxp1
      if(w(k3).ge.1.e28.or.w(k4).ge.1.e28) return
      z1=w(k1)
      z2=w(k2)-z1
      z3=w(k3)-z1
      z4=w(k4)-z1
      d=z3-z2-z4
c
      wx=z2+d*dy
      wy=z4+d*dx
      ww=z1+z2*dx+z4*dy+d*dx*dy
c
      return
      end
c------------------------------------------------------------------
      subroutine xparse(cc,nth,ifo,req,carg,iarg,rarg)
c------------------------------------------------------------------
c
c..... parse command line args
c
      character*(*) cc,req,carg
      character*256 rec
c
      iver=0
      lcc=len(cc)
      na=iargc()
      ith=0
      do 100 ia=1,na
      call getarg(ia,rec)
      if(rec(1:15).eq.'-xparse_verbose') iver=1
      lrec=jlen(rec,256)
      if(lrec.eq.lcc) then
      if(rec(1:lcc).eq.cc(1:lcc)) then
      if(iver.eq.1)write(0,'(a,a)')'xparse ,parsing:',rec(1:lcc)
c
      if(nth.le.0) then
      iarg=1
      return
      endif
c
      ith=ith+1
      if(ith.eq.nth) then
c
      if(ia.eq.na) then
      write(0,*)'parse error :  missing value for ',cc
      endif
      call getarg(ia+1,rec)
      lrec=jlen(rec,256)
      if(iver.eq.1)write(0,'(a,a)')'xparse, string:',rec(1:lrec)
      if(ifo.eq.1) then
      if(iver.eq.1)write(0,'(a,a)')'xparse, character:',rec(1:lrec)
      carg=rec
      else if(ifo.eq.2) then
      s=gets(rec)
      if(s.ge.0.0)iarg=s+0.1
      if(s.lt.0.0)iarg=s-0.1
      if(iver.eq.1)write(0,*)'xparse, integer:',iarg
      else if(ifo.eq.3) then
      rarg=gets(rec)
      if(iver.eq.1)write(0,*)'xparse, real:',rarg
      endif
      return
c
      endif
      endif
      endif
100   continue
c
      if(req(1:8).eq.'required') then
      write(0,*)'parse error :  cant find required arg:  ',cc
      call exit(0)
      endif
      end
c------------------------------------------
      real function gets(cc)
c------------------------------------------
c
c  decodes integer or floating f format
c  from character string 
c
      character*(*) cc
      nn=len(cc)
c
      fak=1.
      ief=0
      l1=0
      l2=0
      do 200 i=1,nn
      if(cc(i:i).eq.'e'.or.cc(i:i).eq.'E')ief=i
      if(cc(i:i).eq.'d'.or.cc(i:i).eq.'D')ief=i
      if(l1.eq.0.and.cc(i:i).ne.' ')l1=i
      if(cc(i:i).ne.' ')l2=i
      if(cc(i:i).eq.' '.and.l2.gt.0) goto 201
200   continue
201   continue
      nn=l2
      if(ief.gt.0) then
      lex=l2-ief
      iex=-9999999
      if(lex.eq.1)read(cc(ief+1:l2),'(i1)',err=900) iex
      if(lex.eq.2)read(cc(ief+1:l2),'(i2)',err=900) iex
      if(lex.eq.3)read(cc(ief+1:l2),'(i3)',err=900) iex
      if(lex.eq.4)read(cc(ief+1:l2),'(i4)',err=900) iex
      if(lex.eq.5)read(cc(ief+1:l2),'(i5)',err=900) iex
      if(iex.gt.-999999) then
      if(iex.lt.0)fak=1./( 10.**(-iex) )
      if(iex.gt.0)fak=10.**iex
      else
      write(0,*)'gets: cannot read ',cc
      endif
      nn=ief-1
      endif
c
      sig=1.
      ss=0.
      tt=1.
      ip=0
      do 100 l=1,nn
      if(cc(l:l).ne.' ') then
      if(cc(l:l).eq.'.') then
      ip=1
      else if(cc(l:l).eq.'-') then
      sig=-1.
      else
c      read(cc(l:l),'(i1)',err=900) ii
      ii=ichar(cc(l:l))-48
      if(ii.lt.0.or.ii.gt.9) goto 109
      if(ip.eq.0) then
      ss=10.*ss+float(ii)
      else
      tt=0.1*tt
      ss=ss+tt*float(ii)
      endif
      endif
      endif
100   continue
109   continue
      gets=ss*sig*fak
      return
900   continue
      write(0,*)' gets: error reading formatted integer:'
      write(0,*)nn
      write(0,'(a,a,a)')'$',cc,'$'
      call exit(0)
      end
'eeeeee'

pgf77 tmp.f

a.out -fin /mnt/md0/ftp/pub/tfr/b017.iof.072000

/bin/rm tmp.f

/bin/rm a.out