       PROGRAM PLOT_LIAP
       implicit real*4 (a-h,o-z)
       real*4 yplot(10000),expl(10000)

       character*5 runlbl
       character*7 modellbl
       character*10 enrgylbl
       character*25 filename
       
       print*,' Input model number '
       read*,imod

       if(imod.eq.1)modellbl = 'model01'
       if(imod.eq.2)modellbl = 'model02'
       if(imod.eq.3)modellbl = 'model03'
       if(imod.eq.6)modellbl = 'model06'
       if(imod.eq.10)modellbl = 'model10'
       if(imod.eq.11)modellbl = 'model11'
       if(imod.eq.12)modellbl = 'model12'
       if(imod.eq.20)modellbl = 'model20'
       if(imod.eq.21)modellbl = 'model21'
       if(imod.eq.22)modellbl = 'model22'
       if(imod.eq.30)modellbl = 'model30'
       if(imod.eq.40)modellbl = 'model40'
       if(imod.eq.50)modellbl = 'model50'
       if(imod.eq.51)modellbl = 'model51'
       if(imod.eq.60)modellbl = 'model60'
       if(imod.eq.61)modellbl = 'model61'
       if(imod.eq.62)modellbl = 'model62'
       if(imod.eq.63)modellbl = 'model63'
       if(imod.eq.64)modellbl = 'model64'

       l = 0

       print*,' Input range of energy levels: lower upper '
       read*,ielo,ieup

       ymin = 100.0
       ymax = -100.0
       ysm = 0.0
       ysq = 0.0

       do i = ielo,ieup

          if(i.eq.1)then
               runlbl = 'run01'
               iendrun = 205
          endif

          if(i.eq.2)then
               runlbl = 'run02'
               iendrun = 309
          endif

          if(i.eq.3)then
               runlbl = 'run03'
               iendrun = 314
          endif

          if(i.eq.4)then
               runlbl = 'run04'
               iendrun = 319
          endif

          if(i.eq.5)then
               runlbl = 'run05'
               iendrun = 566
          endif

          if(i.eq.6)then
               runlbl = 'run06'
               iendrun = 1039
          endif

          if(i.eq.7)then
               runlbl = 'run07'
               iendrun = 1142
          endif

          if(i.eq.8)then
               runlbl = 'run08'
               iendrun = 631
          endif

          filename = modellbl//runlbl//'liap_init.dat'
   
          open(99,file=filename,status='unknown')

          do k = 1,7
            read(99,*)
          enddo

          do j = 1,10000

            read(99,*,err=3)norbit,x,z,xdot,ydot,zdot,expl(j)
            l = l+1
            yplot(l) = expl(j)

            if(yplot(l).lt.ymin)ymin = yplot(l)
            if(yplot(l).gt.ymax)ymax = yplot(l)
            ysm = ysm + yplot(l)
            ysq = ysq + yplot(l)*yplot(l)

            print*,l,yplot(l)

          enddo

3         close(99)

       enddo

       ymean = ysm/float(l)
       ydisp = sqrt(ysq/float(l) - ymean*ymean)

       print*,' Total number of orbits = ',l
       print*,' Mean Liap expt = ',ymean
       print*,' Dispersion    = ',ydisp
       print*,' Minimum value  = ',ymin
       print*,' Maximum value  = ',ymax


       call pgbegin(0,'?',1,1)

       print*,' '
 911   print*,' Input min,max values for plots (0,0 to end)'
       read*,ypmin,ypmax
       if(ypmin.eq.0.0.and.ypmax.eq.0.0)goto 999
       print*,' Input max number in bin '
       read*,binmax
       print*,' Input number of bins '
       read*,nbin
       call pgenv(ypmin,ypmax,0.0,binmax,0,0)
       call pghist(l,yplot,ypmin,ypmax,nbin,1)
       call pglabel('Liapunov exponent','Number',modellbl)
       goto 911
 999   call pgend
       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc













