program source1d
      integer j,n,nk,k,is
      parameter(n=100)
      double precision t,t1,dt,a,cfl,s,adtx,pi,delta,xj,h
      double precision uold(0:n),unew(0:n)
      common/valupi/pi
      parameter(cfl=0.2D0,a=1.D0)
      data t/0.D0/
      print *,'enter nk:'
      read *,nk
      pi = 4.D0*datan(1.D0)
      h = 1.D0/dfloat(n)
      dt = cfl*+h/a
      adtx = a*dt/h
c      nk = int(t1/dt)
      call initia(uold,n)
  
      uold(0) = 1.D0
      uold(n) = 0.D0

      do 10 k=1,nk
         t = t+dt
         call identu(uold,n,s)
         is = int(s/h)
         do 20 j=is-9,is+10
            xj = dfloat(j)*h
            unew(j) = uold(j) - adtx*(uold(j)-uold(j-1))
     +                + delta(xj-s,h)
 20      continue
         do 30 j=is-9,is+10
            uold(j) = unew(j)
 30      continue
 10   continue
      do 40 j=0,n
         xj = dfloat(j)*h
         write(11,101) xj,uold(j)
 40   continue
 101  format(2(f14.8,2x))
      end

      subroutine initia(u,n)
      integer j,n
      double precision s0,h,xs,pi
      double precision u(0:n)
      common/valupi/pi
      h = 1.D0/dfloat(n)
      s0 = 0.50001D0
      do 10 j=0,n
         xs = dfloat(j)*h - s0
         if(xs .le. -4.D0*h) then
            u(j) = 1.D0
         elseif(xs .ge. 4.D0*h) then
            u(j) = 0.D0
         else
            u(j) = 0.5D0*(1.D0-dsin(pi*xs/8.D0/h))
         endif
         write(10,101) dfloat(j)*h,u(j)
 10   continue
 101  format(2(f14.8,2x))
      end

      
      subroutine identu(u,n,s)
      integer j,n
      double precision s,h
      double precision u(0:n)
      h = 1.D0/dfloat(n)
      do 10 j=1,n
         if(u(j-1)*u(j) .gt. 0.D0) goto 10
         s = dfloat(j)*h - dabs(u(j))/(dabs(u(j-1))+dabs(u(j)))
         return
 10   continue
      end

      function delta(x,h)
      double precision x,h,pi,delta
      common/valupi/pi
      if ( dabs(x) .ge. 2.D0*h ) then
         delta = 0.D0
      else
         delta = 0.25D0*(1.D0-dcos(0.5D0*pi*x/h))/h
      endif
      end