c
c
c
c----------------------------------------------
      program iocl4
c----------------------------------------------
c
c...... computes open/closed flux areas in the ionosphere
c       produces io_opcl variable ionosphere field
c       io_opcl=-1 for closed field, io_opcl=1 for
c       open (polar cap) field.
c       JR, 5/1995
c       add epar,vfl integrals, JR, 2006-07-26
c
c
      common /aaa1bx/bx(92000000)
      common /aaa1by/by(92000000)
      common /aaa1bz/bz(92000000)
      common /aaa1ex/ex(92000000)
      common /aaa1ey/ey(92000000)
      common /aaa1ez/ez(92000000)
      common /aaa2/gx(1000),gy(1000),gz(1000)
      common /aaa2ex/gxex(1000),gyex(1000),gzex(1000)
      common /aaa2ey/gxey(1000),gyey(1000),gzey(1000)
      common /aaa2ez/gxez(1000),gyez(1000),gzez(1000)
      common /aaa2bx/gxbx(1000),gybx(1000),gzbx(1000)
      common /aaa2by/gxby(1000),gyby(1000),gzby(1000)
      common /aaa2bz/gxbz(1000),gybz(1000),gzbz(1000)
      common /aaa3/mx,my,mz,nx,ny,nz,np,nt,istep
      common /aaa5/dpi,dti,the1,the2
      common /aaa6opcl/opcl(121,361) 
      common /aaa6absepar/absepar(121,361) 
      common /aaa6eparop/eparop(121,361) 
      common /aaa6eparcl/eparcl(121,361) 
      common /aaa6epar1/epar1(121,361) 
      common /aaa6vfl/vfl(121,361) 
      common /ca4/rad,pi,pih,pi2,deg
      common /aaa9/ f1,f2,rec,rec2,cdate,l1,l2
      character*80 f1,f2,rec,rec2,cdate,l1,l2
      common /bbb9/idate(5),istat(3)
      integer ida(6)
      deg=45.0/atan(1.0)
c
      nsamp=1
      call xparse('-sample',1,2,'optional',rec,nsamp,rdum)
c
c.... dipole time
      call xparse('-dip',1,1,'required',rec,idum,rdum)
      call crcol(rec,80,bx,nb)
      iy=bx(1)+0.1
      mo=bx(2)+0.1
      id=bx(3)+0.1
      ih=bx(4)+0.1
      mi=bx(5)+0.1
      se=bx(6)
      call timstr(rec,iy,mo,id,ih,mi,se)
      call cotr_set(iy,mo,id,ih,mi,se)
      write(0,*)'iocl4 dipole: ',rec(1:jlen(rec,80))
c
      call xparse('-time',1,1,'required',rec,idum,rdum)
      call ccfield(':',rec,1,rec2)
      ida(1)=int(gets(rec2)+0.1)
      call ccfield(':',rec,2,rec2)
      ida(2)=int(gets(rec2)+0.1)
      call ccfield(':',rec,3,rec2)
      ida(3)=int(gets(rec2)+0.1)
      call ccfield(':',rec,4,rec2)
      ida(4)=int(gets(rec2)+0.1)
      call ccfield(':',rec,5,rec2)
      ida(5)=int(gets(rec2)+0.1)
      call ccfield(':',rec,6,rec2)
      ida(6)=int(gets(rec2)+0.1)
      write(0,'(a,6i5)')'iocl4:date: ',ida
c
c.... get grid
      call xparse('-grid',1,1,'required',f1,idum,rdum)
      open(13,file=f1,status='old')
      l1='any'
      call getf11(13,gx,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid x ',nnx,gx(1),gx(nnx)
      if(nnx.le.0) stop
      nx=nnx
      l1='any'
      call getf11(13,gy,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid y ',nny,gy(1),gy(nny)
      if(nny.le.0) stop
      ny=nny
      l1='any'
      call getf11(13,gz,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid z ',nnz,gz(1),gz(nnz)
      if(nnz.le.0) stop
      nz=nnz
      l1='gx-bx'
      call getf11(13,gxbx,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid gx-bx ',nnx,gxbx(1),gxbx(nnx)
      l1='gy-bx'
      call getf11(13,gybx,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid gy-bx ',nny,gybx(1),gybx(nny)
      l1='gz-bx'
      call getf11(13,gzbx,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid gz-bx ',nnz,gzbx(1),gzbx(nnz)
      l1='gx-by'
      call getf11(13,gxby,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid gx-by ',nnx,gxby(1),gxby(nnx)
      l1='gy-by'
      call getf11(13,gyby,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid gy-by ',nny,gyby(1),gyby(nny)
      l1='gz-by'
      call getf11(13,gzby,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid gz-by ',nnz,gzby(1),gzby(nnz)
      l1='gx-bz'
      call getf11(13,gxbz,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid gx-bz ',nnx,gxbz(1),gxbz(nnx)
      l1='gy-bz'
      call getf11(13,gybz,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid gy-bz ',nny,gybz(1),gybz(nny)
      l1='gz-bz'
      call getf11(13,gzbz,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid gz-bz ',nnz,gzbz(1),gzbz(nnz)
      l1='gx-ex'
      call getf11(13,gxex,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid gx-ex ',nnx,gxex(1),gxex(nnx)
      l1='gy-ex'
      call getf11(13,gyex,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid gy-ex ',nny,gyex(1),gyex(nny)
      l1='gz-ex'
      call getf11(13,gzex,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid gz-ex ',nnz,gzex(1),gzex(nnz)
      l1='gx-ey'
      call getf11(13,gxey,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid gx-ey ',nnx,gxey(1),gxey(nnx)
      l1='gy-ey'
      call getf11(13,gyey,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid gy-ey ',nny,gyey(1),gyey(nny)
      l1='gz-ey'
      call getf11(13,gzey,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid gz-ey ',nnz,gzey(1),gzey(nnz)
      l1='gx-ez'
      call getf11(13,gxez,nnx,l1,l2,idum)
      write(0,*)'iocl4 got grid gx-ez ',nnx,gxez(1),gxez(nnx)
      l1='gy-ez'
      call getf11(13,gyez,nny,l1,l2,idum)
      write(0,*)'iocl4 got grid gy-ez ',nny,gyez(1),gyez(nny)
      l1='gz-ez'
      call getf11(13,gzez,nnz,l1,l2,idum)
      write(0,*)'iocl4 got grid gz-ez ',nnz,gzez(1),gzez(nnz)
      close(13)
c
c.... 3d file with field
      call xparse('-f',1,1,'required',f1,idum,rdum)
      open(13,file=f1,status='old')
c
c....  get fields , face centered B and edge centered E
      l1='bx1'
      call getf31(13,bx,mx,my,mz,l1,l2,istep)
      write(0,*)'iocl4 got bx ',mx,my,mz,istep 
      l1='by1'
      call getf31(13,by,mx,my,mz,l1,l2,istep)
      write(0,*)'iocl4 got by ',mx,my,mz,istep 
      l1='bz1'
      call getf31(13,bz,mx,my,mz,l1,l2,istep)
      write(0,*)'iocl4 got bz ',mx,my,mz,istep 
      l1='eflx'
      call getf31(13,ex,mx,my,mz,l1,l2,istep)
      write(0,*)'iocl4 got ex ',mx,my,mz,istep 
      l1='efly'
      call getf31(13,ey,mx,my,mz,l1,l2,istep)
      write(0,*)'iocl4 got ey ',mx,my,mz,istep 
      l1='eflz'
      call getf31(13,ez,mx,my,mz,l1,l2,istep)
      write(0,*)'iocl4 got ez ',mx,my,mz,istep 
      close(13)
c
c....... trace from ionosphere
c
      fnorth=0.0
      fsouth=0.0
      iop=0
      dphi=2.0*(3.1415927)/float(121-1)
      dtheta=(3.1415927)/float(361-1)
      do 550 it=1,361
      do 550 ip=1,121
      opcl(ip,it)=0.0 
      absepar(ip,it)=0.0 
      eparop(ip,it)=0.0 
      eparcl(ip,it)=0.0 
      epar1(ip,it)=0.0 
      vfl(ip,it)=0.0 
      opcl(ip,it)=-1.0
550   continue
      rr=3.7
      dp=2.0*(3.1415927)/float(nsamp*(121-1))
      dt=(3.1415927)/float(nsamp*(361-1))
      do 500 it=1,361
      write(0,'(a,$)')'.'
      do 500 ip=1,121
c
c..... subsampling, the max of E_par counts
      xepar=0.0
      xvfl=0.0
      do 505 ipp=0,nsamp-1
      do 505 itt=0,nsamp-1
      t=dt*float(nsamp*(it-1)+itt)
      p=dp*float(nsamp*(ip-1)+ipp)-(3.1415927)
c
      r=1.0
      r0=1.0/(sin(t)**2)
      arg=sqrt(rr/r0)
      if(arg.ge.0.95) goto 501
      t1=asin(arg)
      if(t.gt.((3.1415927)/2.0))t1=(3.1415927)-t1
      call radxyz(rr,p,t1,xsm,ysm,zsm)
      call cotr('sm ','mhd',xsm,ysm,zsm,xst,yst,zst)
      if(it.le.(361/2))dir=-0.01
      if(it.ge.(361/2))dir= 0.01
      call trace4(xst,yst,zst,xe,ye,ze,dir,iend,xxepar,xxvfl)
      if(iend.ne.0) then
      xvfl=xvfl+xxvfl
      if(xxepar.lt.0.0)xepar=min( -abs(xepar), xxepar)
      if(xxepar.gt.0.0)xepar=max( abs(xepar), xxepar)
c      write(0,'(4i6,8(1x,f9.2))')iend,nsamp,it,ip,deg*t,deg*p,xst,yst,zst,xe,ye,ze
      endif
505   continue
c
      epar1(ip,it)=xepar
      vfl(ip,it)=xvfl/float(nsamp*nsamp)
      absepar(ip,it)=abs(xepar)
      if(iend.eq.1) eparcl(ip,it)=xepar
      if(iend.eq.2) eparop(ip,it)=xepar
      if(iend.ne.1) then
      opcl(ip,it)=1.0
      iop=iop+1
      ff=(6371040.0)*dphi*(6371040.0)*dtheta*sin(t+0.5*dtheta)
      if(it.lt.(361/2))fnorth=fnorth+ff
      if(it.ge.(361/2))fsouth=fsouth+ff
      endif
501   continue
500   continue
c
      write(0,*)iop
      l1='io_opcl'
      call datf21(6,opcl,121,361,l1,l2,istep) 
      l1='io_absepar'
      call datf21(6,absepar,121,361,l1,l2,istep) 
      l1='io_eparop'
      call datf21(6,eparop,121,361,l1,l2,istep) 
      l1='io_eparcl'
      call datf21(6,eparcl,121,361,l1,l2,istep) 
      l1='io_epar1'
      call datf21(6,epar1,121,361,l1,l2,istep) 
      l1='io_vfl'
      call datf21(6,vfl,121,361,l1,l2,istep) 
      call system('/bin/rm -f tmp.iocl4.fluxarea')
      open(8,file='tmp.iocl4.fluxarea')
      write(8,'(6i5,1x,g12.5,1x,a)')ida,fnorth,' pc_area_n'
      write(8,'(6i5,1x,g12.5,1x,a)')ida,fsouth,' pc_area_s'
      write(0,'(6i5,1x,g12.5,1x,a)')ida,fnorth,
     *' open_flux_area_north'
      write(0,'(6i5,1x,g12.5,1x,a)')ida,fsouth,
     *' open_flux_area_south'
c      write(8,*)' open flux, northern hemisphere (m**2)',fnorth
c      write(8,*)' open flux, southern hemisphere (m**2)',fsouth
c      write(0,*)' open flux, northern hemisphere (m**2)',fnorth
c      write(0,*)' open flux, southern hemisphere (m**2)',fsouth
      close(8)
c
      stop
      end
c------------------------------------------------------------------
      subroutine trace4(x1,y1,z1,x2,y2,z2,dir,iend,xepar,xvfl)
c------------------------------------------------------------------
      common /aaa1bx/bx(92000000)
      common /aaa1by/by(92000000)
      common /aaa1bz/bz(92000000)
      common /aaa1ex/ex(92000000)
      common /aaa1ey/ey(92000000)
      common /aaa1ez/ez(92000000)
      common /aaa2/gx(1000),gy(1000),gz(1000)
      common /aaa2ex/gxex(1000),gyex(1000),gzex(1000)
      common /aaa2ey/gxey(1000),gyey(1000),gzey(1000)
      common /aaa2ez/gxez(1000),gyez(1000),gzez(1000)
      common /aaa2bx/gxbx(1000),gybx(1000),gzbx(1000)
      common /aaa2by/gxby(1000),gyby(1000),gzby(1000)
      common /aaa2bz/gxbz(1000),gybz(1000),gzbz(1000)
      common /aaa3/mx,my,mz,nx,ny,nz,np,nt,istep
      common /aaa5/dpi,dti,the1,the2
      common /aaa6opcl/opcl(121,361) 
      common /aaa6absepar/absepar(121,361) 
      common /aaa6eparop/eparop(121,361) 
      common /aaa6eparcl/eparcl(121,361) 
      common /aaa6epar1/epar1(121,361) 
      common /aaa6vfl/vfl(121,361) 
      common /ca4/rad,pi,pih,pi2,deg
      common /aaa9/ f1,f2,rec,rec2,cdate,l1,l2
      character*80 f1,f2,rec,rec2,cdate,l1,l2
      common /bbb9/idate(5),istat(3)
c..... trace till we hit
      iend=0
      xepar=0.0
      xvfl=0.0
      adir=abs(dir)*(6371040.0)
      deg=45.0/atan(1.0)
      x=x1
      y=y1
      z=z1
      itt=0
      do 100 itr=1,50000
      itt=itt+1
      r1=sqrt(x*x+y*y+z*z)
      call ipol3a(bx,nx,ny,nz,bbx,gxbx,gybx,gzbx,x,y,z,iout1)
      call ipol3a(ex,nx,ny,nz,eex,gxex,gyex,gzex,x,y,z,iout2)
      call ipol3a(by,nx,ny,nz,bby,gxby,gyby,gzby,x,y,z,iout1)
      call ipol3a(ey,nx,ny,nz,eey,gxey,gyey,gzey,x,y,z,iout2)
      call ipol3a(bz,nx,ny,nz,bbz,gxbz,gybz,gzbz,x,y,z,iout1)
      call ipol3a(ez,nx,ny,nz,eez,gxez,gyez,gzez,x,y,z,iout2)
      if(iout1+iout2.ne.0) goto 999
      r1=sqrt(x*x+y*y+z*z)
      if(r1.le.3.5) then
      x2=x
      y2=y
      z2=z
      iend=1
      return
      endif
      bb=sqrt(bbx*bbx+bby*bby+bbz*bbz)
      if(bb.le.0.0) goto 999
      s=dir/bb
      x=x+bbx*s
      y=y+bby*s
      z=z+bbz*s
      if(x.le.gx(3).or.x.ge.gx(nx-2)) iend=2
      if(y.le.gy(3).or.y.ge.gy(ny-2)) iend=2
      if(z.le.gz(3).or.z.ge.gz(nz-2)) iend=2
      xvfl=xvfl+adir/bb
      xepar=xepar+1.0e-3*adir*(eex*bbx+eey*bby+eez*bbz)/bb 
      if(iend.ne.0) then
      x2=x
      y2=y
      z2=z
      return
      endif
100   continue
999   continue
c      write(0,*)' no hit ',itt,bb,x,y,z,s
      return
      end
c
c--------------------------------------------
      subroutine tbdec(iu,a,n)
c--------------------------------------------
      real a(*)
      character*72 rec
      read(iu,*,end=999,err=999)n
      aa1=alog((1.0e-8))
      aa2=alog((1.0e10))
      daa=45.0*90.0*90.0/(aa2-aa1)
      daa=1.0/daa
      naa=45*90*90
      i=0
      do 200 k=1,9999999
      read(iu,'(a)',end=999,err=999) rec
      ll=0
      do 210 l=1,24
      i=i+1
      if(i.gt.n) return
      ll=ll+1
      ib1=ichar(rec(ll:ll))-34
      ll=ll+1
      ib2=ichar(rec(ll:ll))-34
      ll=ll+1
      ib3=ichar(rec(ll:ll))-34
      iaa=ib1+90*(ib2+90*ib3)
      ineg=0
      if(iaa.ge.naa)then
      ineg=1
      iaa=iaa-naa
      endif
      aa=iaa
      aa=(aa*daa)+aa1
      aa=exp(aa)
      if(ineg.gt.0)aa=-aa
      a(i)=aa
210   continue
200   continue
      return
999   continue
      n=-1
      return
      end
c--------------------------------------------
      subroutine tbenc(iu,a,n)
c--------------------------------------------
      real a(n)
      write(iu,*)n
      aa1=alog((1.0e-8))
      aa2=alog((1.0e10))
      daa=45.0*90.0*90.0/(aa2-aa1)
      naa=45*90*90
      do 200 i=1,n
      aa=abs(a(i))
      aa=amax1(aa,1.0e-12)
      aa=daa*(alog(aa)-aa1)
      iaa=int(aa+0.5)
      iaa=max0(0,min0(naa-1,iaa))
      if(a(i).lt.0.0)iaa=iaa+naa
      ib1=mod(iaa,90)
      iaa=iaa/90
      ib2=mod(iaa,90)
      ib3=iaa/90
      write(iu,'(a,a,a,$)')char(34+ib1),char(34+ib2),char(34+ib3)
      iw10=0
      if(mod(i,24).eq.0)then
      write(iu,'(a,$)')char(10)
      iw10=1
      endif
200   continue
      if(iw10.eq.0)write(iu,'(a,$)')char(10)
      return
      end
c-------------------------------------------------------------------
      integer function iifield(c,r,n)
c-------------------------------------------------------------------
      character*80 r
      character*1 c
      character*80 rec
      call ccfield(c,r,n,rec)
      xx=gets(rec)
      ii=int(abs(xx+0.5))
      if(xx.lt.0)ii=-ii
      iifield=ii
      return
      end
c-------------------------------------------------------------------
      real function rrfield(c,r,n)
c-------------------------------------------------------------------
      character*80 r
      character*1 c
      character*80 rec
      call ccfield(c,r,n,rec)
      xx=gets(rec)
      rrfield=xx
      return
      end
c-------------------------------------------------------------------
      subroutine ccfield(c,r,n,d)
c-------------------------------------------------------------------
      character*(*) r,d
      character*1 c
      d=' '
      ic=1
      k=0
      do 100 i=1,len(r)
      if(ic.eq.n.and.r(i:i).ne.c) then
      k=k+1
      d(k:k)=r(i:i)
      if(d(k:k).eq.'~')d(k:k)=' '
      endif
      if(i.gt.1) then
      if(r(i:i).eq.c.and.r(i-1:i-1).ne.c)ic=ic+1
      endif
c      write(0,*)'cf: ',i,k,'|',r(i:i),'|',ic
      if(ic.gt.n) goto 190
100   continue
190   continue
      return
      end
c-------------------------------------------------------------------
      subroutine cfield(c,r,n,d)
c-------------------------------------------------------------------
c.... get n'th occurence of 'c' seperated string in r, return in d
      character*80 r,d
      character*1 c
      d=' '
      ic=1
      k=0
      do 100 i=1,80
      if(ic.eq.n.and.r(i:i).ne.c) then
      k=k+1
      d(k:k)=r(i:i)
      if(d(k:k).eq.'~')d(k:k)=' '
      endif
      if(r(i:i).eq.c)ic=ic+1
100   continue
      return
      end
c----------------------------------------------------------------
      subroutine epoch1966(dsecs,iy,mo,id,ih,mi,sec,iimod)
c----------------------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-epoch1966
cdoc-------------------------------------------------------------------
cdoc-
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine epoch1966(dsec,iy,mo,id,ih,mi,se,iimod)
cdoc-  
cdoc-FUNCTION:
cdoc-  converts seconds since JAN 1, 1966 0UT to year,month,day,hour,
cdoc-  minute and second or vice versa
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  dsec (in/out): seconds of epoch 1966
cdoc-  iy,mo,id,ih,mi,se (in):  time in UT
cdoc-  iimod (in,integer): flag
cdoc-                     iimod=0: convert iy,... to dsec
cdoc-                     iimod!=0: convert dsec to iy,...
cdoc-
cdoc-NOTE:
cdoc-  works only for times after JAN 1, 1966, i.e. dsec positive
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      real*8 dsecs,ss,ds,dso
      integer rdays(0:11),jdays(0:11)
      data rdays /31,28,31,30,31,30,31,31,30,31,30,31/
      data jdays /31,29,31,30,31,30,31,31,30,31,30,31/
c
      if(iimod.eq.0) then
      jj=iy
      if(jj.lt.100)jj=jj+1900
      ss=0.0
      do 100 i=1966,jj-1
      if(mod(i,4).eq.0) ss=ss+366.0d0*24.0d0*3600.0d0
      if(mod(i,4).ne.0) ss=ss+365.0d0*24.0d0*3600.0d0
100   continue
      do 110 i=1,mo-1
      if( mod(jj,4).eq.0 ) ss=ss+dble(jdays(i-1))*24.0d0*3600.0d0
      if( mod(jj,4).ne.0 ) ss=ss+dble(rdays(i-1))*24.0d0*3600.0d0
110   continue
      ss=ss+dble(id-1)*3600.0d0*24.0d0
      ss=ss+dble(ih)*3600.0d0
      ss=ss+dble(mi)*60.0d0
      ss=ss+dble(sec)
      dsecs=ss
c
      else
      ss=dsecs
      ds=0.0d0
      dso=0.0d0
      jj=1966
      do 200 i=0,200
      if(mod(jj,4).eq.0) then
      ds=ds+(366.0d0*24.0d0*3600.0d0)
      else
      ds=ds+(365.0d0*24.0d0*3600.0d0)
      endif
      if(ds.gt.ss) goto 210
      dso=ds
      jj=jj+1
200   continue
210   continue
      ss=ss-dso
      iy=jj-1900
      ij=0
      if(mod(jj,4).eq.0)ij=1
      ds=0.0d0
      dso=0.0d0
      jj=0
      do 220 i=0,11
      if(ij.eq.0) then
      ds=ds+dble(rdays(i))*24.0d0*3600.0d0
      else
      ds=ds+dble(jdays(i))*24.0d0*3600.0d0
      endif
      if(ds.gt.ss) goto 230
      dso=ds
      jj=jj+1
220   continue
230   continue
      ss=ss-dso
      mo=jj+1
      ds=ss/(24.0d0*3600.0d0)
      id=int(ds)
      ss=ss-dble(id)*24.0d0*3600.0d0
      id=id+1
      ds=ss/3600.0d0
      ih=int(ds)
      ss=ss-dble(ih)*3600.0d0
      ds=ss/60.0d0
      mi=int(ds)
      ss=ss-dble(mi)*60.0d0
      sec=ss
      endif
c
      return
      end
c---------------------------------------------------------------
      double precision function djul(iy,mo,id,ih,mi,se)
c---------------------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-djul
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  double precision function djul(iy,mo,id,ih,mi,se)
cdoc-  
cdoc-FUNCTION:
cdoc-  converts time to fractional julian days
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  iy,mo,id,ih,mi,se (in):  time in UT
cdoc-  djul (out,real*8):   fractional julian day
cdoc-
cdoc-NOTE:
cdoc-  expected to work for any time after 4500 BC
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      implicitdoubleprecision(a-h,o-z)
      real*4 se
      jy=iy
      if(jy.lt.100)jy=jy+1900
      ier=0
      if(jy.lt.1901) ier=1
      if(jy.gt.2049) ier=7
      if(mo.lt.1.or.mo.gt.12) ier=2
      if(id.lt.1.or.id.gt.31) ier=3
      if(ih.lt.0.or.ih.gt.23) ier=4
      if(mi.lt.0.or.mi.gt.63) ier=5
      if(se.lt.0.0.or.se.gt.60.0) ier=6
      if(ier.ne.0) then
      write(0,*)'error djul, ier= ',ier
      stop
      endif
      uthour=dble(float(ih))+(dble(float(mi))/60.0d0)+dble(se)/
     *3600.0d0
      i1=367*jy
      i2=(mo+9)/12
      i2=(i2+jy)*7
      i2=i2/4
      i3=(275*mo)/9
      i4=100*jy+mo
      d1=1.0d0
      zero=0.0d0
      if((dble(i4)-190002.5d0).lt.zero)d1=(-d1)
      xjd=dble(i1) -dble(i2) +dble(i3) +dble(id) +1721013.5d0 +uthour
     */24.0d0 -0.5d0*d1 +0.5
      djul=xjd
      return
      end
c---------------------------------------------------------------
      subroutine mhdmlt(iy,mo,id,ih,mi,se,pmhd,tmhd,xmlt,colat)
c---------------------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-mhdmlt
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine mhdmlt(iy,mo,id,ih,mi,se,pmhd,tmhd,xmlt,colat)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert MHD logitude (-180,180) and colatitude (0,180)
cdoc-  to magnetic local time (0-24) and colatitude
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  iy,mo,id,ih,mi,se (in):  time in UT
cdoc-  pmhd (in,real):  deg longitude in MHD system (-180.0,180.0)
cdoc-  tmhd (in,real):  deg colatitude in MHD system (0.0,180.0)
cdoc-  xmlt (out,real): hours magnetic local time (0.0-24.0)
cdoc-  colat (out,real): deg magnetic colatitude (0.0,180.0)
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      call cotr_set(iy,mo,id,ih,mi,se)
      call degxyz(1.0,pmhd,tmhd,x1,y1,z1)
      call cotr('mhd','sm ',x1,y1,z1,x2,y2,z2)
      call xyzdeg(x2,y2,z2,rr,pp,tt)
      colat=tt
      xmlt=(pp+180.0)/15.0
      return
      end
c---------------------------------------------------------
      subroutine diptil(iy,mo,id,ih,mi,se,degsun,degy)
c---------------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-diptil
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine diptil(iy,mo,id,ih,mi,se,degsun,degy)
cdoc-  
cdoc-FUNCTION:
cdoc-  compute dipole tilt angles for a given time
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  iy,mo,id,ih,mi,se (in):  time in UT
cdoc-  degsun (out,real): deg tilt in GSE x-z plane
cdoc-  degy (out,real):   deg tilt in GSE y-z plane
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      data deg /57.29577951d0/
      data xsm,ysm,zsm /0.0,0.0,1.0/
      call cotr_set(iy,mo,id,ih,mi,se)
      call cotr('sm ','gse',xsm,ysm,zsm,xgse,ygse,zgse)
      degsun=deg*atan2(xgse,zgse)
      degy=deg*atan2(ygse,zgse)
      return
      end
c----------------------------------------------------
      subroutine xyzrad(x,y,z,r,p,t)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-xyzrad
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine xyzrad(x,y,z,r,p,t)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert from cartesian to polar coordinates, angles in radian
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  x,y,z (in,real): cartesian coordinates
cdoc-  r,p,t (out,real): polar coordinates, -pi<=p<=pi, 0<=t<=pi
cdoc-                    i.e. colatitude system
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      r=sqrt(x*x+y*y+z*z)
      p=0.0
      if(abs(x)+abs(y).gt.0.0)p=atan2(y,x)
      t=0.0
      if(r.gt.0.0)t=acos(z/r)
      return
      end
c----------------------------------------------------
      subroutine xyzdeg(x,y,z,r,p,t)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-xyzdeg
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine xyzdeg(x,y,z,r,p,t)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert from cartesian to polar coordinates, angles in degrees
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  x,y,z (in,real): cartesian coordinates
cdoc-  r,p,t (out,real): polar coordinates, 
cdoc-                    -180.0<=p<=180.0, 0<=t<=180.0
cdoc-                    i.e. colatitude system
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      data deg /57.29577951d0/
      r=sqrt(x*x+y*y+z*z)
      p=0.0
      if(abs(x)+abs(y).gt.0.0)p=deg*atan2(y,x)
      t=0.0
      if(r.gt.0.0)t=deg*acos(z/r)
      return
      end
c----------------------------------------------------
      subroutine radxyz(r,p,t,x,y,z)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-radxyz
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine radxyz(r,p,t,x,y,z)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert from polar coordinates (angles in radian) to cartesian
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  r,p,t (in,real): polar coordinates, -pi<=p<=pi, 0<=t<=pi
cdoc-                    i.e. colatitude system
cdoc-  x,y,z (out,real): cartesian coordinates
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      st=sin(t)
      x=r*cos(p)*st
      y=r*sin(p)*st
      z=r*cos(t)
      return
      end
c----------------------------------------------------
      subroutine degxyz(r,p,t,x,y,z)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-degxyz
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine degxyz(r,p,t,x,y,z)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert from polar coordinates (angles in degrees) to cartesian
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  r,p,t (in,real): polar coordinates,
cdoc-                    -180.0<=p<=180.0, 0<=t<=180.0
cdoc-                    i.e. colatitude system
cdoc-  x,y,z (out,real): cartesian coordinates
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      data rad /17.45329252d-3/
      pp=p*rad
      tt=t*rad
      st=sin(tt)
      x=r*cos(pp)*st
      y=r*sin(pp)*st
      z=r*cos(tt)
      return
      end
c----------------------------------------------------
      subroutine daymon(iy,jd,mo,id)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-daymon
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine daymon(iy,jd,mo,id)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert date from julian day system to day - month system
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  iy (in,integer):  year (nn or nnnn notation)
cdoc-  jd (in,integer):  julian day (1<=jd<=366), january 1 <> jd=1
cdoc-  mo (out,integer): month
cdoc-  id (out,integer): day of month
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      integer mdays(12,2)
      data mdays / 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 
     *334, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 /
      leap=1
      if(mod(iy,4).eq.0)leap=2
      k=1
      do 100 i=1,12
      if(jd.le.mdays(i,leap)) goto 190
      k=i
100   continue
190   continue
      mo=k
      id=jd-mdays(k,leap)
      return
      end
c----------------------------------------------------
      integer function julday(iy,mo,id)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-julday
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  integer function julday(iy,mo,id)
cdoc-  
cdoc-FUNCTION:
cdoc-  convert date from day - month system to julian day system
cdoc-
cdoc-CALL SEQUENCE:
cdoc-
cdoc-PARAMETERS:
cdoc-  iy (in,integer):  year (nn or nnnn notation)
cdoc-  mo (in,integer): month
cdoc-  id (in,integer): day of month
cdoc-  julday (out,integer): julian day of year
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
      integer mdays(12,2)
      data mdays / 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 
     *334, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 /
      jmo=max0(1,min0(12,mo))
      leap=1
      if(mod(iy,4).eq.0)leap=2
      julday=mdays(jmo,leap)+id
      return
      end
c----------------------------------------------------
      subroutine cotr_set(iy,mo,id,ih,mi,se)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-cotr_set
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine cotr_set(iy,mo,id,ih,mi,se)
cdoc-  
cdoc-FUNCTION:
cdoc-  computes coordinate transformation matrices
cdoc-  for cotr() for a given time
cdoc-
cdoc-CALL SEQUENCE:
cdoc-  must be called before cotr()
cdoc-
cdoc-PARAMETERS:
cdoc-  iy (in,integer):  year (nn or nnnn notation) (1900<iy<2000)
cdoc-  mo (in,integer):  month             (1<=mo<=12)
cdoc-  id (in,integer):  day of month      (1<=id<=31)
cdoc-  ih (in,integer):  hour of day       (0<=ih<24)
cdoc-  mi (in,integer):  minute of hour    (0<=mi<60)
cdoc-  se (in,real):     seconds of minute (0.0<=se<60.0)
cdoc-
cdoc-NOTES:
cdoc-  1.  no output, transformation matrices are stored internally
cdoc-      for use by cotr()
cdoc-  2.  dipole coordinates are computed from IGRF g10,g11,h11
cdoc-      coefficients between 1945 and 1990 with linear time
cdoc-      interpolation.  Before 1945 coefficients for 1945 are
cdoc-      used and after 1990 coefficients for 1990 are used.
cdoc-  3.  internal arithmetic is handled in double precision
cdoc-
cdoc-  4.  for iy=1967 mo=1 id=1 ih=0 mi=0 se=any all transformations
cdoc-      are set to identity
cdoc-
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
cdoc-REFERENCES:
cdoc-
cdoc-        Almanac for Computers 1991
cdoc-        Nautical Almanac Office
cdoc-        US Naval Observatory
cdoc-        Washington, DC
cdoc-
cdoc-        M A Hapgood
cdoc-        Space physics coordinate transformations: A user guide
cdoc-        Planet. Space Sci., 40, 711, 1992
cdoc-
      implicitdoubleprecision(a-h,o-z)
      common /cotrd1/t0,iscalled
      common /cotrd2/tmpmat(3,3,2),cmat(3,3),pmat(3,3,-7:7)
      common /cotrd3/ccfr,ccto
      character*3 ccfr,ccto
      integer iy,mo,id,ih,mi
      real*4 se
      dimension tig(12),g10(12),g11(12),h11(12)
      integer mdays(12,2)
      data mdays / 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 
     *334, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 /
      data tig( 1),g10( 1),g11( 1),h11( 1)/ 1000.00,-30594.00, -
     *2285.00, 5810.00/
      data tig( 2),g10( 2),g11( 2),h11( 2)/ 1945.00,-30594.00, -
     *2285.00, 5810.00/
      data tig( 3),g10( 3),g11( 3),h11( 3)/ 1950.00,-30554.00, -
     *2250.00, 5815.00/
      data tig( 4),g10( 4),g11( 4),h11( 4)/ 1955.00,-30500.00, -
     *2215.00, 5820.00/
      data tig( 5),g10( 5),g11( 5),h11( 5)/ 1960.00,-30421.00, -
     *2169.00, 5791.00/
      data tig( 6),g10( 6),g11( 6),h11( 6)/ 1965.00,-30334.00, -
     *2119.00, 5776.00/
      data tig( 7),g10( 7),g11( 7),h11( 7)/ 1970.00,-30220.00, -
     *2068.00, 5737.00/
      data tig( 8),g10( 8),g11( 8),h11( 8)/ 1975.00,-30100.00, -
     *2013.00, 5675.00/
      data tig( 9),g10( 9),g11( 9),h11( 9)/ 1980.00,-29992.00, -
     *1956.00, 5604.00/
      data tig(10),g10(10),g11(10),h11(10)/ 1985.00,-29873.00, -
     *1905.00, 5500.00/
      data tig(11),g10(11),g11(11),h11(11)/ 1990.00,-29775.00, -
     *1851.00, 5411.00/
      data tig(12),g10(12),g11(12),h11(12)/ 3000.00,-29775.00, -
     *1851.00, 5411.00/
      data pih /1.570796372d0/
      data twopi /6.283185307179586476925287/
      data rad /17.45329252d-3/
      data deg /57.29577951d0/
      data two /2.0d0/
      data one /1.0d0/
      data zero /0.0d0/
c
c..... check input consistency
c
      jy=iy
      if(jy.lt.149)jy=jy+1900
      ier=0
      if(jy.lt.1901) ier=1
      if(jy.gt.2049) ier=7
      if(mo.lt.1.or.mo.gt.12) ier=2
      if(id.lt.1.or.id.gt.31) ier=3
      if(ih.lt.0.or.ih.gt.23) ier=4
      if(mi.lt.0.or.mi.gt.63) ier=5
      if(se.lt.0.0.or.se.gt.60.0) ier=6
      if(ier.ne.0) then
      write(0,*)'error cotr_set, ier= ',ier
      write(0,*)iy,jy,mo,id,ih,mi,se
      stop
      endif
c
      iscalled=1
c
c.....  universal time of day in hours (UT1)
c
      uthour=dble(float(ih))+(dble(float(mi))/60.0d0)+dble(se)/
     *3600.0d0
c
c..... full julian day
c
      i1=367*jy
      i2=(mo+9)/12
      i2=(i2+jy)*7
      i2=i2/4
      i3=(275*mo)/9
      i4=100*jy+mo
      d1=one
      if((dble(i4)-190002.5d0).lt.zero)d1=(-d1)
      xjd=dble(i1) -dble(i2) +dble(i3) +dble(id) +1721013.5d0 +uthour
     */24.0d0 -0.5d0*d1 +0.5
c
c..... julian day at 0 UT
c
      xjd0=xjd-(uthour/24.0d0)
c
c.....  modified julian day
c
      xmjd=xjd-2400000.5d0
      mjd=int(xmjd)
c
c.....  set to identity transformation if iy=1967,....
c       except for the MHD transformation which still flips x- y- axes
c
      if(jy.le.1967) then
      do 300 k=-7,7
      do 300 j=1,3
      do 300 i=1,3
      pmat(i,j,k)=0.0d0
      if(i.eq.j)pmat(i,j,k)=1.0d0
300   continue
      do 310 i=1,3
      do 310 j=1,3
      pmat(j,i,6)=zero
310   continue
      pmat(1,1,6)=(-1.0d0)
      pmat(2,2,6)=(-1.0d0)
      pmat(3,3,6)=(1.0d0)
      do 320 i=1,3
      do 320 j=1,3
      pmat(j,i,-6)=pmat(i,j,6)
320   continue
      return
c      
      endif
c
c
c..... solar coordinates
c
      t=(xjd-2451545.0d0)/36525.0d0
      t0=(xjd0-2451545.0d0)/36525.0d0
      xm=357.528d0+35999.050d0*t
      xl=280.460d0+36000.772d0*t
      xm=dmod(xm+360.0d3,360.0d0)
      xl=dmod(xl+360.0d3,360.0d0)
      xls=xl +((1.915d0-0.0048*t0)*dsin(rad*xm)) +(0.020d0*dsin(2.0d0
     **rad*xm))
      eps=23.439d0-0.013d0*t
      eps=dmod(eps+360.0d3,360.0d0)
c
c..... greenwich mean siderial time
c
      gmsth=6.69737456d0 +2400.051336d0*t0 +0.0000258622d0*t0*t0 +
     *1.002737909d0*uthour
      gmsth=dmod(gmsth+24.0d3,24.0d0)
      gmstd=15.0d0*gmsth
      gang=rad*gmstd
c
      eps=rad*eps
      xls=rad*xls
c
c.... gei to geo
c
      do 94971 i=1,3
      do 94971 j=1,3
      pmat(j,i,1)=zero
94971 continue
      tmpsi=dsin(gang)
      tmpco=dcos(gang)
      pmat(1,1,1)=tmpco
      pmat(2,2,1)=tmpco
      pmat(3,3,1)=one
      pmat(2,1,1)=-tmpsi
      pmat(1,2,1)=tmpsi
      do 94913 i=1,3
      do 94913 j=1,3
      pmat(j,i,-1)=pmat(i,j,1)
94913 continue
c
c.... gei to gse
c
      do 94884 i=1,3
      do 94884 j=1,3
      tmpmat(j,i,1)=zero
94884 continue
      tmpsi=dsin(xls)
      tmpco=dcos(xls)
      tmpmat(1,1,1)=tmpco
      tmpmat(2,2,1)=tmpco
      tmpmat(3,3,1)=one
      tmpmat(2,1,1)=-tmpsi
      tmpmat(1,2,1)=tmpsi
c
      do 94826 i=1,3
      do 94826 j=1,3
      tmpmat(j,i,2)=zero
94826 continue
      tmpsi=dsin(eps)
      tmpco=dcos(eps)
      tmpmat(1,1,2)=one
      tmpmat(2,2,2)=tmpco
      tmpmat(3,3,2)=tmpco
      tmpmat(3,2,2)=-tmpsi
      tmpmat(2,3,2)=tmpsi
c
      do 94768 i=1,3
      do 94768 j=1,3
      pmat(i,j,2)=tmpmat(i,1,1)*tmpmat(1,j,2) +tmpmat(i,2,1)*tmpmat(2
     *,j,2) +tmpmat(i,3,1)*tmpmat(3,j,2)
94768 continue
      do 94739 i=1,3
      do 94739 j=1,3
      pmat(j,i,-2)=pmat(i,j,2)
94739 continue
c
c
c .... gse to gsm
c
      xy=100.0d0*(20.0d0+t) 
      do 100 k=1,12-1
      if(xy.ge.tig(k).and.xy.le.tig(k+1)) then
      tmp=(xy-tig(k))/(tig(k+1)-tig(k))
      gg10=(one-tmp)*g10(k)+tmp*g10(k+1)
      gg11=(one-tmp)*g11(k)+tmp*g11(k+1)
      hh11=(one-tmp)*h11(k)+tmp*h11(k+1)
      goto 190
      endif
100   continue
190   continue
      xxl=datan(hh11/gg11)
      xxp=pih-dasin( (gg11*dcos(xxl)+hh11*dsin(xxl))/gg10 )
c
c      write(0,*)' dipole lambda  ',xxl*45.0d0/datan(1.0d0)
c      write(0,*)' dipole phi     ',xxp*45.0d0/datan(1.0d0)
c
c
c
      qxg=dcos(xxp)*dcos(xxl)
      qyg=dcos(xxp)*dsin(xxl)
      qzg=dsin(xxp)
      tmpx=pmat(1,1,-1)*qxg+pmat(1,2,-1)*qyg+pmat(1,3,-1)*qzg
      tmpy=pmat(2,1,-1)*qxg+pmat(2,2,-1)*qyg+pmat(2,3,-1)*qzg
      tmpz=pmat(3,1,-1)*qxg+pmat(3,2,-1)*qyg+pmat(3,3,-1)*qzg
      qx=pmat(1,1,2)*tmpx+pmat(1,2,2)*tmpy+pmat(1,3,2)*tmpz
      qy=pmat(2,1,2)*tmpx+pmat(2,2,2)*tmpy+pmat(2,3,2)*tmpz
      qz=pmat(3,1,2)*tmpx+pmat(3,2,2)*tmpy+pmat(3,3,2)*tmpz
      psi=datan2(qy,qz)
c      write(0,*)' dipole psi     ',psi*45.0d0/datan(1.0d0)
      do 94652 i=1,3
      do 94652 j=1,3
      pmat(j,i,-3)=zero
94652 continue
      tmpsi=dsin(psi)
      tmpco=dcos(psi)
      pmat(1,1,-3)=one
      pmat(2,2,-3)=tmpco
      pmat(3,3,-3)=tmpco
      pmat(3,2,-3)=-tmpsi
      pmat(2,3,-3)=tmpsi
      do 94594 i=1,3
      do 94594 j=1,3
      pmat(j,i,3)=pmat(i,j,-3)
94594 continue
c
c.... gsm to sm
c
      xmu=datan(qx/sqrt(qy*qy+qz*qz))
      do 94565 i=1,3
      do 94565 j=1,3
      pmat(j,i,-4)=zero
94565 continue
      tmpsi=dsin(xmu)
      tmpco=dcos(xmu)
      pmat(1,1,-4)=tmpco
      pmat(2,2,-4)=one
      pmat(3,3,-4)=tmpco
      pmat(3,1,-4)=-tmpsi
      pmat(1,3,-4)=tmpsi
      do 94507 i=1,3
      do 94507 j=1,3
      pmat(j,i,4)=pmat(i,j,-4)
94507 continue
c
c.... geo to mag
c
      w1=xxp-pih
      do 94478 i=1,3
      do 94478 j=1,3
      tmpmat(j,i,1)=zero
94478 continue
      do 94449 i=1,3
      do 94449 j=1,3
      tmpmat(j,i,2)=zero
94449 continue
      tmpsi=dsin(w1)
      tmpco=dcos(w1)
      tmpmat(1,1,1)=tmpco
      tmpmat(2,2,1)=one
      tmpmat(3,3,1)=tmpco
      tmpmat(3,1,1)=-tmpsi
      tmpmat(1,3,1)=tmpsi
      tmpsi=dsin(xxl)
      tmpco=dcos(xxl)
      tmpmat(1,1,2)=tmpco
      tmpmat(2,2,2)=tmpco
      tmpmat(3,3,2)=one
      tmpmat(2,1,2)=-tmpsi
      tmpmat(1,2,2)=tmpsi
      do 94362 i=1,3
      do 94362 j=1,3
      pmat(i,j,5)=tmpmat(i,1,1)*tmpmat(1,j,2) +tmpmat(i,2,1)*tmpmat(2
     *,j,2) +tmpmat(i,3,1)*tmpmat(3,j,2)
94362 continue
      do 94333 i=1,3
      do 94333 j=1,3
      pmat(j,i,-5)=pmat(i,j,5)
94333 continue
c
c.....  gse to mhd
c
      do 94304 i=1,3
      do 94304 j=1,3
      pmat(j,i,6)=zero
94304 continue
      pmat(1,1,6)=(-1.0d0)
      pmat(2,2,6)=(-1.0d0)
      pmat(3,3,6)=(1.0d0)
      do 94275 i=1,3
      do 94275 j=1,3
      pmat(j,i,-6)=pmat(i,j,6)
94275 continue
c      
      return
      end
c----------------------------------------------------
      subroutine cotr(cfr,cto,x1,y1,z1,x2,y2,z2)
c----------------------------------------------------
cdoc-
cdoc-------------------------------------------------------------------
cdoc-cotr
cdoc-------------------------------------------------------------------
cdoc-
cdoc-SYNOPSIS:
cdoc-  subroutine cotr(cfr,cto,x1,y1,z1,x2,y2,z2)
cdoc-  
cdoc-FUNCTION:
cdoc-  transform vector (x1,y1,z1) in coordinate system 'cfr'
cdoc-  to vector (x2,y2,z2) in coordinate system 'cto'
cdoc-  coordinate systems 'cfr' and 'cto' can be any of the following:
cdoc-  'gei'  :  geocentric equatorial inertial
cdoc-  'geo'  :  geographic
cdoc-  'gse'  :  geocentric solar ecliptic
cdoc-  'gsm'  :  geocentric solar magnetospheric
cdoc-  'sm '  :  solar magnetic
cdoc-  'mag'  :  geomagnetic
cdoc-  'mhd'  :  global mhd simulation system, like 'gse' with
cdoc-            x and y axes mirrored
cdoc-
cdoc-CALL SEQUENCE:
cdoc-  cotr_set() must be called before the first call to 
cdoc-  cotr() or if the transformation is for a different time
cdoc-
cdoc-PARAMETERS:
cdoc-  cfr (in,character*3)   'from' coordinate system, for ex. 'gse'
cdoc-  cto (in,character*3)   'to' coordinate system, for ex. 'gsm'
cdoc-  x1,y1,z1 (in,real)     input vector
cdoc-  x2,y2,z2 (out,real)    transformed vector
cdoc-
cdoc-NOTES:
cdoc-  see NOTES and REFERENCES for cotr_set()
cdoc-
cdoc-FILES:
cdoc-  none
cdoc-
cdoc-AUTHOR:
cdoc-  Joachim Raeder, IGPP/UCLA   June 1994
cdoc-
c
c
      implicitdoubleprecision(a-h,o-z)
      character*(*) cfr,cto
      real*4 x1,y1,z1,x2,y2,z2
      common /cotrd1/t0,iscalled
      common /cotrd2/tmpmat(3,3,2),cmat(3,3),pmat(3,3,-7:7)
      common /cotrd3/ccfr,ccto
      character*3 ccfr,ccto
      character*3 clfr,clto
      integer imat(7,7,7)
      data icalled /0/
      data (imat(l,1,1),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
      data (imat(l,2,1),l=1,7)/-1, 0, 0, 0, 0, 0, 0/
      data (imat(l,3,1),l=1,7)/-2, 0, 0, 0, 0, 0, 0/
      data (imat(l,4,1),l=1,7)/-3,-2, 0, 0, 0, 0, 0/
      data (imat(l,5,1),l=1,7)/-4,-3,-2, 0, 0, 0, 0/
      data (imat(l,6,1),l=1,7)/-5,-1, 0, 0, 0, 0, 0/
      data (imat(l,7,1),l=1,7)/ 6,-2, 0, 0, 0, 0, 0/
      data (imat(l,1,2),l=1,7)/ 1, 0, 0, 0, 0, 0, 0/
      data (imat(l,2,2),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
      data (imat(l,3,2),l=1,7)/-2, 1, 0, 0, 0, 0, 0/
      data (imat(l,4,2),l=1,7)/-3,-2, 1, 0, 0, 0, 0/
      data (imat(l,5,2),l=1,7)/-4,-3,-2, 1, 0, 0, 0/
      data (imat(l,6,2),l=1,7)/-5, 0, 0, 0, 0, 0, 0/
      data (imat(l,7,2),l=1,7)/ 6,-2, 1, 0, 0, 0, 0/
      data (imat(l,1,3),l=1,7)/ 2, 0, 0, 0, 0, 0, 0/
      data (imat(l,2,3),l=1,7)/-1, 2, 0, 0, 0, 0, 0/
      data (imat(l,3,3),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
      data (imat(l,4,3),l=1,7)/-3, 0, 0, 0, 0, 0, 0/
      data (imat(l,5,3),l=1,7)/-4,-3, 0, 0, 0, 0, 0/
      data (imat(l,6,3),l=1,7)/-5,-1, 2, 0, 0, 0, 0/
      data (imat(l,7,3),l=1,7)/ 6, 0, 0, 0, 0, 0, 0/
      data (imat(l,1,4),l=1,7)/ 2, 3, 0, 0, 0, 0, 0/
      data (imat(l,2,4),l=1,7)/-1, 2, 3, 0, 0, 0, 0/
      data (imat(l,3,4),l=1,7)/ 3, 0, 0, 0, 0, 0, 0/
      data (imat(l,4,4),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
      data (imat(l,5,4),l=1,7)/-4, 0, 0, 0, 0, 0, 0/
      data (imat(l,6,4),l=1,7)/-5,-1, 2, 3, 0, 0, 0/
      data (imat(l,7,4),l=1,7)/ 6, 3, 0, 0, 0, 0, 0/
      data (imat(l,1,5),l=1,7)/ 2, 3, 4, 0, 0, 0, 0/
      data (imat(l,2,5),l=1,7)/-1, 2, 3, 4, 0, 0, 0/
      data (imat(l,3,5),l=1,7)/ 3, 4, 0, 0, 0, 0, 0/
      data (imat(l,4,5),l=1,7)/ 4, 0, 0, 0, 0, 0, 0/
      data (imat(l,5,5),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
      data (imat(l,6,5),l=1,7)/-5,-1, 2, 3, 4, 0, 0/
      data (imat(l,7,5),l=1,7)/ 6, 3, 4, 0, 0, 0, 0/
      data (imat(l,1,6),l=1,7)/ 1, 5, 0, 0, 0, 0, 0/
      data (imat(l,2,6),l=1,7)/ 5, 0, 0, 0, 0, 0, 0/
      data (imat(l,3,6),l=1,7)/-2, 1, 5, 0, 0, 0, 0/
      data (imat(l,4,6),l=1,7)/-3,-2, 1, 5, 0, 0, 0/
      data (imat(l,5,6),l=1,7)/-4,-3,-2, 1, 5, 0, 0/
      data (imat(l,6,6),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
      data (imat(l,7,6),l=1,7)/ 6,-2, 1, 5, 0, 0, 0/
      data (imat(l,1,7),l=1,7)/ 2,-6, 0, 0, 0, 0, 0/
      data (imat(l,2,7),l=1,7)/-1, 2, 6, 0, 0, 0, 0/
      data (imat(l,3,7),l=1,7)/-6, 0, 0, 0, 0, 0, 0/
      data (imat(l,4,7),l=1,7)/-3,-6, 0, 0, 0, 0, 0/
      data (imat(l,5,7),l=1,7)/-4,-3,-6, 0, 0, 0, 0/
      data (imat(l,6,7),l=1,7)/-5,-1, 2,-6, 0, 0, 0/
      data (imat(l,7,7),l=1,7)/ 0, 0, 0, 0, 0, 0, 0/
c
      if(icalled.eq.0.or.iscalled.ne.0) then
      icalled=1
      ccfr='XXX'
      ccto='xyz'
      endif
      iscalled=0
c
      xx1=dble(x1)
      yy1=dble(y1)
      zz1=dble(z1)
      clfr=cfr
      clto=cto
c
      if(clfr.ne.ccfr.or.clto.ne.ccto) then
      ifr=0
      if(clfr.eq.'GEI'.or.clfr.eq.'gei')ifr=1
      if(clfr.eq.'GEO'.or.clfr.eq.'geo')ifr=2
      if(clfr.eq.'GSE'.or.clfr.eq.'gse')ifr=3
      if(clfr.eq.'GSM'.or.clfr.eq.'gsm')ifr=4
      if(clfr.eq.'SM '.or.clfr.eq.'sm ')ifr=5
      if(clfr.eq.'MAG'.or.clfr.eq.'mag')ifr=6
      if(clfr.eq.'MHD'.or.clfr.eq.'mhd')ifr=7
      ito=0
      if(clto.eq.'GEI'.or.clto.eq.'gei')ito=1
      if(clto.eq.'GEO'.or.clto.eq.'geo')ito=2
      if(clto.eq.'GSE'.or.clto.eq.'gse')ito=3
      if(clto.eq.'GSM'.or.clto.eq.'gsm')ito=4
      if(clto.eq.'SM '.or.clto.eq.'sm ')ito=5
      if(clto.eq.'MAG'.or.clto.eq.'mag')ito=6
      if(clto.eq.'MHD'.or.clto.eq.'mhd')ito=7
      if(ifr.lt.1.or.ito.lt.1) then
      write(0,*)' error cotr: ifr,ito ',ifr,ito
      stop
      endif
c
      ccfr=clfr
      ccto=clto
      do 100 i=1,3
      do 100 j=1,3
      cmat(i,j)=0.0d0
      if(i.eq.j)cmat(i,j)=1.0d0
100   continue
      do 200 k=1,7
      kmat=imat(k,ifr,ito)
      if(kmat.eq.0) goto 290
      do 210 i=1,3
      do 210 j=1,3
c      tmpmat(i,j,1)=cmat(i,1)*pmat(1,j,kmat)
c     *             +cmat(i,2)*pmat(2,j,kmat)
c     *             +cmat(i,3)*pmat(3,j,kmat)
      tmpmat(i,j,1)=pmat(i,1,kmat)*cmat(1,j) +pmat(i,2,kmat)*cmat(2,j
     *) +pmat(i,3,kmat)*cmat(3,j)
210   continue
      do 220 i=1,3
      do 220 j=1,3
      cmat(i,j)=tmpmat(i,j,1)
220   continue
200   continue
290   continue
      endif
c
c      write(0,*)'cotr1: ',xx1,yy1,zz1
c      write(0,*)'cotr3: ',ifr,ito
c      write(0,*)'cotr3: ',cmat
c
      xx2=cmat(1,1)*xx1+cmat(1,2)*yy1+cmat(1,3)*zz1
      yy2=cmat(2,1)*xx1+cmat(2,2)*yy1+cmat(2,3)*zz1
      zz2=cmat(3,1)*xx1+cmat(3,2)*yy1+cmat(3,3)*zz1
      x2=sngl(xx2)
      y2=sngl(yy2)
      z2=sngl(zz2)
c
c      write(0,*)'cotr2: ',xx2,yy2,zz2
c
      return
      end
c----------------------------------------------------------------
      subroutine crcol(c,nc,g,ng)
c----------------------------------------------------------------
      character*(*) c
      real g(*)
      character*256 r
      character*256 r2
c      write(0,*)'crcol: ',c
      r=c
      ng=0
100   continue
      k1=0
      k2=-1
      do 110 i=1,nc
      if(r(i:i).ne.' ') then
      if(k1.eq.0) k1=i
      k2=i
      endif
      if(r(i:i).eq.':') goto 111
110   continue
111   continue
c      write(0,*)' k1,k2 ',k1,k2
      if(k2.lt.k1) return
      if(r(k1:k1).eq.':')r(k1:k1)=' '
      if(r(k2:k2).eq.':')r(k2:k2)=' '
      ng=ng+1
      r2=' '
      r2=r(k1:k2)
c      write(0,*)'crcol r2: ',r(k1:k2)
      g(ng)=gets(r2)
      do 120 i=k1,k2
      r(i:i)=' '
120   continue
      goto 100
      end
c-----------------------------------------------------------
      subroutine skipii(iu,istep,isdim)
c-----------------------------------------------------------
      character*80 rec
      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,err=190) rec
      if(rec(1:10).eq.'FIELD-2D-1') then
      read(iu,1000,end=190,err=190) rec
      if((isany.eq.0).and.l1(1:m).ne.rec(1:m)) goto 100
      l1=rec
      read(iu,1000,end=190,err=190) l2
      read(iu,*,end=190,err=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.---------------------------------------
c character encoding of real array into 2 bytes
      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
      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.lt.0)ir=ir+256
      if(ir.eq.32) goto 190
      if(ir.gt.127) then
      ic=ir-170
      j=j+1
      ir=ichar(c(j:j))
      if(ir.lt.0)ir=ir+256
      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.---------------------------------------
      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
      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------------------------------------------
      real function gets(cc)
c------------------------------------------
c
c  decodes integer or floating f format
c  from character string 
c
      character*(*) cc
      nn=len(cc)
c
      gets=0.0
      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,'$'
      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 ipol2g(a,nx,ny,b,gx,gy,x,y,iout)
c.--------------------------------------------------------------
      real a(nx,ny),gx(nx),gy(ny)
      data ixl,iyl/1,1/ 
c
c
      iout=0
      ix=ixl
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94969 
      ix=min0(ixl+1,nx-1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94969 
      ix=max0(ixl-1,1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94969 
      do 94961 i=1,nx-1
      if(x.ge.gx(i).and.x.le.gx(i+1)) then
      ix=i
      goto 94969 
      endif
94961 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94969 continue
      ixl=ix
      iy=iyl
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94959 
      iy=min0(iyl+1,ny-1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94959 
      iy=max0(iyl-1,1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94959 
      do 94951 i=1,ny-1
      if(y.ge.gy(i).and.y.le.gy(i+1)) then
      iy=i
      goto 94959 
      endif
94951 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94959 continue
      iyl=iy
c
      dxi=(x-gx(ix))/(gx(ix+1)-gx(ix))
      dyi=(y-gy(iy))/(gy(iy+1)-gy(iy))
      x00=a(ix,iy)+dxi*(a(ix+1,iy)-a(ix,iy))
      x10=a(ix,iy+1)+dxi*(a(ix+1,iy+1)-a(ix,iy+1))
      b=x00+dyi*(x10-x00)
      return
      end
c.--------------------------------------------------------------
      subroutine ipol2n(a,nx,ny,b,gx,gy,x,y,iout)
c.--------------------------------------------------------------
      real a(nx,ny),gx(nx),gy(ny)
      data ixl,iyl/1,1/ 
c
c
      iout=0
      ix=ixl
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94949 
      ix=min0(ixl+1,nx-1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94949 
      ix=max0(ixl-1,1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94949 
      do 94941 i=1,nx-1
      if(x.ge.gx(i).and.x.le.gx(i+1)) then
      ix=i
      goto 94949 
      endif
94941 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94949 continue
      ixl=ix
      iy=iyl
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94939 
      iy=min0(iyl+1,ny-1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94939 
      iy=max0(iyl-1,1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94939 
      do 94931 i=1,ny-1
      if(y.ge.gy(i).and.y.le.gy(i+1)) then
      iy=i
      goto 94939 
      endif
94931 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94939 continue
      iyl=iy
c
      kx=ix
      if( abs(x-gx(ix)) .gt. abs(x-gx(ix+1)) )kx=ix+1
      ky=iy
      if( abs(y-gy(iy)) .gt. abs(y-gy(iy+1)) )ky=iy+1
      b=a(kx,ky)
      return
      end
c--------------------------------------------------------------
      subroutine ipol3a(a,nx,ny,nz,b,gx,gy,gz,x,y,z,iout)
c--------------------------------------------------------------
      real a(nx,ny,nz),gx(nx),gy(ny),gz(nz)
      data ixl,iyl,izl /1,1,1/ 
c
c
      iout=0
      ix=ixl
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94929 
      ix=min0(ixl+1,nx-1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94929 
      ix=max0(ixl-1,1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94929 
      do 94921 i=1,nx-1
      if(x.ge.gx(i).and.x.le.gx(i+1)) then
      ix=i
      goto 94929 
      endif
94921 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94929 continue
      ixl=ix
      iy=iyl
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94919 
      iy=min0(iyl+1,ny-1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94919 
      iy=max0(iyl-1,1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94919 
      do 94911 i=1,ny-1
      if(y.ge.gy(i).and.y.le.gy(i+1)) then
      iy=i
      goto 94919 
      endif
94911 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94919 continue
      iyl=iy
      iz=izl
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94909 
      iz=min0(izl+1,nz-1)
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94909 
      iz=max0(izl-1,1)
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94909 
      do 94901 i=1,nz-1
      if(z.ge.gz(i).and.z.le.gz(i+1)) then
      iz=i
      goto 94909 
      endif
94901 continue
      iout=1
c      write(0,*)' ipol2g miss ',gA(1),A,gA(nA),iAl
      return
94909 continue
      izl=iz
c
      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 ipol3g(a,nx,ny,nz,b,gx,gy,gz,x,y,z,iout)
c.--------------------------------------------------------------
      real a(nx,ny,nz),gx(nx),gy(ny),gz(nz)
      data ixl,iyl,izl /1,1,1/ 
c
c
      iout=0
      ix=ixl
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94899 
      ix=min0(ixl+1,nx-1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94899 
      ix=max0(ixl-1,1)
      if(x.ge.gx(ix).and.x.le.gx(ix+1)) goto 94899 
      do 94891 i=1,nx-1
      if(x.ge.gx(i).and.x.le.gx(i+1)) then
      ix=i
      goto 94899 
      endif
94891 continue
      iout=1
      return
94899 continue
      ixl=ix
      iy=iyl
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94889 
      iy=min0(iyl+1,ny-1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94889 
      iy=max0(iyl-1,1)
      if(y.ge.gy(iy).and.y.le.gy(iy+1)) goto 94889 
      do 94881 i=1,ny-1
      if(y.ge.gy(i).and.y.le.gy(i+1)) then
      iy=i
      goto 94889 
      endif
94881 continue
      iout=1
      return
94889 continue
      iyl=iy
      iz=izl
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94879 
      iz=min0(izl+1,nz-1)
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94879 
      iz=max0(izl-1,1)
      if(z.ge.gz(iz).and.z.le.gz(iz+1)) goto 94879 
      do 94871 i=1,nz-1
      if(z.ge.gz(i).and.z.le.gz(i+1)) then
      iz=i
      goto 94879 
      endif
94871 continue
      iout=1
      return
94879 continue
      izl=iz
c
      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-------------------------------------
      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 lineq1(a,n,b,rc)
c-------------------------------------------------------------
c
c...... simple interface to solve lineae equations
c
      real a(n,n),b(n)
      common /dfdfa/ip(10000)
      common /dccfa/z(10000)
      call sgeco(a,n,n,ip,rc,z)
      job=0
      call sgesl(a,n,n,ip,b,job)
      return 
      end
c-------------------------------------------------------------
      subroutine sgeco(a,lda,n,ipvt,rcond,z)
c-------------------------------------------------------------
      integer lda,n,ipvt(1)
      real a(lda,1),z(1)
      real rcond
c
c     sgeco factors a real matrix by gaussian elimination
c     and estimates the condition of the matrix.
c
c     if  rcond  is not needed, sgefa is slightly faster.
c     to solve  a*x = b , follow sgeco by sgesl.
c     to compute  inverse(a)*c , follow sgeco by sgesl.
c     to compute  determinant(a) , follow sgeco by sgedi.
c     to compute  inverse(a) , follow sgeco by sgedi.
c
c     on entry
c
c        a       real(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        rcond   real
c                an estimate of the reciprocal condition of  a .
c                for the system  a*x = b , relative perturbations
c                in  a  and  b  of size  epsilon  may cause
c                relative perturbations in  x  of size  epsilon/rcond .
c                if  rcond  is so small that the logical expression
c                           1.0 + rcond .eq. 1.0
c                is true, then  a  may be singular to working
c                precision.  in particular,  rcond  is zero  if
c                exact singularity is detected or the estimate
c                underflows.
c
c        z       real(n)
c                a work vector whose contents are usually unimportant.
c                if  a  is close to a singular matrix, then  z  is
c                an approximate null vector in the sense that
c                norm(a*z) = rcond*norm(a)*norm(z) .
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     linpack sgefa
c     blas saxpy,sdot,sscal,sasum
c     fortran abs,amax1,sign
c
c     internal variables
c
      real sdot,ek,t,wk,wkm
      real anorm,s,sasum,sm,ynorm
      integer info,j,k,kb,kp1,l
c
c
c     compute 1-norm of a
c
      anorm = 0.0e0
      do 10 j = 1, n
      anorm = amax1(anorm,sasum(n,a(1,j),1))
10    continue
c
c     factor
c
      call sgefa(a,lda,n,ipvt,info)
c
c     rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
c     estimate = norm(z)/norm(y) where  a*z = y  and  trans(a)*y = e .
c     trans(a)  is the transpose of a .  the components of  e  are
c     chosen to cause maximum local growth in the elements of w  where
c     trans(u)*w = e .  the vectors are frequently rescaled to avoid
c     overflow.
c
c     solve trans(u)*w = e
c
      ek = 1.0e0
      do 20 j = 1, n
      z(j) = 0.0e0
20    continue
      do 100 k = 1, n
      if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k))
      if (abs(ek-z(k)) .le. abs(a(k,k))) go to 30
      s = abs(a(k,k))/abs(ek-z(k))
      call sscal(n,s,z,1)
      ek = s*ek
30    continue
      wk = ek - z(k)
      wkm = -ek - z(k)
      s = abs(wk)
      sm = abs(wkm)
      if (a(k,k) .eq. 0.0e0) go to 40
      wk = wk/a(k,k)
      wkm = wkm/a(k,k)
      go to 50
40    continue
      wk = 1.0e0
      wkm = 1.0e0
50    continue
      kp1 = k + 1
      if (kp1 .gt. n) go to 90
      do 60 j = kp1, n
      sm = sm + abs(z(j)+wkm*a(k,j))
      z(j) = z(j) + wk*a(k,j)
      s = s + abs(z(j))
60    continue
      if (s .ge. sm) go to 80
      t = wkm - wk
      wk = wkm
      do 70 j = kp1, n
      z(j) = z(j) + t*a(k,j)
70    continue
80    continue
90    continue
      z(k) = wk
100   continue
      s = 1.0e0/sasum(n,z,1)
      call sscal(n,s,z,1)
c
c     solve trans(l)*y = w
c
      do 120 kb = 1, n
      k = n + 1 - kb
      if (k .lt. n) z(k) = z(k) + sdot(n-k,a(k+1,k),1,z(k+1),1)
      if (abs(z(k)) .le. 1.0e0) go to 110
      s = 1.0e0/abs(z(k))
      call sscal(n,s,z,1)
110   continue
      l = ipvt(k)
      t = z(l)
      z(l) = z(k)
      z(k) = t
120   continue
      s = 1.0e0/sasum(n,z,1)
      call sscal(n,s,z,1)
c
      ynorm = 1.0e0
c
c     solve l*v = y
c
      do 140 k = 1, n
      l = ipvt(k)
      t = z(l)
      z(l) = z(k)
      z(k) = t
      if (k .lt. n) call saxpy(n-k,t,a(k+1,k),1,z(k+1),1)
      if (abs(z(k)) .le. 1.0e0) go to 130
      s = 1.0e0/abs(z(k))
      call sscal(n,s,z,1)
      ynorm = s*ynorm
130   continue
140   continue
      s = 1.0e0/sasum(n,z,1)
      call sscal(n,s,z,1)
      ynorm = s*ynorm
c
c     solve  u*z = v
c
      do 160 kb = 1, n
      k = n + 1 - kb
      if (abs(z(k)) .le. abs(a(k,k))) go to 150
      s = abs(a(k,k))/abs(z(k))
      call sscal(n,s,z,1)
      ynorm = s*ynorm
150   continue
      if (a(k,k) .ne. 0.0e0) z(k) = z(k)/a(k,k)
      if (a(k,k) .eq. 0.0e0) z(k) = 1.0e0
      t = -z(k)
      call saxpy(k-1,t,a(1,k),1,z(1),1)
160   continue
c     make znorm = 1.0
      s = 1.0e0/sasum(n,z,1)
      call sscal(n,s,z,1)
      ynorm = s*ynorm
c
      if (anorm .ne. 0.0e0) rcond = ynorm/anorm
      if (anorm .eq. 0.0e0) rcond = 0.0e0
      return
      end
      integer function isamax(n,sx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c     modified 3/93 to return if incx .le. 0.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(*),smax
      integer i,incx,ix,n
c
      isamax = 0
      if( n.lt.1 .or. incx.le.0 ) return
      isamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      smax = abs(sx(1))
      ix = ix + incx
      do 10 i = 2,n
      if(abs(sx(ix)).le.smax) go to 5
      isamax = i
      smax = abs(sx(ix))
5     ix = ix + incx
10    continue
      return
c
c        code for increment equal to 1
c
20    smax = abs(sx(1))
      do 30 i = 2,n
      if(abs(sx(i)).le.smax) go to 30
      isamax = i
      smax = abs(sx(i))
30    continue
      return
      end
      real function sasum(n,sx,incx)
c
c     takes the sum of the absolute values.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 3/93 to return if incx .le. 0.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(*),stemp
      integer i,incx,m,mp1,n,nincx
c
      sasum = 0.0e0
      stemp = 0.0e0
      if( n.le.0 .or. incx.le.0 )return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
      stemp = stemp + abs(sx(i))
10    continue
      sasum = stemp
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
20    m = mod(n,6)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
      stemp = stemp + abs(sx(i))
30    continue
      if( n .lt. 6 ) go to 60
40    mp1 = m + 1
      do 50 i = mp1,n,6
      stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2)) + 
     *abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5))
50    continue
60    sasum = stemp
      return
      end
      real function sdot(n,sx,incx,sy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(*),sy(*),stemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      stemp = 0.0e0
      sdot = 0.0e0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
      stemp = stemp + sx(ix)*sy(iy)
      ix = ix + incx
      iy = iy + incy
10    continue
      sdot = stemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
20    m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
      stemp = stemp + sx(i)*sy(i)
30    continue
      if( n .lt. 5 ) go to 60
40    mp1 = m + 1
      do 50 i = mp1,n,5
      stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + sx(i + 2)*
     *sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
50    continue
60    sdot = stemp
      return
      end
      subroutine sscal(n,sa,sx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to 1.
c     jack dongarra, linpack, 3/11/78.
c     modified 3/93 to return if incx .le. 0.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sa,sx(*)
      integer i,incx,m,mp1,n,nincx
c
      if( n.le.0 .or. incx.le.0 )return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
      sx(i) = sa*sx(i)
10    continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
20    m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
      sx(i) = sa*sx(i)
30    continue
      if( n .lt. 5 ) return
40    mp1 = m + 1
      do 50 i = mp1,n,5
      sx(i) = sa*sx(i)
      sx(i + 1) = sa*sx(i + 1)
      sx(i + 2) = sa*sx(i + 2)
      sx(i + 3) = sa*sx(i + 3)
      sx(i + 4) = sa*sx(i + 4)
50    continue
      return
      end
      subroutine sgefa(a,lda,n,ipvt,info)
      integer lda,n,ipvt(1),info
      real a(lda,1)
c
c     sgefa factors a real matrix by gaussian elimination.
c
c     sgefa is usually called by sgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for sgeco) = (1 + 9/n)*(time for sgefa) .
c
c     on entry
c
c        a       real(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that sgesl or sgedi will divide by zero
c                     if called.  use  rcond  in sgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas saxpy,sscal,isamax
c
c     internal variables
c
      real t
      integer isamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
      kp1 = k + 1
c
c        find l = pivot index
c
      l = isamax(n-k+1,a(k,k),1) + k - 1
      ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
      if (a(l,k) .eq. 0.0e0) go to 40
c
c           interchange if necessary
c
      if (l .eq. k) go to 10
      t = a(l,k)
      a(l,k) = a(k,k)
      a(k,k) = t
10    continue
c
c           compute multipliers
c
      t = -1.0e0/a(k,k)
      call sscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
      do 30 j = kp1, n
      t = a(l,j)
      if (l .eq. k) go to 20
      a(l,j) = a(k,j)
      a(k,j) = t
20    continue
      call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
30    continue
      go to 50
40    continue
      info = k
50    continue
60    continue
70    continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0e0) info = n
      return
      end
      subroutine sgedi(a,lda,n,ipvt,det,work,job)
      integer lda,n,ipvt(1),job
      real a(lda,1),det(2),work(1)
c
c     sgedi computes the determinant and inverse of a matrix
c     using the factors computed by sgeco or sgefa.
c
c     on entry
c
c        a       real(lda, n)
c                the output from sgeco or sgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from sgeco or sgefa.
c
c        work    real(n)
c                work vector.  contents destroyed.
c
c        job     integer
c                = 11   both determinant and inverse.
c                = 01   inverse only.
c                = 10   determinant only.
c
c     on return
c
c        a       inverse of original matrix if requested.
c                otherwise unchanged.
c
c        det     real(2)
c                determinant of original matrix if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. abs(det(1)) .lt. 10.0
c                or  det(1) .eq. 0.0 .
c
c     error condition
c
c        a division by zero will occur if the input factor contains
c        a zero on the diagonal and the inverse is requested.
c        it will not occur if the subroutines are called correctly
c        and if sgeco has set rcond .gt. 0.0 or sgefa has set
c        info .eq. 0 .
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas saxpy,sscal,sswap
c     fortran abs,mod
c
c     internal variables
c
      real t
      real ten
      integer i,j,k,kb,kp1,l,nm1
c
c
c     compute determinant
c
      if (job/10 .eq. 0) go to 70
      det(1) = 1.0e0
      det(2) = 0.0e0
      ten = 10.0e0
      do 50 i = 1, n
      if (ipvt(i) .ne. i) det(1) = -det(1)
      det(1) = a(i,i)*det(1)
c        ...exit
      if (det(1) .eq. 0.0e0) go to 60
10    if (abs(det(1)) .ge. 1.0e0) go to 20
      det(1) = ten*det(1)
      det(2) = det(2) - 1.0e0
      go to 10
20    continue
30    if (abs(det(1)) .lt. ten) go to 40
      det(1) = det(1)/ten
      det(2) = det(2) + 1.0e0
      go to 30
40    continue
50    continue
60    continue
70    continue
c
c     compute inverse(u)
c
      if (mod(job,10) .eq. 0) go to 150
      do 100 k = 1, n
      a(k,k) = 1.0e0/a(k,k)
      t = -a(k,k)
      call sscal(k-1,t,a(1,k),1)
      kp1 = k + 1
      if (n .lt. kp1) go to 90
      do 80 j = kp1, n
      t = a(k,j)
      a(k,j) = 0.0e0
      call saxpy(k,t,a(1,k),1,a(1,j),1)
80    continue
90    continue
100   continue
c
c        form inverse(u)*inverse(l)
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 140
      do 130 kb = 1, nm1
      k = n - kb
      kp1 = k + 1
      do 110 i = kp1, n
      work(i) = a(i,k)
      a(i,k) = 0.0e0
110   continue
      do 120 j = kp1, n
      t = work(j)
      call saxpy(n,t,a(1,j),1,a(1,k),1)
120   continue
      l = ipvt(k)
      if (l .ne. k) call sswap(n,a(1,k),1,a(1,l),1)
130   continue
140   continue
150   continue
      return
      end
      subroutine saxpy(n,sa,sx,incx,sy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(*),sy(*),sa
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
      sy(iy) = sy(iy) + sa*sx(ix)
      ix = ix + incx
      iy = iy + incy
10    continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
20    m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
      sy(i) = sy(i) + sa*sx(i)
30    continue
      if( n .lt. 4 ) return
40    mp1 = m + 1
      do 50 i = mp1,n,4
      sy(i) = sy(i) + sa*sx(i)
      sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
      sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
      sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
50    continue
      return
      end
      subroutine sswap (n,sx,incx,sy,incy)
c
c     interchanges two vectors.
c     uses unrolled loops for increments equal to 1.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(*),sy(*),stemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c       code for unequal increments or equal increments not equal
c         to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
      stemp = sx(ix)
      sx(ix) = sy(iy)
      sy(iy) = stemp
      ix = ix + incx
      iy = iy + incy
10    continue
      return
c
c       code for both increments equal to 1
c
c
c       clean-up loop
c
20    m = mod(n,3)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
      stemp = sx(i)
      sx(i) = sy(i)
      sy(i) = stemp
30    continue
      if( n .lt. 3 ) return
40    mp1 = m + 1
      do 50 i = mp1,n,3
      stemp = sx(i)
      sx(i) = sy(i)
      sy(i) = stemp
      stemp = sx(i + 1)
      sx(i + 1) = sy(i + 1)
      sy(i + 1) = stemp
      stemp = sx(i + 2)
      sx(i + 2) = sy(i + 2)
      sy(i + 2) = stemp
50    continue
      return
      end
c-------------------------------------------------------
      subroutine sgesl(a,lda,n,ipvt,b,job)
c-------------------------------------------------------
      integer lda,n,ipvt(1),job
      real a(lda,1),b(1)
c
c     sgesl solves the real system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by sgeco or sgefa.
c
c     on entry
c
c        a       real(lda, n)
c                the output from sgeco or sgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from sgeco or sgefa.
c
c        b       real(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if sgeco has set rcond .gt. 0.0
c        or sgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call sgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call sgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas saxpy,sdot
c
c     internal variables
c
      real sdot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
      if (nm1 .lt. 1) go to 30
      do 20 k = 1, nm1
      l = ipvt(k)
      t = b(l)
      if (l .eq. k) go to 10
      b(l) = b(k)
      b(k) = t
10    continue
      call saxpy(n-k,t,a(k+1,k),1,b(k+1),1)
20    continue
30    continue
c
c        now solve  u*x = y
c
      do 40 kb = 1, n
      k = n + 1 - kb
      b(k) = b(k)/a(k,k)
      t = -b(k)
      call saxpy(k-1,t,a(1,k),1,b(1),1)
40    continue
      go to 100
50    continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
      do 60 k = 1, n
      t = sdot(k-1,a(1,k),1,b(1),1)
      b(k) = (b(k) - t)/a(k,k)
60    continue
c
c        now solve trans(l)*x = y
c
      if (nm1 .lt. 1) go to 90
      do 80 kb = 1, nm1
      k = n - kb
      b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1)
      l = ipvt(k)
      if (l .eq. k) go to 70
      t = b(l)
      b(l) = b(k)
      b(k) = t
70    continue
80    continue
90    continue
100   continue
      return
      end
c--------------------------------------------------
      subroutine mcmap(rec,ncmap,ir,ig,ib)
c--------------------------------------------------
      character*80 rec
      integer ir(*),ig(*),ib(*)
c
      ic1=1
      ic2=ncmap
      ich=ncmap/2
c
      if(rec(1:4).eq.'lin5') then
      do 439 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      hue1=1.75
      hue2=0.95
      sat1=0.99
      sat2=0.99
      val1=0.99
      val2=0.99
      h=hue1+t*(hue2-hue1)
      if(h.gt.1.0)h=h-1.0
      s=sat1+t*(sat2-sat1)
      v=val1+t*(val2-val1)
      v=v*0.5*(1.6+0.4*sin(16.0*6.28*t))
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
439   continue
      endif
c
      if(rec(1:4).eq.'lin1') then
      do 409 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      hue1=1.75
      hue2=0.95
      sat1=0.99
      sat2=0.99
      val1=0.99
      val2=0.99
      h=hue1+t*(hue2-hue1)
      if(h.gt.1.0)h=h-1.0
      s=sat1+t*(sat2-sat1)
      v=val1+t*(val2-val1)
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
409   continue
      endif
c
      if(rec(1:4).eq.'lin2') then
      do 408 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      hue1=1.75
      hue2=0.95
      sat1=0.05
      sat2=0.99
      val1=0.99
      val2=0.99
      h=hue1+t*(hue2-hue1)
      if(h.gt.1.0)h=h-1.0
      s=sat1+t*(sat2-sat1)
      v=val1+t*(val2-val1)
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
408   continue
      endif
c
      if(rec(1:4).eq.'lin3') then
      do 417 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      hue1=1.75
      hue2=0.95
      sat1=0.05
      sat2=0.99
      val1=0.99
      val2=0.79
      h=hue1+t*(hue2-hue1)
      if(h.gt.1.0)h=h-1.0
      s=sat1+(6.0*t)*(sat2-sat1)
      s=amax1(sat1,amin1(sat2,s))
      v=val1+t*(val2-val1)
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
417   continue
      endif
c
      if(rec(1:4).eq.'sym1') then
      do 410 ic=ic1,ic2
      s=float(iabs(ic-ich))/float(ich+1)
      v=1.0-0.2*s
      if(ic.lt.ich) h=0.5+s*0.2
      if(ic.ge.ich) h=0.15-s*0.15
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
410   continue
      endif
c
      if(rec(1:4).eq.'log1') then
      do 411 ic=ic1,ic2
      s=float(iabs(ic-ic1))/float(ic2-ic1+2)
      v=1.0-0.2*s
      h=0.8+s*0.4
      if(h.ge.1.0)h=h-1.0
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
411   continue
      endif
c
      if(rec(1:4).eq.'log2') then
      do 412 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      v=1.0-0.5*s
      h=0.5+t*0.6
      if(h.ge.1.0)h=h-1.0
      s=0.1+0.9*t
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
412   continue
      endif
c
      if(rec(1:4).eq.'log3') then
      do 413 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      v=1.0-0.2*t
      v=1.0-0.4*t*t
      h=1.75-0.8*t
      if(h.ge.1.0)h=h-1.0
      s=0.05+0.95*t
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
413   continue
      endif
c
c
      if(rec(1:4).eq.'log4') then
      do 414 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      v=1.0
      h=0.1+1.00*t
      if(h.gt.1.0)h=h-1.0 
      h=1.0-h
      s=0.1+0.4*t+0.60*t*(sin(t*0.5*pi*7.2)**2)
      s=amax1(0.0,amin1(1.0,s))
      call hsvrgb(h,s,v,rr,gg,bb)
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
414   continue
      endif
c
      if(rec(1:5).eq.'logbw') then
      do 415 ic=ic1,ic2
      t=float(iabs(ic-ic1))/float(ic2-ic1+2)
      rr=1.0-t*t
      gg=1.0-t*t
      bb=1.0-t*t
      jjr=int(rr*255.0)
      ir(ic)=max0(1,min0(255,jjr)) 
      jjg=int(gg*255.0)
      ig(ic)=max0(1,min0(255,jjg)) 
      jjb=int(bb*255.0)
      ib(ic)=max0(1,min0(255,jjb)) 
415   continue
      endif
c
      lrec=max(jlen(rec,80),5)
      if(rec(lrec-4:lrec).eq.'noneg')then
      ir(1)=255
      ig(1)=255
      ib(1)=255
      endif
c
c      do 500 i=1,ncmap
c      write(0,*)i,ir(i),ig(i),ib(i)
c500   continue
c
      return
      end
c-----------------------------------------------------
      subroutine minmax(x,n,xmin,xmax)
c-----------------------------------------------------
c
c  utility routine to compute minimum and maximum
c  of a real array
c
      real x(n)
      xmin=x(1)
      xmax=x(1)
      do 100 i=2,n
      if(x(i).gt.xmax) xmax=x(i)
      if(x(i).lt.xmin) xmin=x(i)
100   continue
      return
      end
c----------------------------------------
      subroutine movlev(a,b,n)
c----------------------------------------
      real a(n),b(n)
      do 100 i=1,n
100   b(i)=a(i)
      return
      end
c----------------------------------------
      subroutine movliv(a,b,n)
c----------------------------------------
      integer a(n),b(n)
      do 100 i=1,n
100   b(i)=a(i)
      return
      end
c--------------------------------------------
      real function ran2(idm)
c--------------------------------------------
c
c   portable random number generator
c   after knuth "numerical recipes"
c
      parameter(m=714027,ia=1366,ic=150889,rm=1.400507264e-6)
      dimension ir(0:100)
      save
      data idum /-1123456/
      if(idum.lt.0) then
      idum=mod(ic-idum,m)
      do 11 j=1,97
      idum=mod(ia*idum+ic,m)
      ir(j)=idum
11    continue
      idum=mod(ia*idum+ic,m)
      iy=idum
      endif
      j=1+(97*iy)/m
      iy=ir(j)
      ran2=iy*rm
      idum=mod(ia*idum+ic,m)
      ir(j)=idum
      return
      end
c---------------------------------------------
      integer function mindex(c1,c2)
c---------------------------------------------
c..... find position of c2 in c1
      character*(*) c1,c2
      l1=len(c1)
      l2=len(c2)
      do 100 i=1,l1-l2+1
      do 200 j=1,l2
      ij=i+j-1
      if(c1(ij:ij).ne.c2(j:j)) goto 210
200   continue
      mindex=i
      return
210   continue
100   continue
      mindex=0
      return
      end
c--------------------------------------
      subroutine hsvrgb(h,s,v,r,g,b)
c--------------------------------------
      implicit real (a-z)
c
c   converts hue/saturation/value to red/green/blue
c
      hh=6.*h
      i=aint(hh)
      f=hh-i
      m=v*(1.-s)
      n=v*(1.-s*f)
      k=v*(1.-s*(1.-f))
      goto (1,2,3,4,5) int(i)
      r=v
      g=k
      b=m
      return
1     r=n
      g=v
      b=m
      return
2     r=m
      g=v
      b=k
      return
3     r=m
      g=n
      b=v
      return
4     r=k
      g=m
      b=v
      return
5     r=v
      g=m
      b=n
      return
      end
c---------------------------------------------------
      subroutine timstr(r,iy,mo,id,ih,mi,se)
c---------------------------------------------------
      character*(*) r
      character*80 rec
      character*3 cmo(12)
      data cmo /'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP'
     *,'OCT','NOV','DEC'/
      imo=0
      jd=julday(iy,mo,id)
      jy=iy
      if(jy.lt.150)jy=jy+1900
      if(imo.eq.0) then
      ise1=se
      ise2=(1000.0*se)
      ise2=mod(ise2,1000)
      write(rec,1000)jy,':',mo,':',id,':',ih,':',mi,':',ise1,'.',ise2
1000  format(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i3.3)
      else if(imo.eq.1) then
      write(rec,1001)jy,':',cmo(mo),'(',mo,'):',id,'(',jd,'):',ih,':'
     *,mi,':',se
1001  format(i4.4,a,a,a,i2.2,a,i2.2,a,i3.3,a,i2.2,a,i2.2,a,f6.3)
      endif
      r=rec
      return
      end
c-----------------------------------------------
      subroutine print_t1966(bt,cc)
c-----------------------------------------------
      character*(*) cc
      real*8 bt
      character*80 rec
      call epoch1966(bt,iy,mo,id,ih,mi,se,1)
      call timstr(rec,iy,mo,id,ih,mi,se)
      write(0,'(a,a,a)')cc,' ',rec(1:jlen(rec,80))
      return
      end
c-----------------------------------------------
      subroutine get_steptime1966(bt)
c-----------------------------------------------
      real*8 bt
      character*80 rec
      istep=-9999
      call xparse('-step',1,2,'required',rec,istep,rdum)
      if(istep.eq.-9999) then
      write(0,*)'get_steptime:no step on cammand line'
      stop
      endif
      call get_basetime1966(bt)
      bt=bt+dble(float(istep))
      return
      end
c-----------------------------------------------
      subroutine get_basetime1966(bt)
c-----------------------------------------------
      real*8 bt
      real a(10)
      character*80 rec
c
c.... get basetime from command line
c
      rec='none'
      call xparse('-basetime',1,1,'required',rec,idum,rdum)
      if(rec.eq.'none') then
      write(0,*)'get_basetime:no basetime on cammand line'
      stop
      endif
      call crcol(rec,80,a,nb)
      iy=a(1)+0.1
      mo=a(2)+0.1
      id=a(3)+0.1
      ih=a(4)+0.1
      mi=a(5)+0.1
      se=a(6)
      call epoch1966(bt,iy,mo,id,ih,mi,se,0)
      return
      end
c------------------------------------------------------------------
      subroutine trace2(nx,ny,nz,gx,gy,gz,bx,by,bz,x1,y1,z1,x2,y2,z2,
     *dir)
c------------------------------------------------------------------
      real gx(*)
      real gy(*)
      real gz(*)
      real bx(nx,ny,nz)
      real by(nx,ny,nz)
      real bz(nx,ny,nz)
c..... trace one step
c
      r1=sqrt(x1*x1+y1*y1+z1*z1)
      if(r1.gt.3.7) then
      call ipol3a(bx,nx,ny,nz,bbx,gx,gy,gz,x1,y1,z1,iout)
      if(iout.ne.0) goto 999
      call ipol3a(by,nx,ny,nz,bby,gx,gy,gz,x1,y1,z1,iout)
      if(iout.ne.0) goto 999
      call ipol3a(bz,nx,ny,nz,bbz,gx,gy,gz,x1,y1,z1,iout)
      if(iout.ne.0) goto 999
      else
      call cotr('gse','sm ',x1,y1,z1,xsm,ysm,zsm)
      xbbx=-3.0*xsm*zsm
      xbby=-3.0*ysm*zsm
      xbbz=r1*r1-3.0*zsm*zsm
      call cotr('sm ','gse',xbbx,xbby,xbbz,bbx,bby,bbz)
      endif
c      write(0,'(a,3(1x,g12.5))')' trace2: ',bbx,bby,bbz
c
      bb=sqrt(bbx*bbx+bby*bby+bbz*bbz)
      if(bb.le.0.0) goto 999
      s=dir/bb
      x2=x1+bbx*s
      y2=y1+bby*s
      z2=z1+bbz*s
c
      return
999   continue
      x2=1.e6
      y2=1.e6
      z2=1.e6
      return
      end
c---------------------------------------------
      subroutine wform(rec,c,nn)
c---------------------------------------------
      character*(*) rec
      character r*16
c
      cc=abs(c)
c
      if(cc.ge.100.and.cc.lt.10000000.) then
      ic=cc
      write(r,'(i8)') ic
c
      else if(cc.ge.10..and.cc.lt.100.) then
      write(r,'(f7.1)') cc
      else if(cc.ge.1..and.cc.lt.10.) then
      write(r,'(f7.2)') cc
      else if(cc.ge.0.1.and.cc.lt.1.) then
      write(r,'(f7.3)') cc
      else if(cc.ge.0.01.and.cc.lt.0.1) then
      write(r,'(f8.4)') cc
      else if(cc.ge.0.001.and.cc.lt.0.01) then
      write(r,'(f9.5)') cc
      else if(cc.ge.0.0001.and.cc.lt.0.001) then
      write(r,'(f12.6)') cc
      else if(cc.ge.0.00001.and.cc.lt.0.0001) then
      write(r,'(f12.7)') cc
c
      else if(cc.eq.0.) then
      r='0'
c
      else
      write(r,'(1pe12.3)') cc
      endif
c
      ie=0
      ip=0
      n=0
      do 100 i=1,15
      if(r(i:i).eq.'e'.or.r(i:i).eq.'E') ie=i
      if(r(i:i).eq.'.') ip=i
      if(r(i:i).ne.' ') n=i
100   continue
c------  remove trailing zeroes in mantissa
      if(ip.gt.0) then
      m=n
      if(ie.ne.0) m=ie-1
      is=0
      do 110 i=m,ip,-1
      if(r(i:i).ne.'0') is=1
      if(is.eq.0) r(i:i)=' '
110   continue
      endif
c-----  remove leading zeroes in exponent
      if(ie.gt.0) then
      is=0
      do 120 i=ie+1,n
      if(r(i:i).eq.'+') r(i:i)=' '
      ic=ichar(r(i:i))
      ic1=ichar('1')
      ic9=ichar('9')
      if(ic.ge.ic1.and.ic.le.ic9) is=1
      if(is.eq.0.and.r(i:i).eq.'0') r(i:i)=' '
120   continue
      endif
c------  compress
      k=0
      if(c.lt.0.) then
      k=1
      rec(k:k)='-'
      endif
      do 130 i=1,n
      if(r(i:i).ne.' ') then
      k=k+1
      rec(k:k)=r(i:i)
      endif
130   continue
c----- remove trailing point
      if(ie.eq.0.and.rec(k:k).eq.'.') k=k-1
      nn=k
c
      return
      end
c
c------------------------------------------------------------------
      subroutine xparse(cc,nth,ifo,req,carg,iarg,rarg)
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
      stop
      endif
      end
c-------------------------------------------------------------------
      subroutine noblanks(r,n)
c-------------------------------------------------------------------
      character r*(*)
      k=0
      do 100 i=1,n
      if(r(i:i).ne.' ') then
      k=k+1
      if(i.ne.k)r(k:k)=r(i:i)
      endif
100   continue
      do 110 i=k+1,n
      r(i:i)=' '
110   continue
      return
      end
c-------------------------------------------------------------
      subroutine stricat(c,ii,f)
c-------------------------------------------------------------
      character*(*) c,f
      character*32 r
      write(r,f)ii
      call strcat(c,r)
      return
      end
c-------------------------------------------------------------
      subroutine strcat(c,d)
c-------------------------------------------------------------
      character*(*) c,d
      lc=len(c)
      ld=len(d)
      kc=0
      do 100 i=1,lc
      if(c(i:i).ne.' ')kc=i
100   continue
      kd=0
      do 110 i=1,ld
      if(d(i:i).ne.' ')kd=i
110   continue
      do 120 i=1,kd
      kc=kc+1
      if(kc.le.lc) c(kc:kc)=d(i:i)
120   continue
      return
      end
c-------------------------------------------------------------
      subroutine strcatb(c,d)
c-------------------------------------------------------------
c.... puts a separating blank in
      character*(*) c,d
      lc=len(c)
      ld=len(d)
      kc=0
      do 100 i=1,lc
      if(c(i:i).ne.' ')kc=i
100   continue
      kc=kc+1
      kd=0
      do 110 i=1,ld
      if(d(i:i).ne.' ')kd=i
110   continue
      do 120 i=1,kd
      kc=kc+1
      if(kc.le.lc) c(kc:kc)=d(i:i)
120   continue
      return
      end
c.include ipol3a.for