PROGRAM BurgersTX

c 
c BurgersTX.f
c Modified April 2004 to correct array bounds error
c     reported by Kent Gee at Penn State
c Version 1.1 Modified September 1999
c Obtained from http://people.bu.edu/robinc/kzk
c
c Robin Cleveland
c Department of Aerospace and Mechanical Engineering
c Boston University
c Boston, MA 02215
c robinc@bu.edu
c
c Mark Hamilton
c Department of Mechanical Engineering 
c The University of Texas at Austin
c Austin, TX 
c hamilton@mail.utexas.edu
c
c   The fortran source code for a time-domain computer code developed at
c The University of Texas at Austin to to solve a very general form of 
c the Burgers equation in the time domain. It accounts for nonlinear
c distortion, thermoviscous absorption, absorption and dispersion due 
c to multiple relaxation processes, geometrical spreading, and a 
c stratified (inhomogeneous) medium.
c
c  The principal references for this code are:
c
c R. O. Cleveland, M. F. Hamilton, and D. T. Blackstock,  "Time-domain
c modeling of finite-amplitude sound in relaxing fluids," J. Acoust. Soc.
c Am., Vol. 99, 3312-3318 (1996).
c
c R. O. Cleveland, "Propagation of sonic booms through a real stratified 
c atmosphere," Ph.D. dissertation, The University of Texas at Austin (1995).
c
c   A variant of this code (THOR) was used to model sonic boom propagation 
c in the atmosphere:
c
c Robin O. Cleveland, James P. Chambers, Henry E. Bass, Richard
c Raspet, Mark F. Hmailton, David T. Blackstock, "Comparison
c of computer codes for the propagation of sonic booms through
c the atmosphere", J. Acoust. Soc. Am., vol 97, pp 906--917 (1995).
c
c   The algorithm for the Burgers code was based on code to solve the KZK
c equation (the KZK code is also available from this web site).   The 
c references for the KZK code are:
c 
c Y.-S. Lee, "Numerical solution of the KZK equation for pulsed finite
c amplitude sound beams in thermoviscous fluids," Ph.D. dissertation, The
c University of Texas at Austin (1993).
c 
c Y.-S. Lee and M. F. Hamilton, "Time-domain modeling of pulsed
c finite-amplitude sound beams," J. Acoust. Soc. Am., Vol. 97, pp. 906-917
c (1995).
c
c   V1.1   September 1999
c      Made the input routine more consistent.  However it is no longer
c        compatible with old input files.
c      Fixed an error in the internal waveform generation algorithm
c   V 1.0   February 1996
c   Beta version created December 1995
c
c  Code Development
c
c  KZKpulse is the mother of all time domain codes
c     A code for solving the KZK equation in the time domain. Accounts for
c     nonlinearity, diffraction, and thermoviscous effects.
c
c  THOR - father of BurgersTX
c     This code was developed to model sonic boom propagation in the
c     atmosphere. Routines were added to cope with multiple relaxation
c     processes. The diffraction was removed and in its place routines
c     were added to account for geometrical spreading, and a stratified
c     (inhomogeneous) medium.
c
c  BurgersTX 
c    This code is really a refinement of THOR. The basic algorithm is 
c    the same and sonic boom specific routines were removed. The major 
c    variation is that spreading and stratification are accounted for 
c    by scaling the pressure and range variables. Also all variables are 
c    nondimensionalised by the source conditions.  BurgersTX can 
c    internally generate simple sinusoidal pulses. 
c
c  The output from file is a large array of the form
c     N      sigma0       sigma1     ... sigmaM
c   tau0  P(sig0,tau0)  P(sig1,tau0) ... P(sigM,tau0)
c   ...
c   tauN  P(sig0,tauN)  P(sig1,tauN) ... P(sigM,tauN)
c
c   where N is the number of time points
c   sigma0 and P(sig0,tau) is all the initial location 
c      and waveform
c
c  For an arbitrary initial waveform the file must be of the form:
c     N     
c   P(0) 
c   P(1)
c   ...
c   P(N)
c
c  The code will squeeze the N points into the number of cycles
c  you specify in the input file
c
c

       implicit none
       integer outfile, rtfile, paramfile, mdmfile, wavefile, maxTauPts
       integer maxRelax, maxMdm
       integer maxRanges
       real fudge
       real dSnl, dStv, dSmr, dSsp, dSmd
       real pi
       common pi

c These parameters are used to set the file numbers, dimension the 
c arrays, and control the step size of the code.
c
       parameter (outfile=10)
       parameter (mdmfile=11)
       parameter (paramfile=12)
       parameter (rtfile=13)
       parameter (wavefile=14)
       parameter (maxTauPts=65536)
       parameter (maxRanges=50)
       parameter (maxRelax=5)
       parameter (maxMdm=10)
       parameter (fudge=1e-7)
       parameter (dSnl=0.2)
       parameter (dStv=0.1)
       parameter (dSmr=0.1)
       parameter (dSsp=0.05)
       parameter (dSmd=0.02)

c VARIABLE DEFINITIONS
c   The following is a fairly complete list of variable definitions.
c   note that most of these vars are only really used in Exec.
c
c    nlflag  = flag which is 1 if nonlinear effects are on
c    tvflag  = flag which is 1 if thermoviscous effects are on
c    mrflag  = flag which is 1 if relaxation effects are on
c    mdflag = flag which is 1 for an inhomogeneous medium, that is,
c                the user provides medium data in another file.
c
c    P	    	= array holding REAL pressure vs. tau and rho 
c                 at a specific sigma
c    tau	= array giving tau values
c
c    NumTauPts 	= number of divisions along tau axis
c    SigMax	= maximum range sigma for which calculations 
c                 are desired
c    sigma 	= range variable
c    sigma0     = initial range variable
c    dSigma	= step size in sigma direction
c    dSigCN	= INITIAL sigma step, Crank-Nicolson method
c    dTau	= grid step size in tau direction

c    spread	= variable controlling spreading

c    scale    = scaling function to make the amplitude physical
c    phinl    = phi function for the nonlinear process
c    TV       = phitv/Gamma - where phitv os the phi function for tv
c    numRelax = number of relaxation processes
c    RxD      = array containing the Dnu*phi_rc for relaxation process
c    RxT      = array with thetanu/phi_rt for relaxation

c    NumRanges   = number of ranges at which beam patterns are desired
c    watchPoints = array containing points at which to output data.
c    Target     = range of the next output point
c    TargetID   = index of next output point
c    toTarget   = distance left to go to next ouput point
c    ThisRange     = the desired pattern range at ThisRangeID.
 
c    dSnl, dStv, dSmr, dSsp, dSmd  these parameters control the maximum step
c                 size that each effect (nonlinearity, thermoviscous, 
c                 relaxation, spreading, and medium (straification) 
c                 respectively) will let the code take. See Robin's disertation
c                 Sec. 4.6 for a more complete description.

c    rtstep	= step between different runtime calculations. Runtime is 
c             a subroutine that lets the user output data on the fly

c    Note that the variables are declared HERE only to reserve memory.
c    How the memory is actually used is determined by the subroutines
c    and how the variables are redefined there.
c
c    The variables declared here are mainly for workspace and to
c    pre-allocate storage space.
c
       integer nlflag, tvflag, mrflag, mdflag
       integer numTauPts, numRelax, numRanges
       real dTau, sigma, dSigCN, SigMax, rtstep, spread
       real scale, phinl, TV
       real RxD(1:maxRelax), RxT(1:maxRelax)
       real WatchPoints(1:maxRanges)

c WORKSPACES
c for Exec:
       real P(0:maxTauPts), Ptarget(1:maxTauPts,0:maxRanges)
       real tau(0:maxTauPts)
c for TriDiag:
       real a(0:maxTauPts),b(0:maxTauPts),c(0:maxTauPts)
       real rhs(0:maxTauPts)
       real gamma1(0:maxTauPts),gamma2(0:maxTauPts)
c for Nonlinear
c      use a and b from TriDiag for Pold and t 
c for Medium
       integer numMdm
       real mdmlast(1:maxMdm), mdmnext(1:maxMdm), mdmddz(1:maxMdm)
c

       pi=3.14159265358979323841217556

       Call Initialise (maxTauPts, maxRelax, maxRanges, maxMdm,
     +               outfile, P, tau, dTau, numTauPts, wavefile,
     +               scale, phinl, TV, RxD, RxT, numRelax,
     +	             mdmlast, mdmnext, mdmddz, numMdm, mdmfile,
     +               sigma, dSigCN, sigMax, rtstep, spread,
     +               fudge, numRanges, WatchPoints, dSnl, dStv, 
     +               nlflag, tvflag, mrflag, mdflag, paramfile)

       Call Exec  (P, Ptarget, tau, dTau, numTauPts,
     +              scale, phinl, TV, RxD, RxT, numRelax,
     +	            mdmlast, mdmnext, mdmddz, numMdm, mdmfile,
     +              sigma, dSigCN, sigMax, rtstep, spread, fudge,
     +              numRanges, WatchPoints, outfile, rtfile,
     +              a,b,c,rhs,gamma1,gamma2,
     +              nlflag, tvflag,mrflag, mdflag, paramfile,
     +              dSnl, dStv, dSmr, dSmd)

       stop
       end
c
       SUBROUTINE Initialise (maxTauPts, maxRelax, maxRanges, maxmdm, 
     +                   outfile, P, tau, dTau, numTauPts, wavefile,
     +                   scale, phinl, TV, RxD, RxT, numRelax,
     +	                 mdmlast, mdmnext, mdmddz, numMdm, mdmfile,
     +                   sigma, dSigCN, sigMax, rtstep, spread, 
     +                   fudge, numRanges, WatchPoints, dSnl, dStv,
     +                   nlflag, tvflag, mrflag, mdflag, paramfile)

c
c This routine accepts user input to set up the programme and
c sets all the required variables
c
       implicit none

c subroutinue's input arguments
       integer maxTauPts,maxRelax,maxRanges,maxMdm,outfile,wavefile

c subroutine's output arguments
       real P(0:maxTauPts), tau(0:maxTauPts)
       real scale, phinl, TV
       real RxD(1:maxRelax), RxT(1:maxRelax)
       real mdmlast(1:maxMdm), mdmnext(1:maxMdm), mdmddz(1:maxMdm) 
       integer mdmfile, numMdm
       integer numTauPts, numRelax, numRanges
       real dTau, sigma, dSigCN, sigMax, rtstep, spread, fudge
       real Watchpoints(1:maxRanges), dSnl, dStv
       integer nlflag, tvflag, mrflag, mdflag, paramfile

c internal variables
       character*40 waveform, filename
       integer loadflag,mexpo
       real cycles, taumin, taumax
       real pi, rTmp1, rTmp2, cyclepts
       common pi
       integer i
       
       write(*,*) 'Burgers Pulse Programme'

       write(*,*) 'Nonlinear effects (1=yes):'
       read(*,*) nlflag
       write(*,*) 'Thermoviscous effects (1=yes):'
       read(*,*) tvflag
       write(*,*) 'Relaxation effects (1=yes):'
       read(*,*) mrflag
       write(*,*) 'Inhomogeneous medium (1=yes):'
       read(*,*) mdflag

       if (mdflag.EQ.0) then
         write(*,*) 'Spreading effects (0:plane, 0.5:cylindrical,',
     +		  	' 1.0:spherical) :'
         read(*,*) spread
       endif

       write(*,*) 'LOAD the initial waveform (1) or calculate it (0)?'
       read(*,*) loadflag

       if (loadflag.EQ.0) then	
          write(*,*) 'ENTER the Exponent (m) for the envelope:'
          read(*,*) mexpo
          write(*,*) 'ENTER the number of cycles (n) in the envelope:'
          read(*,*) cycles
          write(*,*) 'RISE-TIME CYCLES = ',(cycles/mexpo)

          write(*,*) 'ENTER number of cycles to pad before the wave:'
          read(*,*) rTmp1
          tauMin=-pi*(2.0*rTmp1+cycles)

          write(*,*) 'ENTER number of cycles to pad after the wave:'
          read(*,*) rTmp2
          tauMax=pi*(cycles+2.0*rTmp2)

          write(*,*) 'ENTER the number of points per cycle:'
          read(*,*) cyclepts
          dTau=2*pi/cyclepts
          numTauPts=(rTmp1+cycles+rTmp2)*cyclepts

         if (numTauPts.GT.maxTauPts) then
           write(*,*) 'ERROR in ReadData - too many tau points.'
           write(*,*) 'User requires ',numTauPts,' but maximum is ',
     +                   maxTauPts
           stop
         endif

          Call MakeWave(P, tau, numTauPts, dTau,
     +                   tauMin, cycles, mexpo)

        else
          write(*,*) 'ENTER the file name of the input data:'
          read(*,'(A40)') waveform
          write(*,*) waveform
          write(*,*) 'ENTER the number of cycles (n) in the waveform:'
          read(*,*) cycles
          write(*,*) 'ENTER number of cycles to pad before the wave:'
          read(*,*) rTmp1
          write(*,*) 'ENTER number of cycles to pad after the wave:'
          read(*,*) rTmp2

          write(*,*) 'ENTER the number of points per cycle:'
          read(*,*) cyclepts

          open(unit=wavefile,file=waveform,status='old')
          numTauPts=maxTauPts
          call ReadWave(P, tau, numTauPts, dTau, cycles,
     +                   rTmp1, rTmp2, wavefile)
          close(wavefile)

        endif

c       write(*,*) 'ENTER the file for the output waveforms:'
c       read(*,'(A40)') filename
c       open(unit=outfile,file=filename,status='unknown')
       open(unit=outfile,file='pulse.out',status='unknown')
       write(outfile,'(i13,$)') numTauPts

 7      format(1x,a,i3,1x,a)

        if (mdflag .EQ. 1) then
          write(*,*) 'ENTER the file containing the medium data:'
          read(*,'(A30)') filename
       	  open(unit=mdmfile,file = filename, status = 'old')

          read(mdmfile,*) numMdm
          numRelax=(numMdm-4)/2

          read(mdmfile,*) (mdmlast(i),i=1,numMdm)
          read(mdmfile,*) (mdmnext(i),i=1,numMdm)

          rTmp1=mdmnext(1)-mdmlast(1)
          do i=2,numMdm
            mdmddz(i)=(mdmnext(i)-mdmlast(i))/rTmp1
          enddo

          scale=mdmlast(2)
          phinl=mdmlast(3)
          TV=mdmlast(4)
          do i=1,numRelax
            RxD(i)=mdmlast(i*2+3)
            RxT(i)=mdmlast(i*2+4)
          enddo
	  	
        else
          scale=1.0
          phinl=1.0

          write(*,*) 'ENTER the TV attenuation A (=1/Gamma):'
          read(*,*) TV

       	  write(*,*) 'ENTER the number of relaxation processes desired:'
          write(*,7) '[max=',maxRelax,']'
          read(*,*) numRelax

          if (numRelax.GT.maxRelax) numRelax=maxRelax

          write(*,7) 'ENTER the',numRelax,'relaxation processes.'
          write(*,*) '  Dnu   thetanu '
          do 70 i=1,numRelax
            read(*,*) RxD(i), RxT(i)
 70       continue
        endif

        write(*,*) 'ENTER the start position :'
        read(*,*) sigma

        write(*,*) ' '
        write(*,*) 'ENTER the number of ranges desired:'
        write(*,7) '[max=',maxRanges,']'
        read(*,*) numRanges
        if (numRanges.LT.1) numRanges=1
        if (numRanges.GT.maxRanges) numRanges=maxRanges

        write(*,7) 'ENTER the ',numRanges,
     +               'sigma range(s), in ASCENDING order.'
        sigMax = sigma
        do 90 i=1,numRanges
 85       read(*,*) rTmp1
          if (i.EQ.1) then
            rTmp2=0.0
            if (rTmp1.LT.sigma) rTmp2=sigma
          endif
          rTmp1=rTmp1+rTmp2
          if (rTmp1.LT.sigMax) then
            write(*,*) 'Enter values in ASCENDING order. RE-ENTER:'
            goto 85
          endif
c adjust by small amount to guard against computational
c roundoff-error when making comparisons against these values.
          watchpoints(i) = rTmp1-fudge
 90     continue
        sigMax = rTmp1-fudge

c Get the step size
        write(*,*) 'ENTER step size :'
        read(*,*) dSigCN

        rtstep=0.0
        write(*,*) 'ENTER step between runtime calcs:'
        read(*,*,END=95) rtstep

95      open(unit=paramfile, file = 'Parameters', status='unknown')
        write(paramfile,305) 'Nonlinear:',dSnl*nlflag,'  TV:',
     +            dStv*tvflag,'  Relaxation:',mrflag,'  Medium:',mdflag
 305	format(a,f6.3,a,f6.3,a,i2,a,i2)

       if (loadflag.EQ.1) then
         write(paramfile,310) waveform, numTauPts
       else
         write(paramfile,*) 'cycles=',cycles, ' m=',mexpo
       endif
 310    format('Waveform:',a18,' numTauPts=',i5)

        if (mdflag.EQ.1) write(paramfile,320) filename
 320      format(' Medium:',a18)

        write(paramfile,330) sigma, dSigCN, spread
 330    format('Sigma=',f6.2,' step=',1p,e7.1,' spread=',0p,f3.1)

        write(paramfile,340) TV, numRelax
 340    format(' TV=',1p,e9.3,' numRelax=',i2)
        write(paramfile,*) 
     +             '    Dnu      thetanu      r/ds      q'
        do 350 i=1, numRelax
          write(paramfile,'(1p,e11.3,e11.3,e11.3,e11.3,e8.2)') Rxd(i), 
     +       RxT(i), RxD(i)*RxT(i)/(2*dTau*dTau),RxT(i)/(2*dTau) 
350     continue

        write(paramfile,360)
 360    format('Output points:',$)
        write(paramfile,380) ( watchpoints(i), i=1, numRanges )
 380    format(1p,e11.3,$)

        write(paramfile,*) ' '
        close (paramfile)

        return
        end
c
c
c
       SUBROUTINE MakeWave(P, tau, NumTauPts, dTau, tauMin, 
     +                       cycles, mexpo)

c This subroutine initializes the P and tau arrays.

       implicit none
       integer numTauPts
       real P(0:numTauPts), tau(0:numTauPts)
       real dTau, cycles, tauMin
       integer mexpo
       integer i
       real pi
       common pi

c VARIABLE DEFINITIONS: (see Main)

       do 500 i=0,numTauPts
         tau(i) = tauMin + i*dTau
         P(i)=sin(tau(i))
         if (mexpo.GT.0) then
            P(i)=P(i)*exp(-(tau(i)/(cycles*pi))**(2*mexpo))
         elseif (abs(tau(i)).GT.pi*cycles) then
            P(i) = 0
         endif
500    continue

       return
       end
c
c
c
       subroutine ReadWave(P, tau, numTauPts, dTau, cycles,
     +                       before, after, wavefile)

c This subroutine initializes the P and tau arrays.

       implicit none
       integer numTauPts
       real P(0:numTauPts), tau(0:numTauPts)
       real dTau, cycles, before, after
       integer wavefile

       real tauMin, ppc, maxPts
       integer i, iTmp
       real pi
       common pi

       maxPts=numTaupts
       read(wavefile,*) numTauPts

       ppc=float(numTauPts)/cycles
       numTauPts=(before+cycles+after)*ppc
       if (numTauPts.GT.maxPts) then
         write(*,*) 'ERROR in ReadData - too many tau points.'
	 write(*,*) 'User requires ',numTauPts,' but maximum is ',maxPts
         stop
       endif

       dTau=2.0*pi/ppc
       tauMin=-pi*(2.0*before+cycles)
       iTmp=nint((before+cycles)*ppc)

       read(wavefile,*) P(0)
       tau(0)=tauMin
       do i=1,nint(before*ppc)+1
         P(i)=P(0)
         tau(i)=tauMin+i*dTau
       end do

       do i=nint(before*ppc)+2,iTmp
         tau(i)=tauMin+i*dTau
         read(wavefile,*) P(i)
       end do

       do i=iTmp+1,numTauPts
         P(i)=P(iTmp)
         tau(i)=tauMin+i*dTau
       end do

      return
      end
c
c
c
       SUBROUTINE Exec (P, Ptarget, tau, dTau, numTauPts,
     +                   scale, phinl, TV, RxD, RxT, numRelax,
     +                   mdmlast, mdmnext, mdmddz, numMdm, mdmfile,
     +                   sigma, dSigCN, sigMax, rtstep, spread, 
     +                   fudge, numRanges, WatchPoints,
     +                   outfile, rtfile, 
     +                   a,b,c,rhs,gamma1,gamma2,
     +                   nlflag, tvflag, mrflag, mdflag, paramfile,
     +                   dSnl, dStv, dSmr, dSmd)

c
c This routine is responsible for the main execution of the program.
c It steps through sigma space and then outputs the waveforms.
c
       implicit none

c passed arguments
       integer numTauPts, numRelax, numMdm, numRanges
       real P(0:numTauPts), Ptarget(1:numTauPts,0:numRanges)
       real tau(0:numTauPts)
       real scale, phinl, TV
       real RxD(1:numRelax), RxT(1:numRelax)
       real mdmlast(1:numMdm), mdmnext(1:numMdm), mdmddz(1:numMdm)
       integer mdmfile
       real dTau, sigma, dSigCN, sigMax, rtstep, spread, fudge
       real watchpoints(1:numRanges)
       integer outfile, rtfile
       integer nlflag, tvflag, mrflag, mdflag, paramfile
       real dSnl, dStv, dSmr, dSmd

c internal vars.
       real thisdS, nextdS, dSigma
       real Target
       integer TargetID
       integer i,j
       real nextrt, totarget
       real dpmax, sigma0
c The following variables may be used to caculate the runtime
       real nlt, mdt, mrt, tvt, rtt, oht
c      real dtime_, dum(2) 

c sub-WORKSPACES:
c for TriDiag:
       real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts)
       real rhs(0:numTauPts)
       real gamma1(0:numTauPts),gamma2(0:numTauPts)
c for Nonlinear
c      use a and b from TriDiag for Pold and t

c variable definitions for internal vars.
c     thisdS	= sigma step size at this step
c     nextdS    = sigma step for next step
c     Target	= next sigma targetted for output by user
c     TargetID  = index of the next target in the watchpoints array
c

       dpmax=0.0
       do i=1,numTauPts
         Ptarget(i,0) = P(i)
         dpmax=max(dpmax,P(i) - P(i-1))
       enddo

       sigma0=sigma

       if (rtstep.GT.0.0) then
         nextrt=sigma
       else
         nextrt=sigMax+fudge
       endif
      
       TargetID=1
       Target = watchpoints(1)
       write(outfile,'(1p,e13.5,$)') sigma
       write(*,*) 'Saved waveform at sigma=',sigma
       totarget=Target-sigma

c Step up initial step size
       nextdS=totarget+fudge
       if (dSigCN.GT.0) nextdS=min(nextdS,dSigCN)
       if (mdflag.EQ.1) nextdS=min(nextdS,mdmnext(1)-mdmlast(1))
       if (spread.GT.0) nextdS=min(nextdS,dSmd*sigma0/spread)
       if (tvflag.GT.0) nextdS=min(nextdS,dStv/(TV*numTauPts))
       if (mrflag.GT.0) then
         do i=1,numRelax
           dSigma=dSmr*(1+numTauPts*RxT(i)**2)/(RxD(i)*numTauPts)
     +              *min(1.0,6.28/(sqrt(float(numTauPts))*RxT(i)))
           nextdS=min(nextdS,dSigma)
         enddo
       endif
       if (nlflag.GT.0) nextdS=min(nextdS,dSnl*dTau/(dpmax*phinl))

       mrt=0.0
       mdt=0.0
       nlt=0.0
       tvt=0.0
       rtt=0.0
c      write(*,*) 'Startup time',dtime_(dum)
       oht=0.0


11000  thisdS=min(nextdS,totarget+fudge)
       nextdS=10.0*thisdS
       if (dSigCN.GT.0) nextdS=min(nextdS,dSigCN)
c      oht=oht+dtime_(dum)

       if (sigma.GE.nextrt) then
          dSigma=thisdS
          call Runtime(P, sigma, dSigma, numTauPts, spread,
     +                  mdmlast, mdmnext, mdmddz, numMdm, rtfile,
     +                  scale, phinl, TV, RxD, RxT, numRelax,
     +	                nlt, mdt, mrt, tvt, rtt, oht,mdflag)
          nextrt=sigma+rtstep
c     	  rtt=rtt+dtime_(dum)
       endif

       if (mdflag.EQ.1) then
         dSigma=thisdS
         call Medium(P, numTauPts, dSigma, sigma, sigma0,
     +               mdmlast, mdmnext, mdmddz, numMdm, mdmfile,
     +               scale,phinl,TV,RxD,RxT,numRelax,mdflag,dSmd)
         nextdS=min(nextdS,dSigma)
c        mdt=mdt+dtime_(dum)
       elseif (spread.GT.0) then
         scale=(sigma0/sigma)**spread
         phinl=scale
         nextdS=min(nextdS,dSmd*sigma/spread)
c        mdt=mdt+dtime_(dum)
       endif

       if (tvflag.GT.0) then
          dSigma=thisdS
          call CNFDtv (P, numTauPts, dSigma, dTau, TV, 
     +                        a,b,c,rhs,gamma1,gamma2, dStv)
          nextdS=min(nextdS,dSigma)
c         tvt=tvt+dtime_(dum)
        endif

       if (mrflag.GT.0) then
         do i=1,numRelax
           dSigma=thisdS
           call CNFDrelax (P, numTauPts,dSigma,dTau,RxD(i),RxT(i),
     +       a,b,c,rhs,gamma1,gamma2,dSmr)
           nextdS=min(nextdS,dSigma)
         enddo
c        mrt=mrt+dtime_(dum)
       endif

       if (nlflag.GT.0)  then
         dSigma=thisdS
         call Nonlinear (P, numTauPts, dSigma, sigma, tau, phinl,
     +                     a, b, dSnl, paramfile)
         nextdS=min(nextdS,dSigma)
c        nlt=nlt+dtime_(dum)
         if (dSnl.EQ.-1) then
           sigMax=sigma
           Target=sigma
         endif
       endif

       sigma = sigma + thisdS

       totarget=Target-sigma

11150  if (totarget.LE.0) then
       write(outfile,'(1p,e13.5,$)') sigma

c  Store waveform to be output later
         do 11200 i=1,numTauPts
           Ptarget(i,TargetID) = P(i)*scale
11200    continue

       write(*,*) 'Saved waveform at sigma=',sigma,'  dsigma=',dsigma

         if (TargetID.LE.numRanges) TargetID = TargetID + 1
         Target = watchpoints(TargetID)
         totarget=Target-sigma

       endif

       if (sigma.LT.sigMax) go to 11000

c       write(*,*) 'Execution time for the subroutines'
c       write(*,*) 'Runtime ',rtt
c       write(*,*) 'Medium ',mdt
c       write(*,*) 'Thermoviscous ',tvt
c       write(*,*) 'Relaxation ',mrt
c       write(*,*) 'Nonlinear',nlt
c       write(*,*) 'Overhead ',oht
c       write(*,*) 'Total ',oht+nlt+mrt+tvt+mdt+rtt

c Write output
       write(outfile,*) ' '
       do 11310 i=1,numTauPts
         write(outfile,'(1p,e13.5,$)') tau(i)
         do 11300 j=0,TargetID-2
           write(outfile,'(1p,e13.5,$)') Ptarget(i,j)
11300    continue
         write(outfile,'(1p,e13.5)') Ptarget(i,TargetID-1)
11310  continue
       close(outfile)

       return
       end
c
c
c
       SUBROUTINE CNFDtv (P, numTauPts, dSigma, dTau,
     +                        TV, a,b,c,rhs,gamma1,gamma2, dStv)
c
c This routine performs the Crank-Nicolson finite difference method 
c for the absorption.  
c
       implicit none
       integer numTauPts
       real P(0:numTauPts)
       real dSigma, dTau, TV, dStv
       real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts)
       real rhs(0:numTauPts)
       real gamma1(0:numTauPts),gamma2(0:numTauPts)
       integer start, finish
       real r
       integer i

c  r = a convenient constant relating dSigma and dTau
c  a,b,c,rhs,soln = temporary workspace used for call to Tridiag

       start = 0
       finish = numTauPts

c initialize the coefficient "bands" in the tridiagonal matrix

       r = TV*dSigma/(2.0*dTau**2)
       a(start)=0.0
       b(start)=1.0
       c(start)=0.0
       do 30010 i=start+1,finish-1
         a(i) = -r
         b(i) = (1.0 + 2.0*r)
         c(i) = -r
30010  continue
       c(finish)=0.0
       b(finish)=1.0
       a(finish)=0.0

       rhs(start)=P(start)
       do 30020 i=start+1,finish-1
         rhs(i) = r*P(i-1) + (1.0-2.0*r)*P(i) + r*P(i+1)
30020  continue
       rhs(finish)=P(finish)

       call Tridiag(start,finish,a,b,c,rhs,P,gamma1,gamma2,numTauPts)
	
       dSigma=dStv/(TV*numTauPts)

       return
       end
c
c
c
       SUBROUTINE CNFDrelax (P, numTauPts, dSigma, dTau, RxD,
     +                 RxT, a,b,c,rhs,gamma1,gamma2,dSmr)
c
c This routine performs the Crank-Nicolson finite difference method
c for the absorption due to molecular relaxation.
c
       implicit none
       integer numTauPts
       real P(0:numTauPts)
       real dSigma, dTau, RxD, RxT, dSmr
       real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts)
       real rhs(0:numTauPts)
       real gamma1(0:numTauPts),gamma2(0:numTauPts)
       integer start, finish
       real r, q
       integer i

c  r, q = convenient constants relating dSigma and dTau
c  a,b,c,rhs,soln = temporary workspace used for call to Tridiag

       start = 1
       finish = numTauPts - 1

c initialize the coefficient "bands" in the tridiagonal matrix

       q = RxT/(2.0*dTau)
       r = RxD*RxT*dSigma/(2.0 * dTau**2.0)

       a(start)=0.0
       b(start)=1.0
       c(start)=0.0
       do i=start+1,finish-1
         a(i) = -r - q
         b(i) = (1.0 + 2.0*r)
         c(i) = q - r
       end do
       c(finish)=0.0
       b(finish)=1.0
       a(finish)=0.0

       rhs(start)=P(start) 
       do i=start+1,finish-1
         rhs(i)=(r-q)*P(i-1) + (1.0-2.0*r)*P(i) + (r+q)*P(i+1)
       end do
       rhs(finish)=P(finish)

       call Tridiag(start,finish,a,b,c,rhs,P,gamma1,gamma2,numTauPts)

       dSigma=dSmr*(1+numTauPts*RxT**2)/(numTauPts*RxD)
       dSigma=dSigma*min(1.0,6.28/(sqrt(float(numTauPts))*RxT))

       return
       end
c
c
c
       SUBROUTINE  Nonlinear (P, numTaupts, dSigma, sigma, tau, 
     +				phinl, Pold, t, dSnl, paramfile)

c
c This routine adds the nonlinear effects to the P array.
c The "t" array is used as temporary workspace to check for shocks.
c
       implicit none
       integer numTauPts
       real P(0:numTauPts)
       integer paramfile
       real sigma, dSigma, phinl, dSnl
       real tau(0:numTauPts)
       real Pold(0:numTauPts), t(0:numTauPts)
       real dpdt, dpm, dpm2
       integer i, j


c VARIABLE DECLARATIONS
c	t		= temp copy of tau array
c	Pold		= old pressure
c
c distort the phase of the wavelet and check that continuity 
c is maintained.
       Pold(0)=P(0)
       Pold(1)=P(1)
       t(1) = tau(1)-P(1)*dSigma*phinl
       do 60060 i=2,(numTauPts-1)
	 Pold(i)=P(i)
         t(i) = tau(i) - P(i)*dSigma*phinl
         if (t(i).LE.t(i-1)) go to 62000
60060  continue
       Pold(numTauPts)=P(numTauPts)

	dpm=0.0
	dpm2=0.0
	i=1
	do j=1,numTauPts
	  do while (tau(j).GT.t(i))
	    i=i+1
	    if (i.GT.numTauPts) then
	      do i=j,numTauPts
	        P(i)=P(numTauPts)
	      enddo
	      go to 60100
	    endif
            dpdt=(Pold(i)-Pold(i-1))/(t(i)-t(i-1))
            if (dpdt.GT.dpm) dpm=dpdt
	  enddo
	  P(j)=Pold(i-1)+dpdt*(tau(j)-t(i-1))
	enddo

60100   dSigma=dSnl/(dpm*phinl)

       return

c The discontinuity error message will halt the program
 
62000  write(*,*) 'Discontinuity: smaller dSigma step required'
       write(*,*) 'sigma=',sigma,' dsigma=',dSigma,' phinl=',phinl
       write(*,*) 'i   ',i-1,i,i+1
       write(*,*) 'P   ',P(i-1), P(i), P(i+1)
       write(*,*) 'tau ',tau(i-1), tau(i), tau(i+1)
       write(*,*) 't   ',t(i-1), t(i), t(i+1)

       open(unit=paramfile,file='Parameters',status='unknown')
       write(paramfile,*) 
     +         'ERROR: Waveform discontinuity at sigma = ', sigma
       write(paramfile,*) 'Position in array, i = ',i
       close(paramfile)

       dSnl=-1
       return
       end
c
c
c
       SUBROUTINE Tridiag(start, finish, a,b,c, rhs,
     +                     soln, gamma1, gamma2, n)

c      This routine solves the tridiagonal system of (FINISH-START+1)
c      equations.  The arrays A,B,C are the diagonals of the matrix
c      defining the tridiagonal system.  START and FINISH are indices
c      for the beginning and ending corners that define a SUBMATRIX
c      within the original matrix.  The routine solves the submatrix
c      using the THOMAS ALGORITHM and stores the results into the
c      SOLN array.
c
c      GAMMA1 and GAMMA2 are used as temporary workspace.
c
       implicit none
       integer start, finish, n, i
       real a(0:n), b(0:n), c(0:n), rhs(0:n), soln(0:n)
       real gamma1(0:n), gamma2(0:n)

       gamma1(start) = b(start)
       gamma2(start)= rhs(start)/gamma1(start)

       do 90000 i=(start+1),finish
         gamma1(i) = b(i) - a(i)*c(i-1)/gamma1(i-1)
         gamma2(i) = (rhs(i) - a(i)*gamma2(i-1))/gamma1(i)
90000  continue

       soln(finish) = gamma2(finish)
       do 90020 i= (finish-1),start,-1
         soln(i) = gamma2(i) - c(i)*soln(i+1)/gamma1(i)
90020  continue

       return
       end

c
c
c
       SUBROUTINE Medium(P, numTauPts, sigma, sigma0, dSigma,
     +                  mdmlast, mdmnext, mdmddz, numMdm, mdmfile, 
     +	                scale, phinl, TV, RxD, RxT, numRelax,
     +                  mdflag, dSmd)

c      This routine changes the absorption coefficients to take into
c     account of an inhomegeneuos atmopshere.
       implicit none

       integer numTauPts
       real P(0:numTauPts)
       integer numMdm, mdmfile, numRelax, mdflag
       real dSigma, sigma, sigma0, dSmd
       real mdmlast(1:numMdm), mdmnext(1:numMdm), mdmddz(1:numMdm)
       real scale, phinl, TV, RxD(1:numRelax), RxT(1:numRelax)

c Internal variables
       real oldscale, dSigMD	
       integer i

       oldscale=scale
       do while (sigma.GT.mdmnext(1))
         do i=1, numMdm
           mdmlast(i)=mdmnext(i)
         enddo
         read(mdmfile,*,end=90120) (mdmnext(i),i=1,numMdm)
       enddo

       dsigMD=mdmnext(1)-mdmlast(1)
       do i=2, numMdm
         mdmddz(i)=(mdmnext(i)-mdmlast(i))/dsigMD
       enddo
       goto 90130
90120  mdflag=-1

90130  dsigMD=sigma-mdmlast(1)
       scale=mdmlast(2)+dsigMD*mdmddz(2)
       phinl=mdmlast(3)+dsigMD*mdmddz(3)
       TV = mdmlast(4)+dsigMD*mdmddz(4)
       do i=1, numRelax
         RxD(i)=mdmlast(2*i+3)+dsigMD*mdmddz(2*i+3)
         RxT(i)=mdmlast(2*i+4)+dsigMD*mdmddz(2*i+4)
       enddo
	 
       dSigma=dSmd*dSmd/(abs(1.0-scale/oldscale))

       return
       end
c
c
c
       SUBROUTINE Runtime(P, sigma, dSigma, numTauPts, spread,
     +                  mdmlast, mdmnext, mdmddz, numMdm, rtfile,
     +                  scale, phinl, TV, RxD, RxT, numRelax,
     +                  nlt, mdt, mrt, tvt, rtt, oht, mdflag)

c      This routine allows the user to output various important values
c	as the programme is running.
       implicit none

       integer numTauPts
        real P(0:numTauPts)
        integer numMdm, rtfile, numRelax, mdflag
        real dSigma, sigma, spread
        real mdmlast(1:numMdm), mdmnext(1:numMdm), mdmddz(1:numMdm)
        real scale, phinl, TV, RxD(1:numRelax), RxT(1:numRelax)
        real nlt, mdt, mrt, tvt, rtt, oht

        real rt, pmax, pmin, gPlo, gPhi, lPlo, invm, P10, P90, dummy
        integer i, imax, imin, iglo, ighi, illo, rtx
 
       pmax=0.0
       pmin=0.0
       
       do 91020 i=1,numTauPts
         if (P(i).GT.pmax) then
           pmax=P(i)
           imax=i
         endif
         if (P(i).LT.pmin) then
           pmin=P(i)
           imin=i
         endif
91020  continue
       
       open(rtfile,file='runtime',status='unknown')
91025  read(rtfile,*,END=91026) dummy
       goto 91025
91026  write(rtfile,*) sigma, pmax, imax, pmin, imin
       close(rtfile)

c        write(rtfile,'(1p,e11.3,$)') sigma
c        write(rtfile,'(1p,e11.3,$)') dSigma
c        write(rtfile,'(1p,e11.3,$)') scale
c        write(rtfile,'(1p,e11.3,$)') phinl
c	write(rtfile,*) ' '

c	   rtx=0
c   Calculate 10% to 90% ristime of the tallest shock
c       if (rtx.EQ.1) then
c         iglo=0
c         gPlo=P(0)
c   	     ighi=0
c         gPhi=P(0)
c	     illo=0
c	     lPlo=P(0)
c         do i=1,numTauPts
c           if (P(i).GT.gPhi) then
c c            ighi=i
c	         gPhi=P(i)
c             if (illo.NE.iglo) then
c	           illo=iglo
c	           lPlo=gPlo
c           endif
c	       if (P(i).LT.gPlo) then
c             iglo=i
c             gPlo=P(i)
c           endif
c	     enddo

c	     P10=lPlo+0.1*(gPhi-lPlo)
c	     P90=lPlo+0.9*(gPhi-lPlo)

c         i=illo
c         do while(P(i).LT.P10)
c      	    i=i+1
c           if (i.GT.min(ighi,numTauPts)) then
c             rt=0
cc             go to 91100
c           endif
c         enddo
c         invm=1.0/(numTauPts*(P(i)-P(i-1)))
c         rt=invm*(P10-P(i-1)) + float(i-1)/numTauPts

c         i=ighi
c         do while (P(i).GT.P90)
c           i=i-1
c           if (i.LT.max(illo,1)) then
c             rt=0
c             go to 91100
c           endif
c         enddo
cc         invm=1.0/(numTauPts*(P(i+1)-P(i)))
cc         rt=invm*(P90-P(i)) + float(i)/numTauPts - rt

c91100    write(rtfile,'(1p,e13.5,$)') rt
c       endif

c  Write execution time for the various algorithms
c       write(rtfile,'($,1p,e13.5)') rtt
c       write(rtfile,'($,1p,e13.5)') mdt
c       write(rtfile,'($,1p,e13.5)') tvt
c       write(rtfile,'($,1p,e13.5)') mrt
c       write(rtfile,'($,1p,e13.5)') nlt
c       write(rtfile,'($,1p,e13.5)') oht

c       if (mdflag.EQ.1) then
c         write(rtfile,*)
c     +          (mdmlast(i)+mdmddz(i)*(sigma-mdmlast(1)), i=2,numMdm)
c       endif

c       close(rtfile)
c       write(*,*) 'Runtime closed2'
       
       return
       end