parameter (N=64, NN=2*N+1)
       implicit real*8(a-h,o-z)
       real x(0:NN), ubar(0:NN,0:3),uapp(0:NN)
       real uapp1(0:NN), xL(0:NN,0:3), ubars(0:NN,0:3), si(0:NN)
       complex abar(0:NN), al(0:NN)
       complex w, A, y, ss, ss1, ss2, F, sum, sum1

       pi=4.0*atan(1.0)
       w=cmplx(0, 1)
       alpha=19.0*log(10.0)

       dx=2.0*pi/NN
       dx1=0.5*dx
       N2=0.0*N
       dt=0.01
       ttime=0.0
       Tmax=0.4

c compute the initial fourier coefficient
       do j=0, NN
          x(j)=j*dx
          if ( (x(j)-dx1) .lt. 0.9 .and. (x(j)+dx1) .gt. 0.9) jp=j
       enddo

       do j=0, NN

          s1=cos(0.5*(x(j)+dx1))
          s2=cos(0.5*(x(j)-dx1))

          if (j .lt. jp) then
             ubar(j, 0)=-2.0*(s1-s2)/dx
           elseif (j .gt. jp) then
             ubar(j, 0)=2.0*(s1-s2)/dx
           else
             p1=-2.0*(cos(0.45)-s2)/dx
             p2=2.0*(s1-cos(0.45))/dx
             ubar(j, 0)=p1+p2
          endif
 
          s1=-cos(x(j)+dx1)
          s2=cos(x(j)-dx1)
c          ubar(j, 0)=0.3+0.7*(s1+s2)/dx
       enddo

 99    do NR=0, 2

       do k=0, NN-1
           k1=k-0.5*(NN-1)
           ss=0.0
           do j=0, NN-1
              p1=k1*x(j)
              s1=cos(p1)
              s2=-sin(p1)
              ss=ss+ubar(j,NR)*cmplx(s1,s2)
           enddo
           abar(k)=ss/NN
       enddo

       do k=0, NN-1
           k1=k-0.5*(NN-1)

           if (k1 .eq. 0) then
           si(k)=1.0
           else
           p1=0.5*k1*dx
           p2=sin(p1)
           si(k)=p2/p1
           endif

           al(k)=abar(k)/si(k)

c           ss1=-w*(k+0.5)*0.9
c           ss1=exp(ss1)
c           ss2=-w*(k-0.5)*0.9
c           ss2=exp(ss2)
c           ss1=(-ss1/(k1+0.5)+ss2/(k1-0.5))
c           print *, al(k), ss1, NR, 'ppppppppppppp'
       enddo

c compute the shock location

         N1=sqrt(1.0*N)
         N1=N1*sqrt(1.0*N1)
         N1=(1.0*N)**0.75
c         N1=28
         ss=(N1+1)*al(N+N1+1)/((N1+2)*al(N+N1+2))
         y=-w*log(ss)
         ss=w*(N1+2)*y
         A=w*(N1+2)*al(N+N1+2)*exp(ss)
         if (NR .eq. 2) print *, 'haaaaaaaaaaa', A, y
         if (y1 .gt. 2.*pi .or. y1 .lt. 0.0) then
             y1=0.0
             y=0.0
         endif
         y1=real(y)
         y=real(y)
         A=real(A)

      do j=0, NN-1
         xj=(j+0.5)*dx
          if (xj .le. y1 ) then
             F=-A*xj
          else
             F=A*(2*pi-xj)
          endif
         
          sum=0.0
          do k=1, N
             if ( k .le. N2) then
                 tau=1.0
             else
                 tau=exp(-alpha*(1.0*(k-N2)/N)**6)
             endif

             ps=k*y1
             s1=cos(ps)
             s2=-sin(ps)
             ss1=cmplx(s1, s2)
             ps=k*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=(al(k+N)+w*A*ss1/k)*ss2
             sum=sum+tau*ss

             k1=-k
             ps=k1*y1
             s1=cos(ps)
             s2=-sin(ps)
             ss1=cmplx(s1, s2)
             ps=k1*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=(al(k1+N)+w*A*ss1/k1)*ss2

             sum=sum+tau*ss
          enddo

          uapp1(j)=real(al(N)-A*(pi-y)+sum+F)
c          uapp1(j)=al(N)+sum+F 

       enddo

       do j=0, NN
          if ( (x(j)-dx1) .lt. y1 .and. (x(j)+dx1) .gt. y1) ip=j
       enddo

      do j=0, NN

         sum1=0.0
         xj=x(j)

          if (j .lt. ip) then
             F=-A*xj
          elseif ( j .gt. ip) then
             F=A*(2*pi-xj)
          elseif (j .eq. ip) then
             s1=y1**2-(xj-dx1)**2
             s2=2.0*pi*(xj+dx1-y1)
             s3=(xj+dx1)**2-y1**2
             F=A*(-0.5*s1+s2-0.5*s3)/dx
          endif

          do k=1, N
             if ( k .le. N2) then
                 tau=1.0
             else
                 tau=exp(-alpha*(1.0*(k-N2)/N)**6)
             endif

             ps=k*y1
             s1=cos(ps)
             s2=-sin(ps)
             ss1=cmplx(s1, s2)
             ps=k*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=(abar(k+N)+w*si(k+N)*A*ss1/k)*ss2
             sum1=sum1+tau*ss

             k1=-k
             ps=k1*y1
             s1=cos(ps)
             s2=-sin(ps)
             ss1=cmplx(s1, s2)
             ps=k1*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=(abar(k1+N)+w*A*si(k1+N)*ss1/k1)*ss2
             sum1=sum1+tau*ss
          enddo

          uapp(j)=real(abar(N)-A*(pi-y1)+sum1+F)

       enddo

c Runge-Kutta

       do j=0, NN-1
          ss1=uapp1(j)
          if (j .eq. 0) then
             ss2=uapp1(NN-1)
             else
             ss2=uapp1(j-1)
          endif
          xL(j, NR)=-0.5*(ss1*ss1-ss2*ss2)/dx
c          ubars(j, NR)=uapp(j)
          ubars(j, NR)=ubar(j, NR)

          if (NR .eq. 0) then
             ubar(j,1)=ubars(j,0)+dt*xL(j,0)
          endif

          if (NR .eq. 1) then
             sa=0.75*ubars(j,0)+0.25*ubars(j,1)+0.25*dt*xL(j,1)
             ubar(j,2)=sa
          endif

          if (NR .eq. 2) then
             sa=ubars(j,0)/3.0+ubars(j,2)*2/3.0+dt*xL(j,2)*2/3.0
             ubar(j,3)=sa
          endif

        enddo

       enddo

       ttime=ttime+dt
       print *, ttime, '    time'
       if (ttime .lt. Tmax) then
          do j=0, NN
            ubar(j,0)=ubar(j, 3)
          enddo
          goto 99
       endif

       write (17, *) 'x=['
       do j=0, NN-1
          write (17, *) x(j)
       enddo
       write (17, *) '];'

        write (18, *) 'y=['
       do j=0, NN-1
          write (18, *) real(ubar(j, 3))
       enddo
       write (18, *) '];'

       print *, A, y, '   A, y'
       
       stop
       end