**********************************************************************
*                                                                    *
       Program Burgers
	   
       complex   k1(100),p(100),p2(100),a(100),k2,j
       complex   sum,absrel,absbl
       real      xmax,tv,omega,d,b,pi
       integer   nhar,count,nout
       common    tv,d,omega,b

*      Get inputs:
       write(*,*) "***********************************************"
       write(*,*) "*    Solution of Spectral Burgers Equation:   *"
       write(*,*) "*             Propagation Curves              *"
       write(*,*) "***********************************************"
       write (*,*) "Enter sigma step size (.01)."
       read(*,*) step
       write (*,*) "Enter number of steps per output."
       read(*,*) nout
       write(*,*) 'Enter number of harmonics in calculation (50).'
       read(*,*) nhar
       write(*,*) "Enter the max sigma range for output."
       read(*,*) xmax
       write(*,*) "Enter the thermoviscous loss coef."
       read(*,*) tv
 
*      Define some constants:
       x = 0.0
       count = 0
       i = 1
       j = cmplx(0.0,1.0)
       omega = 1.0
       pi = 2.0*asin(1.0)

*      Zero out p(n) and get the absorption values:
       do 100 n=1, nhar
         p(n)=cmplx(0,0)
         k1(n)=cmplx(0,0)
         p2(n)=cmplx(0,0)
         a(n)=abstv(n)
  100  continue

*      Initial condition -- Here we start with a pure sinusoid:
       p(1)=cmplx(0.0,-1.0)

*      Integrating with Second order Runge-Kutta:
       write(*,*) 'sigma , fund  , 2nd  ,  3rd, 4th , 5th'
       do 500 i=1, xmax/step+1
         x = x + step
         count = count + 1
         do 300 n=1, nhar
           k1(n)=step*((j*n/4.0)*sum(p,n,nhar)-a(n)*p(n))
           p2(n)=p(n)+k1(n)
  300    continue
         do 400 n=1, nhar
           k2=step*((j*n/4.0)*sum(p2,n,nhar)-a(n)*p2(n))
           p(n)=p(n)+(k1(n)+k2)/2.0
  400    continue
         if (count.eq.nout) then
           write(*,900) x,',', cabs(p(1)),',', cabs(p(2)),',', 
     +     cabs(p(3)),',', cabs(p(4)),',', cabs(p(5))
           count = 0
         end if
  500  continue
  900  format(5(e20.8,A1,2X),e20.8,2X)

       end
*                                                                    *
**********************************************************************
**********************************************************************
*      Subroutine for nonlinearity (discrete convolution)

       complex function sum(p,n,nhar)

       complex p(100)
       integer n,nhar

       sum = cmplx(0.0,0.0)
       do 1400 m=1, n-1
         sum = sum + p(m)*p(n-m)
 1400  continue
       do 1410 m=nhar, n+1, -1
         sum = sum + 2.0*p(m)*conjg(p(m-n))
 1410  continue

       return
       end
*                                                                    *
**********************************************************************
**********************************************************************
*      Absorption coefficients and dispersion relations:
*      Thermoviscious absorption

       real function abstv(n)
       
       common  tv, d, omega, b
       real tv,d,omega,b
       integer n
 
       abstv = (n**2)*tv

       return
       end
*                                                                    *
**********************************************************************
**********************************************************************
*      Absorption coefficients and dispersion relations:
*      thermoviscous absorption and relaxation

       complex function absrel(n)

       complex j, arelax
       common  tv,d,omega,b

       j = cmplx(0.0,1.0)
       relax = d*omega
       atv = (n**2)*tv
       arelax = (n**2)*relax*(1.0-j*n*omega)/(1.0+(n*omega)**2)
       absrel = atv + arelax

       return
       end
*                                                                    *
**********************************************************************
**********************************************************************
*      Absorption coefficients and dispersion relations:
*      Thermoviscous absorption and boundary layer effects

       complex function absbl(n)

       complex j
       common  tv,d,omega,b

       j = cmplx(0.0,1.0)
       atv = (n**2)*tv
       absbl = atv + (1.0+j)*b*sqrt(real(n))

       return
       end
*                                                                    *
**********************************************************************