program compare
c
c compares the density array from *.res files of two different runs:
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'

      integer ntim   ! # of time points
      integer ntimr  ! # of time points in phi (which may be higher res)

      real tim   ! final time when dt is the same for all ky's.
      real dt1   ! time step when dt is the same for all ky's.
      real time(mz)  ! final time for various ky's.
      real dt(mz)  ! time step for various ky's.

      complex density_1(lz,mz,nz,nspecz)  ! from run 1
      complex err_a(lz,mz,nz,nspecz)

      real avgmag,err,avgerr,xxx
      integer l,m,n
      integer ivec(4)
      intrinsic maxloc,sum

      character run1*80, run2*80  !run names
      integer lrun1,lrun2

      integer iargc,icount
      external iargc, getarg

c determine the names of the runs to compare:
      icount=iargc()
      if(icount .ge. 2) then
         call getarg(1,run1)
         call getarg(2,run2)
      else
         write(6,*) 'usage: compare run1 run2'
         stop
      endif
      lrun1=index(run1,' ')-1
      lrun2=index(run2,' ')-1
      write(6,*) 'Comparing ',run1(1:lrun1)//'.res and ',
     >     run2(1:lrun2)//'.res'

c Read data from the first file

      open(unit=7,file=run1(1:lrun1)//'.res',status='old',err=991)
      ntim=0
      ntimr=0
      call wread1(7,ntim,ntimr,tim,dt1,time,dt)
      close(unit=7)

      density_1=density

c Read data from the second file

      open(unit=7,file=run2(1:lrun2)//'.res',status='old',err=992)
      ntim=0
      ntimr=0
      call wread1(7,ntim,ntimr,tim,dt1,time,dt)
      close(unit=7)

      avgmag=0.0
      err=0.0
      do l=1,ldb
         do m=1,md
            do n=1,nd
               avgmag=avgmag+abs(density_1(l,m,n,1))
               err=max(err,abs(density_1(l,m,n,1)-density(l,m,n,1)))
            enddo
         enddo
      enddo
      avgmag=avgmag/(ldb*md*nd*2)
      write(6,*) 'err, relerr, avgmag=',err,err/avgmag,avgmag

      err=0.0
      err_a=abs(density_1-density)
      ivec=maxloc(abs(err_a))
      err=err_a(ivec(1),ivec(2),ivec(3),ivec(4))

      write(6,*) 'Array op maxerr=',err
      write(6,*) 'Occurs at array position = ',ivec

      avgerr=sum(err_a)/(ldb*md*nd)
      write(6,*) 'Avg abs(err) =',avgerr

      stop

 991  write(6,*) 'Error: Unable to open ', run1(1:lrun1)//'.res'
      stop

 992  write(6,*) 'Error: Unable to open ', run2(1:lrun2)//'.res'

      stop
      end