cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc  This code uses the standard RK3, without modifying
cc  the u bar values
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       parameter (N=64, NN=2*N+1)
       implicit real*8(a-h,o-z)
       real x(0:NN), ubar(0:NN,0:3)
       real uapp1(0:NN), xL(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=15.0*log(10.0)

       dx=2.0*pi/NN
       dx1=0.5*dx
       N2=0.0*N
       dt=0.005
       ttime=0.0
       Tmax=2.0
       par=1.9

c compute the initial fourier coefficient
       do j=0, NN
          x(j)=j*dx
          if ( (x(j)-dx1) .lt. par .and. (x(j)+dx1) .gt. par) 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.5*par)-s2)/dx
             p2=2.0*(s1-cos(0.5*par))/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)

       enddo

c compute the shock location

         N1=(1.0*N)**0.6
         ss=(N1+1)*al(N+N1+1)/((N1+2)*al(N+N1+2))
         y=-w*log(ss)
         ss=w*(N1+1)*y
         A=w*(N1+1)*al(N+N1+1)*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)**8)
             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)

       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

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

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

          if (NR .eq. 2) then
             sa=ubar(j,0)/3.0+ubar(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

       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
           al(k)=abar(k)/si(k)
       enddo

      do j=0, NN-1
         xj=x(j)

        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)**4)
             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)
       enddo

       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(uapp1(j))
       enddo
       write (18, *) '];'

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