=================== beginning of wave.f ===============================
c+++++++++++  wave.f (ISSSDEMO.FOR) by K. Watanabe ++++++++++++++
C//simldemo job ,msgclass=d,
C//          region=5m,time=(2,0)
C//     exec fort77,option='nos'
C//sysin dd *
c
c
      parameter(nx=49)
      parameter(nxpt=1000)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com02/
     &          nloop,ntend,nwrpf,nlplw,nlprk,ndelt
      common /com03/
     &          wvnmbr
      common /com04/
     &          mdnmbr
      common /com05/
     &          alpha
      common /com11/
     &          bxa(nx),  bxb(nx)
      common /com12/
     &          bx4(nx),  dbx4(nx),  qbx4(nx)
      common /com21/
     &          bx0(nx,6,3),energy(nxpt,2),timeen(nxpt,2)
      common /com22/
     &          timepr(6)
      common /com23/
     &          energ0
c
      call param(icon)
      if(icon.ne.0) stop
c
c     start of 2-step lax-wendroff
c
      icon = 0
      istep1 = 0
      istep2 = 0
      nloop = 0
c
      call reclw(istep1,istep2,icon)
c
 1000 continue
c
      nloop = nloop + 1
      time = dt*float(nloop)
c
      call step1
c
      call step2
c
      call reclw(istep1,istep2,icon)
c
      if(icon.eq.0.and.nloop.lt.ntend) go to 1000
c
c     start of runge-kutta-gill
c
      icon  = 0
      istep1 = 0
      istep2 = 0
      nloop = 0
      time = 0.
c
      call recrk(istep1,istep2,icon)
c
 2000 continue
c
      nloop = nloop + 1
      time = dt*float(nloop)
c
      call caldlt
      call rkg(1)
c
      call caldlt
      call rkg(2)
c
      call caldlt
      call rkg(3)
c
      call caldlt
      call rkg(4)
c
      call recrk(istep1,istep2,icon)
c
      if(icon.eq.0.and.nloop.lt.ntend) go to 2000
c
      call drwgrp(istep1)
c
      stop
      end
c
c***********************************************************************
      subroutine param(icon)
c***********************************************************************
c
      parameter(nx=49)
      parameter(nxpt=1000)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com02/
     &          nloop,ntend,nwrpf,nlplw,nlprk,ndelt
      common /com03/
     &          wvnmbr
      common /com04/
     &          mdnmbr
      common /com05/
     &          alpha
      common /com11/
     &          bxa(nx),  bxb(nx)
      common /com12/
     &          bx4(nx),  dbx4(nx),  qbx4(nx)
      common /com21/
     &          bx0(nx,6,3),energy(nxpt,2),timeen(nxpt,2)
      common /com23/
     &          energ0
c
      namelist /data01/  mdnmbr , alpha
c
      read(5,data01)
c
      if(mdnmbr.ne.3.and.mdnmbr.ne.4.and.mdnmbr.ne.6
     &              .and.mdnmbr.ne.8.and.mdnmbr.ne.12.and.mdnmbr.ne.16
     &              .and.mdnmbr.ne.24.and.mdnmbr.ne.48) then
      write(6,9999)
 9999 format(1h ,'mode number is wrong for this test program'
     &      /1h ,'   please choose it from the following list'
     &      /1h ,'      mdnmbr = 3, 4, 6, 8, 12, 16, 24, or 48')
      icon = 9999
      return
      endif
c
      dx = 1.
c
      nloop = 0
      time  = 0.
c
      pai   = 3.141593
c
      dt    = dx * alpha
c
      ntend = ifix(10. * float(mdnmbr) * dx / dt + 0.5)
      nwrpf = ifix( float(ntend-1) / 4. + 0.5 )
      ndelt = 1 + ntend / nxpt
c
      dtdx   = dt / dx
      dtdx2  = dt / (2. * dx)
      dtdx12 = dt / (12. * dx)
c
      wvnmbr = 2. * pai / float(mdnmbr)
      do 10  i = 1 , nx - 1
      tht = wvnmbr * float(i - 1)
      bxa(i) = sin(tht)
      bxb(i) = 0.
      bx4(i) = bxa(i)
   10 continue
c
      bxa(nx) = bxa(1)
      bx4(nx) = bx4(1)
c
      energ0 = 0.
      do 20  i = 1 , nx-1
      energ0 = energ0 + bx4(i)**2
   20 continue
c
      do 30  i = 1 , 6
      do 30  j = 1 , 3
      bx0(1,i,j) = 1.e+10
   30 continue
c
      write(6,100) alpha , mdnmbr
  100 format(1h0, 'alpha = ',1pe10.3,'    mode number = ',i2)
c
      write(6,200) dx , dt
  200 format(1h0,'dx = ',1pe10.3,2x,'dt = ',1pe10.3)
c
      return
      end
c
c***********************************************************************
      subroutine step1
c***********************************************************************
c
      parameter(nx=49,nxm1=nx-1)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com11/
     &          bxa(nx),  bxb(nx)
c
      do 100  i = 1 , nxm1
c
      ip1 = i + 1
      im1 = i
c
      bxb(i) =   ( bxa(ip1) + bxa(im1) ) / 2
     &         + ( bxa(ip1) - bxa(im1) ) * dtdx2
c
  100 continue
c
      bxb(nx) = bxb(1)
c
      return
      end
c
c***********************************************************************
      subroutine step2
c***********************************************************************
c
      parameter(nx=49,nxm1=nx-1)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com11/
     &          bxa(nx),  bxb(nx)
c
      do 100  i = 2 , nx
c
      ip1 = i
      im1 = i - 1
c
      bxa(i) =     bxa(i)
     &         + ( bxb(ip1) - bxb(im1) ) * dtdx
c
  100 continue
c
      bxa(1) = bxa(nx)
c
      return
      end
c
c***********************************************************************
      subroutine reclw(istep1,istep2,icon)
c***********************************************************************
c
      parameter(nx=49)
      parameter(nxpt=1000)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com02/
     &          nloop,ntend,nwrpf,nlplw,nlprk,ndelt
      common /com03/
     &          wvnmbr
      common /com04/
     &          mdnmbr
      common /com11/
     &          bxa(nx),  bxb(nx)
      common /com21/
     &          bx0(nx,6,3),energy(nxpt,2),timeen(nxpt,2)
      common /com23/
     &          energ0
c
      if(mod(nloop,ndelt).eq.0) then
c
      istep2 = istep2 + 1
      timeen(istep2,1) = time / float(mdnmbr) / dx
      ampmax = -1.e+10
      energy(istep2,1) = 0.
      do 10  i = 1 , nx-1
      ampmax = amax1(bxa(i),ampmax)
      energy(istep2,1) = energy(istep2,1) + bxa(i)**2
   10 continue
      energy(istep2,1) = energy(istep2,1) / energ0
c
      endif
c
      if(ampmax.gt.2.) icon = 9999
c
      if(mod(nloop,nwrpf).eq.0.and.nloop.ne.ntend.and.icon.eq.0) then
c
      istep1 = istep1 + 1
      do 20  i = 1 , nx
      bx0(i,istep1,3) = bxa(i)
   20 continue
c
      else if(nloop.eq.ntend.or.icon.ne.0) then
c
      nlplw = istep2
      istep1 = istep1 + 1
      do 30  i = 1 , nx
      bx0(i,istep1,3) = bxa(i)
   30 continue
c
      timex = time / float(mdnmbr) / dx
      write(6,9000) nloop , timex , icon
 9000 format(1h0,'end of l-w simulation at loop = ',i7,
     &           '  (time = ',1pe10.3,' )    condition = ',i7)
c
      endif
c
      return
      end
c
c***********************************************************************
      subroutine caldlt
c***********************************************************************
c
      parameter(nx=49,nxm1=nx-1)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com12/
     &          bx4(nx),  dbx4(nx),  qbx4(nx)
c
      do 100  m = 1 , nxm1
      im2 = m - 2
      im1 = m - 1
      ip1 = m + 1
      ip2 = m + 2
c
      if(m.eq.1) then
      im1 = nx - 1
      im2 = nx - 2
      endif
      if(m.eq.2) then
      im2 = nx - 1
      endif
      if(m.eq.nxm1) then
      ip2 = 2
      endif
c
      dbx4(m) =
     & (  -      ( bx4(ip2) - bx4(im2) )
     &    + 8. * ( bx4(ip1) - bx4(im1) ) ) * dtdx12
c
  100 continue
c
      return
      end
c
c***********************************************************************
      subroutine rkg(id)
c***********************************************************************
c
      parameter(nx=49,nxm1=nx-1)
c
      common /com12/
     &          bx4(nx),  dbx4(nx),  qbx4(nx)
c
      if      (id.eq.1) then
c
       c1 =  0.5
       c2 = -1.
       cq = -2.
       c0 =  1.
c
      else if (id.eq.2) then
c
       c1 =  1. - sqrt(0.5)
       c2 = -c1
       cq =  1. - 3.*c1
       c0 =  2.*c1
c
      else if (id.eq.3) then
c
       c1 =  1. + sqrt(0.5)
       c2 = -c1
       cq =  1. - 3.*c1
       c0 =  2.*c1
c
      else if (id.eq.4) then
c
       c1 =  1./6.
       c2 = -1./3.
       cq =  0.
       c0 =  0.
c
      end  if
c
      do 10  m = 1 , nxm1
      bx4(m) = bx4(m) + c1*dbx4(m) + c2*qbx4(m)
   10 continue
c
      do 20  m = 1 , nx - 1
      qbx4(m) = cq*qbx4(m) + c0*dbx4(m)
   20 continue
c
      bx4(nx) = bx4(1)
c
      return
      end
c
c***********************************************************************
      subroutine recrk(istep1,istep2,icon)
c***********************************************************************
c
      parameter(nx=49)
      parameter(nxpt=1000)
c
      common /com01/
     &          dt, dx, dtdx, dtdx2, dtdx12, time
      common /com02/
     &          nloop,ntend,nwrpf,nlplw,nlprk,ndelt
      common /com03/
     &          wvnmbr
      common /com04/
     &          mdnmbr
      common /com12/
     &          bx4(nx),  dbx4(nx),  qbx4(nx)
      common /com21/
     &          bx0(nx,6,3),energy(nxpt,2),timeen(nxpt,2)
      common /com22/
     &          timepr(6)
      common /com23/
     &          energ0
c
      if(mod(nloop,ndelt).eq.0) then
c
      istep2 = istep2 + 1
      timeen(istep2,2) = time / float(mdnmbr) / dx
      ampmx4 = -1.e+10
      energy(istep2,2) = 0.
      do 10  i = 1 , nx-1
      ampmx4 = amax1(bx4(i),ampmx4)
      energy(istep2,2) = energy(istep2,2) + bx4(i)**2
   10 continue
      energy(istep2,2) = energy(istep2,2) / energ0
c
      endif
c
      if(ampmx4.gt.2.) icon = 9999
c
      if(mod(nloop,nwrpf).eq.0.and.nloop.ne.ntend.and.icon.eq.0) then
c
      istep1 = istep1 + 1
      do 20  i = 1 , nx
      timepr(istep1) = time / float(mdnmbr) / dx
      bx0(i,istep1,2) = bx4(i)
      tht = wvnmbr * (float(i - 1) - time)
      bx0(i,istep1,1) = sin(tht)
   20 continue
c
      else if(nloop.eq.ntend.or.icon.ne.0) then
c
      nlprk = istep2
      istep1 = istep1 + 1
      do 30  i = 1 , nx
      timepr(istep1) = time / float(mdnmbr) / dx
      bx0(i,istep1,2) = bx4(i)
      tht = wvnmbr * (float(i - 1) - time)
      bx0(i,istep1,1) = sin(tht)
   30 continue
c
      timex = time / float(mdnmbr) / dx
      write(6,9000) nloop , timex , icon
 9000 format(1h0,'end of r-k simulation at loop = ',i7,
     &           '  (time = ',1pe10.3,' )    condition = ',i7)
c
      endif
c
      return
      end
c
c***********************************************************************
      subroutine drwgrp(istep1)
c***********************************************************************
c
      parameter(nx=49)
      parameter(nxpt=1000)
c
      common /com02/
     &          nloop,ntend,nwrpf,nlplw,nlprk,ndelt
      common /com04/
     &          mdnmbr
      common /com05/
     &          alpha
      common /com21/
     &          bx0(nx,6,3),energy(nxpt,2),timeen(nxpt,2)
      common /com22/
     &          timepr(6)
c
      kishu = 1
c      call gopen(kishu,'isss')
      call plots
c
      amdnm = float(mdnmbr) + 0.1
      anlpr = float(nlprk) + 0.1
      call symbol(4.,16.,0.5,'wave propagation test',0.,21)
      call symbol(6.,13.,0.3,'mode number = ',0.,14)
      call number(10.2,13.,0.3,amdnm,0.,-1)
      call symbol(11.1,13.,0.3,'( meshes for one wave length )',0.,30)
      call symbol(6.,12.,0.3,'v * dt / dx = ',0.,14)
      call number(10.2,12.,0.3,alpha,0.,2)
      call symbol(6.,11.,0.3,'nlprk = ',0.,8)
      call number(10.2,11.,0.3,anlpr,0.,-1)
      call chart
c
      xs  = 2.5
      ys0 = 14.5
      xl = 16.
      yl = 4.
      xsl = xs + xl
      bmin = -1.
      bmax =  1.
c
      dx = xl / float(nx - 1)
      dy = yl / (bmax - bmin)
c
      do 1000  j = 1 , istep1
c
      do 100  k = 1 , 3
      ys = ys0 - (yl + 2.)*float(k-1)
      ysl = ys + yl
      yslh = ys + yl/2.
c
      if(bx0(1,j,k).lt.10.) then
c
      call newpen(2)
      do 10  i = 1 , nx
      xx = xs + dx*float(i-1)
      yy = ys + dy*(bx0(i,j,k)-bmin)
c
      if(i.eq.1) then
      call plot(xx,yy,3)
      else
      call plot(xx,yy,2)
      endif
c
   10 continue
c
      call plot(xs,yslh,3)
c     call dashp(xsl,yslh,0.1)
      call dplot(0.,0.,0)
      call dplot(xs,yslh,3)
      call dplot(xsl,yslh,2)

c
      else
c
      call newpen(1)
      call symbol(xs+3.5,yslh-0.5,1.0,'over flow',0.,9)
c
c     call plot(xs,yslh,3)
c     call dashp(xs+2.,yslh,0.1)
c     call plot(xs+14.,yslh,3)
c     call dashp(xsl,yslh,0.1)
      call dplot(0.,0.,0)
      call dplot(xs,yslh,3)
      call dplot(xs+2.,yslh,2)
      call dplot(xs+14.,yslh,3)
      call dplot(xsl,yslh,2)

c
      endif
c
      call newpen(1)
c
      call plot(xs,ys,3)
      call plot(xsl,ys,2)
      call plot(xsl,ysl,2)
      call plot(xs,ysl,2)
      call plot(xs,ys,2)
c
      do 20  i = 1 , nx , 8
      xx = xs + dx*float(i-1)
      call plot(xx,ys,3)
      call plot(xx,ys-0.2,2)
   20 continue
c
      call plot(xs,ys,3)
      call plot(xs-0.2,ys,2)
      call plot(xs,yslh,3)
      call plot(xs-0.2,yslh,2)
      call plot(xs,ysl,3)
      call plot(xs-0.2,ysl,2)
c
      if(k.eq.1) then
      call symbol(xs+0.5,ysl+0.2,0.3,'analytic result',0.,15)
c
      call symbol(xsl+1.,ysl-1.,0.3,'time = ',0.,7)
      call number(xsl+3.1,ysl-1.,0.3,timepr(j),0.,2)
      else if(k.eq.2) then
      call symbol(xs+0.5,ysl+0.2,0.3,'4th order r-k-g result',0.,22)
      else if(k.eq.3) then
      call symbol(xs+0.5,ysl+0.2,0.3,'2step lax-wendroff result',0.,25)
      endif
c
  100 continue
c
      call chart
c
 1000 continue
c
      xs = 4.
      ys = 4.
      xl = 10.
      yl = 10.
      xsl = xs + xl
      ysl = ys + yl
c
      dx = xl / 10.
      dy = yl
c
      do 2000  k = 1 , 2
      if(k.eq.1) then
      call newpen(1)
      else if(k.eq.2) then
      call newpen(2)
      endif
      if(k.eq.1) then
      nnnn = nlplw
      else
      nnnn = nlprk
      endif
      do 2100  i = 1 , nnnn
      xx = xs + dx*timeen(i,k)
      yy = ys + dy*energy(i,k)
      if(i.eq.1) then
      call plot(xx,yy,3)
      else
      call plot(xx,yy,2)
      endif
 2100 continue
 2000 continue
c
      call newpen(1)
c
      call plot(xs,ys,3)
      call plot(xsl,ys,2)
      call plot(xsl,ysl,2)
      call plot(xs,ysl,2)
      call plot(xs,ys,2)
c
      dx = xl / 10.
      do 2200  i = 1 , 11
      xx = xs + dx*float(i-1)
      call plot(xx,ys,3)
      call plot(xx,ys-0.2,2)
      val = float(i-1) + 0.1
      if(i.le.10) then
      call number(xx-0.1,ys-0.5,0.2,val,0.,-1)
      else
      call number(xx-0.2,ys-0.5,0.2,val,0.,-1)
      endif
 2200 continue
c
      dy = yl / 10.
      do 2300  i = 1 , 11
      yy = ys + dy*float(i-1)
      call plot(xs,yy,3)
      call plot(xs-0.2,yy,2)
      val = 0.1*float(i-1) + 0.01
      call number(xs-0.9,yy,0.2,val,0.,1)
 2300 continue
c
      call symbol(xs+xl/2.-0.6,ys-1.0,0.3,'time',0.,4)
      call symbol(xs-1.1,ys+yl/2.-0.9,0.3,'energy',90.,6)
c
c     call gclose
      call plot(0.,0.,999)
c
      return
      end
c//  exec lnknlp
c//lked.sysprint dd dummy
c//  exec gonlp
c//* mdnmbr : mode number which should be
c//*                3, 4, 6, 8, 12, 16, 24, or, 48.
c//* alpha  = vp*dt/dx, where vp, dt and dx are phase velocity,
c//*                    time step and grid size, respectively.
c//go.sysin  dd  *
c &data01  mdnmbr=16,      alpha=0.80,        &end
c//
c
C  dplot routine is for plotting a dashed line. (Y.O.)
c
=================== end of wave.f ==================================