c                                                                     72
      parameter(nxm=6400,nkomp=2)
      
c     implicit none

      real*8 p,pnew,pold,vel,dx,dt,eps,tmax,ampmin,ampmax,pi,time,dtout
      real*8 timeou,dummy,rho,rho1,rho2,length,delta,plam,xrho1,xrho2
      real*8 peakcu,vel1,vel2,xvel1,xvel2,disamp

      integer*4 ntplot,nx,lgraph,nxm,nkomp,ind,i,ikomp,imeth
      integer*4 ntfrog,lfstar,lfinim,accu,lfavin

      character*65 file

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx
      common/graph/ampmin,ampmax,ntplot,lgraph
      common/const/pi,peakcu
      common/method/imeth
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu
      common/dissip/disamp,lag

      open(11,file='wave.err')
c      open(12,file='wave.puls')
      open(13,file='wave.ref')
      open(14,file='wave.tra')

10    write(*,*) ' Graphics (0/1)?'
      read(*,*,err=10) lgraph
      if (lgraph.eq.1) call initialize

      call defaul
      dtout=tmax

100   continue
      write(*,1500) imeth,accu
      write(*,1510) rho1,rho2,xrho1,xrho2
      write(*,1515) vel1,vel2,xvel1,xvel2
      write(*,1520) nx,dx,tmax,dt
      write(*,*)
     &' 0: calculate   1: i.c. 2: length,nx,tmax,lambda 3: Parameters'
      write(*,*) 
     &' 4: Method         5: ntfrog,lfinim 6: ntplot   7: velocity'
      write(*,*) 
     &' 8: translate      9: A -> File    10: Playback'
      write(*,*) 
     &' 11: Mitschreiben 12: Playback->File 13: Playback<-File'
      write(*,*) 
     &' 14: eps,abstma   15: Grenze A und q 16: Default file'
      write(*,*) 
     &' 17: t=0 setzen   18: dissipation    19: fakdt'
      write(*,*) 
     &' 20: Test         21: Ausgabe fuer Plot 22: Spikes'
      write(*,*)
     &' 23: Additive Stoerung 24: Rdbdg  25: iosz,loszwr'

      read(*,*,err=100) ind

      if (ind.eq.9999) then
        stop
      end if

      if (ind.eq.1) then
        call ic 
        call dense
        time=.0d0
      end if

      if (ind.eq.2) then
        length=nx*dx
        plam=dt/dx
        write(*,2000) length,nx,tmax,plam
        read(*,*,err=100) length,nx,tmax,plam
        if (nx.gt.nxm) then
          write(*,*) ' nx zu gross'
          goto 100
        end if
        dx=length/nx
        dt=plam*dx
        call constan
        dtout=tmax
        call dense
      end if

      if (ind.eq.3) then
        write(*,3000) rho1,rho2,xrho1,xrho2
        read(*,*,err=100) rho1,rho2,xrho1,xrho2
        call dense
      end if

      if (ind.eq.4) then
        write(*,4000) imeth
        read(*,*,err=100) imeth
        if (imeth.eq.4.or.imeth.eq.5.or.imeth.eq.6) then
4001      write(*,4010) accu
          read(*,*) accu
          if (accu.ne.2.and.accu.ne.4) goto 4001
        end if
      end if

      if (ind.eq.5) then
        write(*,5000)  ntfrog,lfinim,lfavin,metbot
        read(*,*,err=100)  ntfrog,lfinim,lfavin,metbot
      end if

      if (ind.eq.6) then
        write(*,6000) ntplot
        read(*,*,err=100) ntplot
      end if

      if (ind.eq.7) then
        write(*,7000)  vel1,vel2,xvel1,xvel2
        read(*,*,err=100)  vel1,vel2,xvel1,xvel2
        call dense
      end if

      if (ind.eq.9) then
        write(*,*) 'To which file?'
        read(*,1000) file
        open(2,file=file)
        write(2,9000) time,dt,nx,dx
        write(2,9000) rho1,rho2,xrho1,xrho2
        write(2,9000) vel1,vel2,xvel1,xvel2
        write(2,9000) dummy
        do 900 i=1,nx
          write(2,*) i,(p(ikomp,i),ikomp=1,nkomp)
900     continue
        close(2)
      end if

      if (ind.eq.14) then
        write(*,1001) eps 
        read(*,*,err=100) eps 
      end if

      if (ind.eq.15) then
        write(*,15500) ampmin,ampmax
        read(*,*,err=100) ampmin,ampmax
      end if

      if (ind.eq.16) call defaul

      if (ind.eq.17) time=.0d0

      if (ind.eq.18) then
        write(*,18000) disamp,lag
        read(*,*) disamp,lag
      end if

      if (ind.eq.21) then
        open(3,file='wave.plot')
        do 21000 ix=1,nx
          write(3,1001) ix*dx,p(1,ix)
21000   continue
        write(3,*) ' p(2,ix):'
        do 21001 ix=1,nx
          write(3,1001) ix*dx,p(2,ix)
21001   continue
        write(3,*) ' rho:'
        do 21002 ix=1,nx
          write(3,1001) ix*dx,rho(ix)
21002   continue
        write(3,*) ' vel:'
        do 21003 ix=1,nx
          write(3,1001) ix*dx,vel(ix)
21003   continue
        write(3,*) ' Nullinie:'
        write(3,1001) .0,.0
        write(3,1001) nx*dx,.0
        close(3)
      end if


      if (ind.ne.0) goto 100

      if (dt.le..1e-15) goto 100

      timeou=time+dtout
      call go(time,timeou)

      goto 100

1000  format(a65)
1001  format(5e12.4)
1002  format(a1)
1500  format(' imeth =',i4,' accu =',i3)
1510  format(' rho =',2f8.4,' xrho =',2f8.4)
1515  format(' vel =',2f8.4,' xvel =',2f8.4)
1520  format(' nx =',i4,' dx =',f8.4,' tmax =',f10.3,' dt =',f8.4)
2000  format(' length =',e10.4,' nx =',i4,' tmax =',e10.4,
     &       ' plam =',e12.4)
3000  format(' rho =',2e12.4,' xrho =',2e12.4)
4000  format(' Method: ',/,
     &  ' -1 = backward        0 = Euler         1 = forward',/,
     &  ' 2  = Lax-Friedrichs  3 = Lax-Wendroff  4 = Leap-Frog',/,
     &  ' 5 = MacCormack       6 = Runge-Kutta   imeth =',i4)
4010  format(' Spatial accuracy (2/4) accu =',i3)
5000  format(' ntfrog =',i8,' lfinim =',i5,' lfavin =',i5,
     &       ' metbot =',i3, ' (1=Lax Friedrich 2=Euler)')
15500 format(' ampmin =',e12.4,' ampmax =',e12.4)
6000  format(' ntplot =',i4)
7000  format(' vel =',2e12.4,' xvel =',2e12.4)
9000  format(8e12.4)
18000 format(' Artificial dissipation in (2,2)-LF: disamp =',e14.6,/,
     &       ' lag =',i4,' (0/1)')
      
      end



      subroutine go(time,t1)
       
      implicit none

      integer*4 ntplot,nx,nxm,nkomp,nt,lgraph
      parameter(nxm=6400,nkomp=2)
     
      real*8 p,pnew,pold,vel,dx,dt,eps,tmax,pi,time,t1
      real*8 ampmin,ampmax,abst,rho,rho1,rho2,delta,xrho1,xrho2
      real*8 peakcu

      real dummy,dtime,tarray(2)

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/discre/dx,dt,eps,tmax,nx
      common/graph/ampmin,ampmax,ntplot,lgraph
      common/const/pi,peakcu

      nt=0

c     dummy=dtime(tarray)

100   continue

      call step(time,nt,abst)
 
      if (lgraph.eq.1.and.mod(nt,ntplot).eq.0) then
        call plot(time)
        call peak(time,peakcu)
        call inquire(*900)
      end if

      if (time.ge.t1.or.abst.lt.eps.or.abst.gt.1e5) then
c       dummy=dtime(tarray)
c        write(*,1000) dummy,tarray(1),tarray(2)
        write(11,*) 'nx =',nx
c        write(12,*) 'nx =',nx
        write(13,*) 'nx =',nx
        write(14,*) 'nx =',nx
        if (lgraph.eq.1) call plot(time)
        call peak(time,peakcu)
        return
      end if

      goto 100

900   continue

      return 

1000  format(' computation time =',3e14.7)
      end

  
      
      subroutine ic
      parameter(nxm=6400,nkomp=2)

c     implicit none

      real*8 p,pnew,pold,vel,dx,dt,eps,tmax,ampmin,ampmax,pi,time
      real*8 dummy,amp,delta,q,rho,rho1,rho2,uex,vex,xi,xrho1,xrho2
      real*8 peakcu

      integer*4 ntplot,nx,nxm,nkomp,ind,ix,nwave,accu,lfavin
      integer*4 ixd,ikomp,lfstar,ntfrog,lgraph,lfinim

      character*65 file

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/discre/dx,dt,eps,tmax,nx
      common/graph/ampmin,ampmax,ntplot,lgraph
      common/const/pi,peakcu
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu

50    write(*,*) ' 1: from file       2: periodic   '
      write(*,*) ' 7: Soliton'
      read(*,*,err=50) ind


      if (ind.eq.1) then
        write(*,*) '  From which file?'
        read(*,1000) file
        open(2,file=file,err=50,status='old')
        read(2,*,err=50) time,dt,nx,dx
        read(2,*,err=50) vel,rho1,rho2,xrho1,xrho2
        read(2,*,err=50) dummy
        do 100 ix=1,nxm
          read(2,*,end=110) ixd,(p(ikomp,ix),ikomp=1,nkomp)
100     continue
110     continue
        nx=ix-1
        call dense
        close(2)
        
        lfstar=0
        return

      end if


      if (ind.eq.2) then

        write(*,*) ' How many wavelengths, what amplitudes?'
        write(*,*) ' Give nwave,amp:'
        read(*,*,err=50) nwave,amp

        waveq=2*pi*nwave/(nx*dx)

        do 200 ix=1,nx

          xi=ix*dx 
          p(1,ix)=amp*cos(waveq*xi)
          p(2,ix)=amp*cos(waveq*xi)
 
200     continue

        lfstar=0
        return

      end if
  

      if (ind.eq.7) then

        write(*,*) '  Amplitude and Width:'
        read(*,*,err=50) amp,delta
 
        do 710 ix=1,nx

          xi=ix*dx 
          p(1,ix)=amp*uex(xi,0)
          p(2,ix)=amp*vex(xi,0)
          pold(1,ix)=amp*uex(xi,-dt)
          pold(2,ix)=amp*vex(xi,-dt)

710     continue
 
        lfstar=0
        return

      end if
      
1000  format(a65)
 
      end 


      subroutine step(time,nt,abst)
      
      implicit none

      integer*4 nx,nxm,nkomp,ikomp,ix,ixm,ixp,imeth,lfstar,iupd

      parameter(nxm=6400,nkomp=2)

      real*8 p,pnew,pold,plam,dt,dx,eps,tmax,time,abst,vel,rho,energy
      real*8 rho1,rho2,delta,ph,xrho1,xrho2,vel1,vel2,xvel1,xvel2,err

      integer*4 nt,ntp,ntfrog,lfinim,accu,lfaver,lfavin,metbot
      integer*4 ixpp,ixmm,ist

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/discre/dx,dt,eps,tmax,nx
      common/method/imeth
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/diag/energy,err

      real*8 ux
      dimension ph(nkomp,nxm),ux(nxm)

      plam=dt/dx

c     downwind Euler

      if (imeth.eq.-1) then

      do 6 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(ikomp,ix)=
     &      p(ikomp,ix)+plam*(p(ikomp,ix)-p(ikomp,ixm))

6     continue

      end if

c     forward Euler

      if (imeth.eq.0) then

      do 11 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(1,ix)=p(1,ix)+
     &             plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix)) 
        pnew(2,ix)=p(2,ix)+
     &             0.5d0*vel(ix)**2*rho(ix)*plam*(p(1,ixp)-p(1,ixm))  
 
11    continue

      end if

c     upwind Euler

      if (imeth.eq.1) then

      do 110 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(ikomp,ix)=
     &      p(ikomp,ix)+plam*(p(ikomp,ixp)-p(ikomp,ix))

110   continue

      end if

      if (imeth.eq.2) then

c     Lax-Friedrichs

      do 210 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(1,ix)=0.5d0*(p(1,ixp)+p(1,ixm))+
     &             plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix)) 
        pnew(2,ix)=0.5d0*(p(2,ixp)+p(2,ixm))+
     &             0.5d0*vel(ix)**2*rho(ix)*plam*(p(1,ixp)-p(1,ixm))  

210   continue

      end if

      if (imeth.eq.3) then

c     Lax-Wendroff

      do 310 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(1,ix)=p(1,ix)+plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix))
     &      +plam**2*(p(1,ixp)-2*p(1,ix)+p(1,ixm))/2
        pnew(2,ix)=p(2,ix)
     &      +plam*vel(ix)**2*rho(ix)*(p(1,ixp)-p(1,ixm))/2
     &      +plam**2*(p(2,ixp)-2*p(2,ix)+p(2,ixm))/2

310   continue

      end if

      if (imeth.eq.4) then

c     Leap Frog

        if (lfstar.eq.0) then

          call lfinit(time)
          lfstar=1
          nt=nt+1
          time=time+dt
          abst=.1d4

          return

        else

          ntp=nt+1
          if (mod(ntp,ntfrog).eq.0) then
            lfaver=1
          else
            lfaver=0
          end if

          iupd=0
          call frog(time,dt,iupd,lfaver)
          
        end if

      end if

      if (imeth.eq.5) then

c     MacCormack

       if (accu.eq.2) then

        if (mod(nt,2).eq.0) then
  
         do 500 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx

          ph(1,ix)=p(1,ix)+plam*(p(2,ixp)-p(2,ix))/rho(ix)
          ph(2,ix)=p(2,ix)+plam*rho(ix)*vel(ix)**2*(p(1,ixp)-p(1,ix))

500      continue

         do 510 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx

          pnew(1,ix)=.5d0*(p(1,ix)+ph(1,ix)
     &                 +plam*(ph(2,ix)-ph(2,ixm))/rho(ix))
          pnew(2,ix)=.5d0*(p(2,ix)+ph(2,ix)
     &                 +plam*rho(ix)*vel(ix)**2*(ph(1,ix)-ph(1,ixm)))

510      continue 

        else

         do 520 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx

          ph(1,ix)=p(1,ix)+plam*(p(2,ix)-p(2,ixm))/rho(ix)
          ph(2,ix)=p(2,ix)+plam*rho(ix)*vel(ix)**2*(p(1,ix)-p(1,ixm))

520      continue

         do 530 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx

          pnew(1,ix)=.5d0*(p(1,ix)+ph(1,ix)
     &                 +plam*(ph(2,ixp)-ph(2,ix))/rho(ix))
          pnew(2,ix)=.5d0*(p(2,ix)+ph(2,ix)
     &                 +plam*rho(ix)*vel(ix)**2*(ph(1,ixp)-ph(1,ix)))

530      continue

        end if

       end if

       if (accu.eq.4) then

        if (mod(nt,2).eq.0) then
  
         do 540 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx
          ixpp=ixp+1
          if (ixpp.gt.nx) ixpp=ixpp-nx
          ixmm=ixm-1
          if (ixmm.lt.1)  ixmm=ixmm+nx

          ph(1,ix)=p(1,ix)+
     &       plam*(-p(2,ixpp)+8*p(2,ixp)-7*p(2,ix))
     &       /(rho(ix)*6)
          ph(2,ix)=p(2,ix)+
     &       plam*rho(ix)*vel(ix)**2*
     &       (-p(1,ixpp)+8*p(1,ixp)-7*p(1,ix))/6

540      continue

         do 550 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx
          ixpp=ixp+1
          if (ixpp.gt.nx) ixpp=ixpp-nx
          ixmm=ixm-1
          if (ixmm.lt.1)  ixmm=ixmm+nx

          pnew(1,ix)=.5d0*(p(1,ix)+ph(1,ix)
     &     +plam*(7*ph(2,ix)-8*ph(2,ixm)+ph(2,ixmm))
     &     /(rho(ix)*6))
          pnew(2,ix)=.5d0*(p(2,ix)+ph(2,ix)
     &     +plam*rho(ix)*vel(ix)**2*
     &     (7*ph(1,ix)-8*ph(1,ixm)+ph(1,ixmm))/6)

550      continue 

        else

         do 560 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx
          ixpp=ixp+1
          if (ixpp.gt.nx) ixpp=ixpp-nx
          ixmm=ixm-1
          if (ixmm.lt.1)  ixmm=ixmm+nx

          ph(1,ix)=p(1,ix)+
     &     plam*(7*p(2,ix)-8*p(2,ixm)+p(2,ixmm))
     &     /(rho(ix)*6)
          ph(2,ix)=p(2,ix)+
     &     plam*rho(ix)*vel(ix)**2*
     &     (7*p(1,ix)-8*p(1,ixm)+p(1,ixmm))/6

560      continue

         do 570 ix=1,nx

          ixp=ix+1
          if (ixp.gt.nx) ixp=ixp-nx
          ixm=ix-1
          if (ixm.lt.1)  ixm=ixm+nx
          ixpp=ixp+1
          if (ixpp.gt.nx) ixpp=ixpp-nx
          ixmm=ixm-1
          if (ixmm.lt.1) ixmm=ixmm+nx

          pnew(1,ix)=.5d0*(p(1,ix)+ph(1,ix)
     &     +plam*(-ph(2,ixpp)+8*ph(2,ixp)-7*ph(2,ix))
     &      /(rho(ix)*6))
          pnew(2,ix)=.5d0*(p(2,ix)+ph(2,ix)
     &     +plam*rho(ix)*vel(ix)**2*
     &      (-ph(1,ixpp)+8*ph(1,ixp)-7*ph(1,ix))/6)

570      continue

        end if

       end if

      end if

      if (imeth.eq.6) then

        call rk4(time,p,ph,.0d0)
        do 610 ikomp=1,nkomp
        do 610 ix=1,nx
          pnew(ikomp,ix)=p(ikomp,ix)+ph(ikomp,ix)/6.d0
610     continue

        call rk4(time+dt/2,p,ph,.5d0)
        do 620 ikomp=1,nkomp
        do 620 ix=1,nx
          pnew(ikomp,ix)=pnew(ikomp,ix)+ph(ikomp,ix)/3.d0
620     continue

        call rk4(time+dt/2,p,ph,.5d0)
        do 630 ikomp=1,nkomp
        do 630 ix=1,nx
          pnew(ikomp,ix)=pnew(ikomp,ix)+ph(ikomp,ix)/3.d0
630     continue

        call rk4(time+dt,p,ph,1.d0)
        do 640 ikomp=1,nkomp
        do 640 ix=1,nx
          pnew(ikomp,ix)=pnew(ikomp,ix)+ph(ikomp,ix)/6.d0
640     continue

      end if

c     Hopskotch

      if (imeth.eq.7) then

      ist=1+mod(nt,2)

      do 71 ix=ist,nx,2
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(1,ix)=p(1,ix)+
     &             plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix)) 
        pnew(2,ix)=p(2,ix)+
     &             0.5d0*vel(ix)**2*rho(ix)*plam*(p(1,ixp)-p(1,ixm))  
 
71    continue

      ist=1+mod(nt+1,2)

      do 72 ix=ist,nx,2
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(1,ix)=p(1,ix)+0.5d0*(
     &             plam*(pnew(2,ixp)-pnew(2,ixm))/(2*rho(ix))+ 
     &             plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix))
     &             ) 
        pnew(2,ix)=p(2,ix)+0.5d0*(
     &        0.5d0*vel(ix)**2*rho(ix)*plam*(pnew(1,ixp)-pnew(1,ixm))+
     &        0.5d0*vel(ix)**2*rho(ix)*plam*(p(1,ixp)-p(1,ixm))
     &             )  
 
72    continue

      end if

      abst=.0d0
      do 900 ikomp=1,nkomp
      do 950 ix=1,nx
     
        abst=abst+(p(ikomp,ix)-pnew(ikomp,ix))**2
        pold(ikomp,ix)=p(ikomp,ix)
        p(ikomp,ix)=pnew(ikomp,ix)

950   continue
900   continue
      abst=sqrt(abst)

      energy=.0d0
      do 960 ix=1,nx
        energy=energy+rho(ix)*vel(ix)**2*p(1,ix)**2+p(2,ix)**2/rho(ix)
960   continue
      energy=energy/nx

      nt=nt+1
      time=time+dt 

      return

      end
    
       
      subroutine rk4(time,p,ph,faktor)

      parameter(nxm=6400,nkomp=2)

c     implicit none

      real*8 time,p,ph,faktor,plam,pl2,vel,delta,temp
      real*8 rho,rho1,rho2,xrho1,xrho2,dx,dt,eps,tmax
      real*8 vel1,vel2,xvel1,xvel2
      
      integer*4 nx,ix,lfstar,lfinim,ntfrog,lfavin,metbot,accu,ikomp
      integer ixp,ixm,ixpp,ixmm

      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu

      dimension p(nkomp,nxm),ph(nkomp,nxm),temp(nkomp,nxm)

      do 50 ikomp=1,nkomp
      do 60 ix=1,nx
        ph(ikomp,ix)=p(ikomp,ix)+faktor*ph(ikomp,ix)
60    continue
50    continue

      pl2=dt/(2*dx)

      if (accu.eq.2) then

      do 110 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        temp(1,ix)=pl2*(ph(2,ixp)-ph(2,ixm))/rho(ix)
        temp(2,ix)=pl2*(ph(1,ixp)-ph(1,ixm))*rho(ix)*vel(ix)**2

110   continue

      end if

      if (accu.eq.4) then

      do 120 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixpp=ix+2
        if (ixpp.gt.nx) ixpp=ixpp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx
        ixmm=ix-2
        if (ixmm.lt.1)  ixmm=ixmm+nx

        temp(1,ix)= 
     &   pl2*(4*(ph(2,ixp)-ph(2,ixm))/3-
     &        (ph(2,ixpp)-ph(2,ixmm))/6)/rho(ix)
        temp(2,ix)= 
     &   pl2*(4*(ph(1,ixp)-ph(1,ixm))/3-
     &        (ph(1,ixpp)-ph(1,ixmm))/6)*rho(ix)*vel(ix)**2

120   continue

      end if

      do 200 ikomp=1,nkomp
      do 210 ix=1,nx
        ph(ikomp,ix)=temp(ikomp,ix)
210   continue
200   continue

      return
      end


      subroutine lfinit(time)
      
      parameter(nxm=6400,nkomp=2)

c     implicit none

      real*8 p,pnew,pold,plam,deltat,dt,eps,tmax,dx,vel,rho,rho1,rho2
      real*8 delta,time,xi,uex,vex,delt,xrho1,xrho2
      real*8 vel1,vel2,xvel1,xvel2

      integer*4 nx,nxm,ikomp,nkomp,ix,ixp,ixm,iupd
      integer*4 lfstar,lfinim,ntfrog,accu,lfavin,metbot

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/discre/dx,dt,eps,tmax,nx
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2

c     Lax-Friedrichs or Euler with dt/lfinim

      if (lfinim.gt.0) then

       deltat=dt/lfinim
       delt=.0d0
       plam=deltat/dx

       do 100 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        if (metbot.eq.1) then
         pnew(1,ix)=(p(1,ixp)+p(1,ixm))/2
     &      +plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix))
         pnew(2,ix)=(p(2,ixp)+p(2,ixm))/2
     &      +plam*rho(ix)*vel(ix)**2*(p(1,ixp)-p(1,ixm))/2
        end if
        if (metbot.eq.2) then
         pnew(1,ix)=p(1,ix)
     &      +plam*(p(2,ixp)-p(2,ixm))/(2*rho(ix))
         pnew(2,ix)=p(2,ix)
     &      +plam*rho(ix)*vel(ix)**2*(p(1,ixp)-p(1,ixm))/2
        end if


          
c          xi=ix*dx 
c          pnew(1,ix)=uex(xi,deltat)
c          pnew(2,ix)=vex(xi,deltat)
c          p(1,ix)=uex(xi,0)
c          p(2,ix)=vex(xi,0)
 
100    continue
c        write(*,*) ' using exact solution in bootstrap'

       do 110 ikomp=1,nkomp
       do 150 ix=1,nx
      
         pold(ikomp,ix)=p(ikomp,ix)
         p(ikomp,ix)=pnew(ikomp,ix)

150    continue
110    continue
       delt=delt+deltat

170    continue

       write(*,1000) deltat,dt,delt

       if (deltat.ge.dt) return

         call frog(time,deltat,0,0)
         do 200 ikomp=1,nkomp
         do 210 ix=1,nx
           p(ikomp,ix)=pnew(ikomp,ix)
210      continue
200      continue
         delt=delt+deltat
         call plot(time)
         deltat=2*deltat

       goto 170

      end if
    
      if (lfinim.le.0) then

       iupd=1
       deltat=dt
       call frog(time,deltat,iupd,0)

       return

      end if
1000  format(' deltat =',e14.6,' dt =',e14.6,' delt =',e14.6)
      end



      subroutine frog(time,deltat,iupd,lfaver)
      
      parameter(nxm=6400,nkomp=2)

c     implicit none

      real*8 p,pnew,pold,plam,dx,dt,eps,tmax,deltat,rho,vel,rho1,rho2
      real*8 delta,time,xrho1,xrho2,vel1,vel2,xvel1,xvel2,disamp

      integer*4 nx,nxm,ikomp,nkomp,ix,ixp,ixm,iupd
      integer*4 accu,lfaver,ixpp,ixmm,lfstar,lfinim,ntfrog,lfavin,metbot

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/discre/dx,dt,eps,tmax,nx
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu
      common/dissip/disamp,lag

      plam=deltat/dx

      if (accu.eq.2) then

      do 110 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx

        pnew(1,ix)=
     &    pold(1,ix)+plam*(p(2,ixp)-p(2,ixm))/rho(ix)
        pnew(2,ix)=
     &    pold(2,ix)+plam*(p(1,ixp)-p(1,ixm))*rho(ix)*vel(ix)**2

c       artificial dissipation

        if (abs(disamp).ge..1e-20) then

          ixpp=ix+2
          if (ixpp.gt.nx) ixpp=ixpp-nx
          ixmm=ix-2
          if (ixmm.lt.1)  ixmm=ixmm+nx
          if (lag.eq.0) then
         pnew(1,ix)=pnew(1,ix)-disamp*dt/(8*dx**4)*
     &    (p(1,ixpp)-4*p(1,ixp)+6*p(1,ix)-4*p(1,ixm)+p(1,ixmm))
         pnew(2,ix)=pnew(2,ix)-disamp*dt/(8*dx**4)*
     &    (p(2,ixpp)-4*p(2,ixp)+6*p(2,ix)-4*p(2,ixm)+p(2,ixmm))
           else
         pnew(1,ix)=pnew(1,ix)-disamp*dt/(8*dx**4)*
     &    (pold(1,ixpp)-4*pold(1,ixp)+6*pold(1,ix)-
     &     4*pold(1,ixm)+pold(1,ixmm))
         pnew(2,ix)=pnew(2,ix)-disamp*dt/(8*dx**4)*
     &    (pold(2,ixpp)-4*pold(2,ixp)+6*pold(2,ix)-
     &     4*pold(2,ixm)+pold(2,ixmm))
           end if
       end if

110   continue

      end if

      if (accu.eq.4) then

      do 120 ix=1,nx
    
        ixp=ix+1
        if (ixp.gt.nx) ixp=ixp-nx
        ixpp=ix+2
        if (ixpp.gt.nx) ixpp=ixpp-nx
        ixm=ix-1
        if (ixm.lt.1)  ixm=ixm+nx
        ixmm=ix-2
        if (ixmm.lt.1)  ixmm=ixmm+nx

        pnew(1,ix)=pold(1,ix)+
     &   plam*(4*(p(2,ixp)-p(2,ixm))/3-
     &        (p(2,ixpp)-p(2,ixmm))/6)/rho(ix)
        pnew(2,ix)=pold(2,ix)+
     &   plam*(4*(p(1,ixp)-p(1,ixm))/3-
     &        (p(1,ixpp)-p(1,ixmm))/6)*rho(ix)*vel(ix)**2

120   continue

      end if

      if (lfaver.eq.1) then

       do 400 ikomp=1,nkomp
       do 410 ix=1,nx
     
         pnew(ikomp,ix)=.5d0*(pnew(ikomp,ix)+p(ikomp,ix))
         p(ikomp,ix)=   .5d0*(p(ikomp,ix)+pold(ikomp,ix))

410    continue
400    continue

       time=time-deltat/2
       write(*,1000) time

      end if

      if (iupd.eq.1) then

        do 200 ikomp=1,nkomp
        do 250 ix=1,nx
     
          pold(ikomp,ix)=p(ikomp,ix)
          p(ikomp,ix)=pnew(ikomp,ix)

250     continue
200     continue

      end if

      return
1000  format(' Averaging at t =',e14.6)
      end

      

      subroutine plot(time)
      parameter(nkomp=2,nxm=6400)
      implicit real*8 (a-h,p-z)
      implicit real*4 (o)
      character*35 text
      integer BLACK, GREEN, RED, BLUE, WHITE, YELLOW
      parameter (BLACK=0,RED=1,GREEN=2,YELLOW=3,BLUE=4,
     &           MAGENTA=5,CYAN=6,WHITE=7)
      common/plotf/x(nxm),empty(nxm)
      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx
      common/graph/ampmin,ampmax,ntplot,lgraph
      common/const/pi,peakcu
      common/diag/energy,err
      common/grenze/xmin,xmax,ymin,ymax,delx,dely,lexx,lexy
      common/delta/xminh,xmaxh,yminh,ymaxh,zminh,zmaxh

      dimension p1(nxm),p2(nxm),ux(nxm),axx(2),axy(2)

      do 100 i=1,nx
        p1(i)=p(1,i)
        p2(i)=p(2,i)
100   continue

      nxmal=nx
      lexx=0
      lexy=1

      xmin=x(1)
      xmax=x(nx)
      delx=xmax-xmin
      ymin=ampmin
      ymax=ampmax
      dely=ymax-ymin

      call color(BLACK)
      call clear
      call color(WHITE)

      axx(1)=.0
      axx(2)=x(nx)
      axy(1)=.0
      axy(2)=.0
      call curve2(axx,axy,2)

      axy(1)=1.
      axy(2)=1.
      call curve2(axx,axy,2)
      axy(1)=-1.
      axy(2)=-1.
      call curve2(axx,axy,2)


      call color(RED)
      call curve2(x,p1,nxmal)

c      call font('futura.m')
c      call textsize(1.5,1.5)

      call color(RED)
      call move2(-.8,.9)
      write(text,1000) yminh
      call drawstr(text)
      call drawstr(' <  p  <')
      write(text,1000) ymaxh
      call drawstr(text)

      call color(YELLOW)
      call curve2(x,p2,nxmal)

      call color(YELLOW)
      call move2(-.8,.84)
      write(text,1000) yminh
      call drawstr(text)
      call drawstr(' <  p  <')
      write(text,1000) ymaxh
      call drawstr(text)

      call color(GREEN)
      call curve2(x,rho,nxmal)

      call color(BLUE)
      call curve2(x,vel,nxmal)

      call color(GREEN)
      call move2(-.8,.78)
      write(text,1000) yminh
      call drawstr(text)
      call drawstr(' <  p  <')
      write(text,1000) ymaxh
      call drawstr(text)

      errq=.0d0
      err0=.0d0
      do 200 ix=1,nx
        xi=ix*dx 
        ux(ix)=uex(xi,time)   
        errq=errq+(p(1,ix)-ux(ix))**2
        err0=err0+uex(xi,time)**2
200   continue
      err=sqrt(errq/err0)

      write(11,1001) time,err

      call color(MAGENTA)
      call curve2(x,ux,nxmal)

      call color(WHITE)
      call move2(.2,.9)
      write(text,1000) time
      call drawstr(text)
      call move2(.2,.84)
      write(text,1000) energy
      call drawstr(text)
      call move2(.2,.78)
      write(text,1000) err
      call drawstr(text)

      call swapbuffers

      return
1000  format(e12.4)
1001  format(2e14.6)
      end




      subroutine initialize
      parameter(nxm=6400)
      implicit real*8 (a-h,o-z)
      real*4 xl,xr
      character*50 device
      integer BLACK, GREEN,RED,YELLOW,BLUE
      parameter (BLACK = 0, RED = 1, GREEN = 2, YELLOW = 3,BLUE = 4)
      integer backbuffer
      common/plotf/x(nxm),empty(nxm)
 
      write(*,*) ' Monitor or Projector? (0/1)'
      read(*,*) lmon

      if (lmon.eq.1) then

       lang=1200
       lhoch=800
       aspekt=float(lang)/float(lhoch)
       call prefsize(lang,lhoch)
       call prefposition(10,10)

      else 

       lang=800
       lhoch=400
       aspekt=float(lang)/float(lhoch)
       call prefsize(lang,lhoch)
       call prefposition(50,50)

      end if

c     Initialisiere Graphik

c      print*,'Enter output device:'
c      read(*,'(a)')device

      device='X11'
      call vinit(device)

      call color(BLACK)
      call clear

        xl=-1.1
        xr=xl+2.2/aspekt
        call ortho2(xl,xr,-1.1,1.1)
        call clipping(0)

      if (lmon.eq.0) then

	call move2(-0.9, 0.1)
	call color(0)
	call draw2(-.8, 0.3)
	call color(1)
	call draw2(-.7, 0.4)
	call color(2)
	call draw2(-.6, 0.5)
	call color(3)
	call draw2(-.5, 0.6)
	call color(4)
	call draw2(.0, 0.8)
	call color(5)
	call draw2(.3, 0.1)
	call color(6)
	call draw2(.8, 1.3)
	call color(7)
	call draw2(1.6, 0.1)
	call color(GREEN)

      end if

	if (backbuffer().lt.0) then
		call vexit
		write(*,*)'Device can''t support doublebuffering'
		stop
	endif

      do 100 i=1,nxm
        x(i)=float(i)
100   continue

      return
      end


        subroutine curve2(x,y,nx)
        parameter(nxm=6400)
        implicit real*8 (a-h,o-z)
        dimension x(nxm),x1(nxm),y(nxm),y1(nxm)
        real*4 xi,yi
        common/grenze/xmi,xma,ymi,yma,delx,dely,lexx,lexy
        common/delta/xminh,xmaxh,yminh,ymaxh,zminh,zmaxh

c
c       Extrema bestimmen
c
        xminh=1.e38
        xmaxh=-1.e38
        yminh=1.e38
        ymaxh=-1.e38
        do 100 i=1,nx
          xminh=min(xminh,x(i))
          xmaxh=max(xmaxh,x(i))
          yminh=min(yminh,y(i))
          ymaxh=max(ymaxh,y(i))
100     continue

        if (lexx.eq.0) then
          xmin=xminh
          xmax=xmaxh
        else
          xmin=xmi
          xmax=xma
        end if

        if (lexy.eq.0) then
          ymin=yminh
          ymax=ymaxh
        else
          ymin=ymi
          ymax=yma
        end if

        xmit=(xmax+xmin)/2.d0
        ymit=(ymax+ymin)/2.d0

        do 200 i=1,nx
          x1(i)=x(i)-xmit
          y1(i)=y(i)-ymit
200     continue

        delx=xmax-xmin
        dely=ymax-ymin

        if (delx.lt..1e-10) delx=.1*xmit
        if (dely.lt..1e-10) dely=.1*ymit

          xi=x1(1)*2./delx
          yi=y1(1)*2./dely

        call move2(xi,yi)
        do 500 i=2,nx
          xi=x1(i)*2./delx
          yi=y1(i)*2./dely
          call draw2(xi,yi)
500     continue

        return
        end



      subroutine defaul
      parameter(nxm=6400,nkomp=2)
      
c     implicit none

      real*8 p,pnew,pold,vel,dx,dt,eps,tmax,ampmin,ampmax,pi
      real*8 rho,rho1,rho2,delta,xrho1,xrho2,peakcu
      real*8 vel1,vel2,xvel1,xvel2

      integer*4 ntplot,nx,nxm,nkomp,imeth,ntfrog,lfstar,lgraph
      integer*4 lfinim,accu,lfaver,lfavin

      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)
      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx
      common/graph/ampmin,ampmax,ntplot,lgraph
      common/const/pi,peakcu
      common/method/imeth
      common/leapfr/lfstar,lfinim,ntfrog,lfavin,metbot,accu

      write(*,1000) nkomp,nxm
         
      open(1,file='wave.def',err=100,status='old')
      read(1,*) dx,nx,tmax,dt
      read(1,*) rho1,rho2,xrho1,xrho2
      read(1,*) vel1,vel2,xvel1,xvel2
      read(1,*) ntplot,ampmin,ampmax
      read(1,*) eps,peakcu
      read(1,*) imeth,ntfrog,lfinim,lfaver,accu,metbot
      close(1)

      if (nx.gt.nxm) then
        write(*,*) ' nx = nxm gesetzt.'
        nx=nxm
      end if 

      call dense 
 
      return

100   continue 

      write(*,*) ' Problem with opening wave.def.'

      return
1000  format('  Compiled with nkomp =',i3,' and nxm =',i4)
      end




      subroutine constan
      
c     implicit none

      real*8 pi,peakcu

      common/const/pi,peakcu
 
      pi=2.d0*acos(.0d0)
 
      return
      end
         

      function uex(x,t)
      parameter(nxm=6400,nkomp=2)

c     implicit none

      real*8 x,t,uex,vel,delta,rho1,rho2,rho,arg,dx,dt,eps,tmax
      real*8 posit,xrho1,xrho2,vel1,vel2,xvel1,xvel2
      integer*4 nx

      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx
     
c     bessere Methode: variiere Argument der Gauss Funktion und
c     berechne das dazugehoerige x anschliessend

c     posit=mod(x-nx*dx/2+vel1*t,nx*dx)
      posit=x-nx*dx/2+vel1*t
100   continue
      if (posit.lt.-nx*dx/2) then
        posit=posit+dx*nx
        goto 100
      end if 
      if (posit.gt.nx*dx/2) then
        posit=posit-dx*nx
        goto 100
      end if

      arg=posit**2/(2*delta**2)
      if (arg.gt.50) then
        uex=0
      else
        uex=exp(-arg)
      end if

      return

      end 


      function vex(x,t)
      parameter(nxm=6400,nkomp=2)

c     implicit none

      real*8 x,t,vex,vel,delta,rho1,rho2,rho,arg,dx,dt,eps,tmax
      real*8 posit,xrho1,xrho2,vel1,vel2,xvel1,xvel2

      integer*4 nx,ix

      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx

c      if (rho1.ne.rho2) write(*,*) ' rho1 and rho2 are different'

      ix=int(x/dx)

c     posit=mod(x-nx*dx/2+vel1*t,nx*dx)
      posit=x-nx*dx/2+vel1*t
100   continue
      if (posit.lt.-nx*dx/2) then
        posit=posit+dx*nx
        goto 100
      end if 
      if (posit.gt.nx*dx/2) then
        posit=posit-dx*nx
        goto 100
      end if

      arg=(posit)**2/(2*delta**2)
      if (arg.gt.50) then
        vex=0
      else
        vex=vel1*rho(ix)*exp(-arg)
      end if

      return

      end 



      subroutine dense
      parameter(nxm=6400,nkomp=2)
c     implicit none

      real*8 rho1,rho2,rho,dx,dt,eps,tmax,vel,delta,xi,xrho1,xrho2
      real*8 vel1,vel2,xvel1,xvel2

      integer*4 nx,ix

      common/parrho/delta,rho(nxm),rho1,rho2,xrho1,xrho2
      common/parvel/vel(nxm),vel1,vel2,xvel1,xvel2
      common/discre/dx,dt,eps,tmax,nx

      do 300 ix=1,nx
        xi=ix*dx
        if (xi.lt.xrho1.or.xi.gt.xrho2) then
          rho(ix)=rho1
        else
          rho(ix)=rho2
        end if
300   continue
 
      do 400 ix=1,nx
        xi=ix*dx
        if (xi.lt.xvel1.or.xi.gt.xvel2) then
          vel(ix)=vel1
        else
          vel(ix)=vel1+vel2*sin((xi-xvel1)*6.28/(xvel2-xvel1))
        end if
400   continue

      return
      end


      subroutine peak(time,peakcu)

      parameter(nxm=6400,nkomp=2)
c     implicit none

      real*8 dx,dt,eps,tmax,p,pnew,pold,time,y1,y2,y3,cex,bex,aex
      real*8 xpuls,ypuls,peakcu,p1h,p1v

      integer*4 nx,ix

      common/discre/dx,dt,eps,tmax,nx
      common/sol/p(nkomp,nxm),pnew(nkomp,nxm),pold(nkomp,nxm)


      do 400 ix=2,nx-1
        if (p(1,ix).gt.p(1,ix-1).and.p(1,ix).gt.p(1,ix+1).or.
     &      p(1,ix).lt.p(1,ix-1).and.p(1,ix).lt.p(1,ix+1)) then
           y1=p(1,ix-1)
           y2=p(1,ix)
           y3=p(1,ix+1)
           cex = -y2+y1/2+y3/2
           bex = -ix*y1-y1/2+2*y2*ix-ix*y3+y3/2
           aex = -y2*ix**2+y2+ix**2*y1/2+ix*y1/2+ix**2*y3/2-ix*y3/2
           xpuls=-dx*bex/(2.d0*cex)
           ypuls=aex-bex**2/(4.d0*cex)
c           if (abs(ypuls).gt.peakcu) write(12,1000) time,xpuls,ypuls
        end if
400   continue

      p1h=p(1,nx/2)
      p1v=p(1,nx/4)
      if (abs(p1h).lt..1e-30) p1h=.0d0
      if (abs(p1v).lt..1e-30) p1v=.0d0
      write(13,1000) time,p1h 
      write(14,1000) time,p1v  

      return

1000  format(5e17.6)
      end





      subroutine inquire(*)
      implicit real*8 (a-h,p-z)
      implicit real*4 (o)
      external checkkey
      integer checkkey
      common/steu/idruck,itanf,itend,itjump,ixanf,ixend 

      iaus=checkkey()

c      write(21,*) iaus

      if (iaus.eq.0) return  

      write(*,*) iaus

      return 1
      
      return
      end
