      subroutine emdlv(r,vp,vs)
c   emstdmdl:
c         set up information on earth model (specified by
c         subroutine call emiasp)
c         set dimension of cpr,rd  equal to number of radial
c         discontinuities in model
      save
      character*(*) name
      character*20 modnam
      dimension cpr(11)
      common/emdlc/np,rd(11)
c
      call mdsmpl(r,rho,vp,vs)
      return
c
      entry emdld(n,cpr,name)
      call mdlget('JB',modnam)
      n=np
      do 1 i=1,np
 1    cpr(i)=rd(i)
      name=modnam
c     print *,'Emdld:  modnam np rd = ',modnam,np,(rd(i),i=1,np)
      return
      end
      subroutine mdlget(mnam,modnam)
c
      implicit double precision (a-h,o-z)
      save
      character*(*) mnam,modnam
      character*3 tag
      character*12 name
      logical log
      real*4 rdd
      dimension w(3,1000)
c   General model.  Allows specification of transverse isotropy and both
c   shear and bulk attenuations.  Elements are:
c     npts number of model points
c     r    radius (km)
c     rho  density (g/cm**3)
c     vpv  velocity of vertically polarized compressional waves (km/s)
c     vph  velocity of horizontally polarized compressional waves (km/s)
c     vsv  velocity of vertically polarized shear waves (km/s)
c     vsh  velocity of horizontally polarized shear waves (km/s)
c     eta  anisotropy parameter
c     qim  inverse q mu
c     qik  inverse q kappa
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
c   Isotropic model.  Elements are:
c     n    number of points in trial model
c     x    radius of trial model (km)
c     vp   velocity of compressional waves (km/s)
c     vs   velocity of shear waves (km/s)
      common/isomod/x(1000),vp(1000),vs(1000),n
c   Polynomial model specification.  Elements are:
c     lyr  number of polynomial segments in radius
c     ri   polynomial normalizing radius (km)
c     rd   outer radius of each segment (km)
c     a    polynomial coefficients
      common/coef/ri,rd(25),a(5,8,25),lyr
c   Logical unit numbers for:
c     in   model input file
c     nout model output file
c     ntrm users listing file
      common/luns/in,nout,ntrm
      common/stdmdl/q(3,1000,3)
      common/emdlc/np,rdd(11)
      data in,nout,ntrm/80,81,10/
c
c   Get the input model file name from the user and open the input model
c   file.
      call query('Model name :',log)
      read(5,*) mnam
      modnam=mnam
*     if(ichar(mnam(1:1)).ge.48.and.ichar(mnam(1:1)).le.57)
*    1 modnam='M'//mnam
      ln=index(mnam,' ')-1
      if(ln.le.0) ln=len(mnam)
      call assign(in,1,'m'//
     1 mnam(1:ln)//'.tmod')
c
c   Read the id line of the input model.
      read(in,103)name,npts,fr,rn
 103  format(5x,a12,i5,d11.3,11x,f9.2)
c   If rn = 0 this is a point-by-point representation.  If rn > 0 this
c   is a polynomial representation.
      iflg=1
      if(rn.gt.0d0) iflg=2
      go to (2,3),iflg
c   Read in the point-by-point model.
 2    read(in,107)tag,r(1),rho(1),vpv(1),vph(1),vsv(1),vsh(1),eta(1),
     1 qim(1),qik(1)
 107  format(1x,a3,f9.2,4f9.4,f8.4,f6.3,2f8.4)
      if(tag.eq.'  1') go to 10
      if(tag.ne.' ') write(ntrm,110)
 110  format(1x,'Mdlget:  inconsistent format.')
      if(tag.eq.' ') write(ntrm,111)
 111  format(1x,'Mdlget:  unexpected file type.')
      call vexit(1)
 10   read(in,107)(tag,r(i),rho(i),vpv(i),vph(i),vsv(i),vsh(i),eta(i),
     1 qim(i),qik(i),i=2,npts)
      go to 4
c   Read in the polynomial coefficients.
 3    call rdpoly(npts,rn)
c   Make a preliminary radial grid.
      call mkgrid
c   Sample the model at the preliminary grid points so that a proper
c   spacing may be determined.
      call mkmodl(iflg)
c   Close the input model file.
 4    call retrns(in)
c   Construct equivalent isotropic velocities.
      call mkiso
c
c   Interpolate the model.
      call drspln(1,npts,r,rho,q(1,1,1),w)
      call drspln(1,npts,r,vp,q(1,1,2),w)
      call drspln(1,npts,r,vs,q(1,1,3),w)
c
c   Find the discontinuities.
      np=0
      do 5 i=2,npts
      if(dabs(r(i)-r(i-1)).gt.1d-7) go to 5
      np=np+1
      rdd(np)=r(i)
 5    continue
      np=np+1
      rdd(np)=r(npts)
      return
      end
      subroutine mdsmpl(r0,ro0,vp0,vs0)
c
      implicit double precision (a-h,o-z)
      save
      real*4 r0,ro0,vp0,vs0
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
      common/isomod/x(1000),vp(1000),vs(1000),n
      common/stdmdl/q(3,1000,3)
c
      rr=r0
      ro0=drsple(1,npts,r,rho,q(1,1,1),rr)
      vp0=drsple(1,npts,r,vp,q(1,1,2),rr)
      vs0=drsple(1,npts,r,vs,q(1,1,3),rr)
      return
      end
      subroutine mkgrid
c
c $$$$$ calls only library routines $$$$$
c
c   Given lyr internal boundaries and their corresponding radii, rd(1) <
c   rd(2) < ... < rd(lyr), mkgrid constructs a spacing which is equal
c   over each interval and allows about 150 km between sucessive
c   samples.
c
      implicit double precision (a-h,o-z)
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
      common/isomod/x(1000),vp(1000),vs(1000),n
      common/coef/ri,rd(25),a(5,8,25),lyr
c   Set up.
      n=0
      r0=0d0
      dr=rd(lyr)/150d0
c   Loop over radial segments.
      do 1 i=1,lyr
c   Determine the number of points in this segment.
      m=max0(int(sngl((rd(i)-r0)/dr+.5d0)),2)-1
c   Set the first point.
      n=n+1
      l=n
      x(n)=r0
      if(m.le.1) go to 2
c   Do the internal points.
      dx=(rd(i)-r0)/m
      do 3 j=2,m
      n=n+1
 3    x(n)=x(n-1)+dx
c   Set the last point.
 2    n=n+1
      x(n)=rd(i)
 1    r0=rd(i)
      return
      end
      subroutine mkiso
c
c $$$$$ calls only library routines $$$$$
c
c   Mkiso combines velocities vpv, vph, vsv, vsh and anisotropy
c   parameter eta at each of npts radial points to form equivalent
c   isotropic velocities vp and vs.  If the medium is already isotropic
c   (vpv = vph, vsv = vsh, and eta = 1) then the result will be vp = vpv
c   = vph and vs = vsv = vsh.
c
      implicit double precision (a-h,o-z)
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts      
      common/isomod/x(1000),vp(1000),vs(1000),n
c   Loop over the discrete radii.
      do 1 i=1,npts
      if(vpv(i).eq.vph(i).and.vsv(i).eq.vsh(i).and.eta(i).eq.1d0)
     1  go to 2
c   Compute the equivalent isotropic velocities if necessary.
      vs(i)=dsqrt(.0666666666666667d0*((1d0-2d0*eta(i))*vph(i)*vph(i)+
     1 vpv(i)*vpv(i)+5d0*vsh(i)*vsh(i)+(6d0+4d0*eta(i))*vsv(i)*vsv(i)))
      vp(i)=dsqrt(.0666666666666667d0*((8d0+4d0*eta(i))*vph(i)*vph(i)+
     1 3d0*vpv(i)*vpv(i)+(8d0-8d0*eta(i))*vsv(i)*vsv(i)))
      go to 1
c   If the velocities at this radius are already isotropic, just copy
c   them into vp and vs.
 2    vs(i)=vsv(i)
      vp(i)=vpv(i)
 1    continue
      return
      end
      subroutine mkmodl(iflg)
c
c $$$$$ calls drsple and drspln $$$$$
c
c   Mkmodl copies grid points from array x into array r and recomputes
c   all model parameters for the new grid points.  If iflg = 1, this is
c   the point-by-point mode and the new values in common /model/ must
c   be determined by b-spline interpolation from the old values.  If
c   iflg = 2, this is the polynomial mode and the new values in common
c   /model/ will be determined directly from the polynomial coefficients
c   a.
c
      implicit double precision (a-h,o-z)
      dimension xmd(1000,8),z(1000),w(1000,8),q(3,1000,8)
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
      common/isomod/x(1000),vp(1000),vs(1000),n
      common/coef/ri,rd(25),a(5,8,25),lyr
      common/luns/in,nout,ntrm
c   Xmd is just a convienent way to manipulate the model parameters in
c   common /model/.
      equivalence (xmd,rho)
      data eps/5d-3/
c   Set up some counters.
      n1=n-1
      m=1
      i=0
      k=0
      go to (1,2),iflg
c
c   Point-by-point mode.  Spline everything.
c
 1    do 3 l=1,8
 3    call drspln(1,npts,r,xmd(1,l),q(1,1,l),w)
c   Copy the old radial grid and the old model parameters into scratch
c   storage.
      do 4 j=1,npts
      z(j)=r(j)
      do 4 l=1,8
 4    w(j,l)=xmd(j,l)
c   Set up some additional counters.
      nn1=npts-1
      mm=1
c   Loop over radial segments bounded by sucessive identical radial
c   samples.
 5    if(k.ge.n) go to 6
c   Find the next boundary in the new grid.
      do 7 k=m,n1
      if(dabs(x(k+1)-x(k)).le.eps) go to 8
 7    continue
      k=n
c   Find the next boundary in the old grid.
 8    if(mm.gt.nn1) go to 22
      do 9 kk=mm,nn1
      if(dabs(z(kk)-z(kk+1)).le.eps) go to 10
 9    continue
      kk=npts
c   Copy model parameters from the old grid to the new for the point
c   just above the lower boundary.
 10   r(m)=z(mm)
      do 11 l=1,8
 11   xmd(m,l)=w(mm,l)
c   Interpolate all model parameters for new internal samples in the
c   segment.
      m1=m+1
      k1=k-1
      if(m1.gt.k1) go to 14
      do 12 j=m1,k1
      r(j)=x(j)
      do 12 l=1,8
 12   xmd(j,l)=drsple(1,npts,z,w(1,l),q(1,1,l),x(j))
c   Copy model parameters from the old grid to the new for the point
c   just below the upper boundary.
 14   r(k)=z(kk)
      do 15 l=1,8
 15   xmd(k,l)=w(kk,l)
c   Go do the next segment.
      m=k+1
      mm=kk+1
      go to 5
c
c   Polynomial mode.  Loop over radial segments.
c
 2    if(k.ge.n) go to 6
c   I is the segment counter.
      i=i+1
c   Find the upper boundary of the current segment.
      do 17 k=m,n1
      if(dabs(x(k+1)-x(k)).le.eps) go to 18
 17   continue
      k=n
c   Determine all model parameters at all points in this segment using
c   the appropriate polynomials.
 18   do 19 j=m,k
      r(j)=x(j)
      s=ri*x(j)
      do 19 l=1,8
 19   xmd(j,l)=a(1,l,i)+s*(a(2,l,i)+s*(a(3,l,i)+s*(a(4,l,i)+
     1 s*a(5,l,i))))
c   Go do the next segment.
      m=k+1
      go to 2
c   Fill in defaults.
 6    do 20 j=1,k
c   If only the vertically or the horizontally polarized p velocity is
c   given assume isotropy.
      if(vpv(j).le.0d0) vpv(j)=vph(j)
      if(vph(j).le.0d0) vph(j)=vpv(j)
c   Do the same for s velocities.
      if(vsv(j).le.0d0) vsv(j)=vsh(j)
      if(vsh(j).le.0d0) vsh(j)=vsv(j)
c   If eta isn't set, assume isotropy.
      if(eta(j).le.0d0) eta(j)=1d0
 20   continue
c   Check that the whole model was used.
      if(iflg.ne.1.or.kk.eq.npts) go to 21
 22   write(ntrm,100)
 100  format(1x,'Mkmodl:  inconsistent grids.')
      call exit
c   Update the number of radii.
 21   npts=k
c   Check that all segments were used.
      if(iflg.eq.1.or.(iflg.ge.2.and.i.eq.lyr)) return
      write(ntrm,101)
 101  format(1x,'Mkmodl:  inconsistent number of layers.')
      call exit
      end
      subroutine mkstep(max,fmax,ppw,itop,ibot,m,ftru)
c
c $$$$$ calls drsple $$$$$
c
c   Subroutine mkstep invents an m point node spacing between radii of
c   r(ibot) < r(itop).  Mkstep tries to make:
c   2*(z(i)-z(i-1))/(vv(i)+vv(i-1))=c(i) <= 1000./(fmax*ppw),
c   i=2,3,...,m.  On exit ftru = max(1000./(c(i)*ppw)), i=1,2,...,m for
c   the last iteration.  This technique produces the minimum number of
c   nodes per model such that there are at least ppw nodes per
c   wavelength for the highest frequency mode of interest (fmax mhz).
c   Max is the most nodes that mkstep may produce.
c
      implicit double precision (a-h,o-z)
      dimension v(1000),x(1000),va(1000)
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
      common/fakvel/vv(1000),qv(3,1000),z(1000)
      data sm,tol,noi/.2d0,3d-5,10/
c   Working from the top of the layer down, use the velocity, vv, at the
c   beginning of each interval to guess the interval length.
      no=0
      fac=1d3/(fmax*ppw)
      z(1)=r(itop)
      v(1)=vv(itop)
      i=2
      z(2)=z(1)-.5d0*fac*v(1)
      go to 10
 1    if(i.ge.max) go to 15
      i=i+1
      z(i)=z(i-1)-fac*v(i-1)
 10   v(i)=drsple(1,npts,r,vv,qv,z(i))
      if(z(i).gt.r(ibot)) go to 1
c   See if one more point is desirable.
      z(i)=.5d0*(z(i)+z(i-1))
      if(z(i).le.r(ibot)) go to 15
c   Add one more point.
 16   z(i)=.333333333333333d0*(z(i-1)+2d0*r(ibot))
      v(i)=drsple(1,npts,r,vv,qv,z(i))
      i=i+1
      go to 9
 15   if(i.le.2) go to 2
c   See if the last point is necessary.
      if(z(i-1)-r(ibot).ge.sm*(z(i-2)-z(i-1))) go to 9
      i=i-1
      go to 16
c   Add the bottom point of the layer.
 9    m=i
      z(m)=r(ibot)
      v(m)=vv(ibot)
      m1=m-1
      if(m.ge.3) go to 8
c   Guarantee that there are at least two points per layer.
 2    i=3
      z(2)=.5d0*(r(itop)+r(ibot))
      v(2)=drsple(1,npts,r,vv,qv,z(2))
      go to 9
c   Calculate the length and mean velocity of each interval.
 8    do 3 i=2,m
      x(i)=z(i)-z(i-1)
 3    va(i)=.5d0*(v(i)+v(i-1))
c   Make half intervals on each end.
      va(m)=.5d0*va(m)
 11   va(2)=.5d0*va(2)
c   Ensure that the sum of the interval lengths will equal the layer
c   thickness.
      c=0d0
      do 4 i=2,m
 4    c=c+va(i)
      c=(r(ibot)-r(itop))/c
      del=0d0
c   Calculate new interval lengths.
      do 5 i=2,m
      dx=c*va(i)-x(i)
      del=del+dx*dx
 5    x(i)=x(i)+dx
c   Construct a new spacing.
      do 6 i=2,m1
 6    z(i)=z(i-1)+x(i)
c   Exit if the rms change to the spacing was negligible or if noi
c   iterations have been done.
      del=sqrt(del/m1)
      no=no+1
      if(del.le.tol.and.no.le.noi) go to 12
c   Resample the velocity, vv, on the new spacing.
      do 7 i=2,m1
      v(i)=drsple(1,npts,r,vv,qv,z(i))
 7    va(i)=.5d0*(v(i)+v(i-1))
      va(m)=.25d0*(v(m)+v(m1))
      go to 11
c   The nodes will be positioned at the center of each interval found by
c   mkstep.
 12   m=m-1
      m1=m-1
      do 14 i=2,m1
 14   z(i)=z(i)+.5d0*x(i+1)
      z(m)=z(m+1)
      ftru=-1d3/(c*ppw)
      return
      end
      subroutine rdpoly(nlyr,rn)
c
c $$$$ calls no other routine $$$$$
c
c   Rdpoly reads the coefficients specifying a polynomial mode model.
c   Radial segments must appear in order of increasing radius.
c   Polynomials for all model parameters defined must immediately
c   follow (one per line) the corresponding "layer" line specifying the
c   upper radius of the segment.  Polynomials for each segment are
c   identified by a character tag and may be in any order.  Polynomials
c   which are not defined are assumed to be all zeros.
c
      implicit double precision (a-h,o-z)
      character*3 tag,ids(8)
c6    integer*4 tag,ids(8)
      logical fnd(8)
      dimension w(5)
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
      common/coef/ri,rd(25),a(5,8,25),lyr
      common/luns/in,nout,ntrm
c   These (upper case only) tags are recognized.
      data ids/'RHO','VPV','VPH','VSV','VSH','ETA','AMU','AKA'/
c   Allow up to max segments.
      data max/25/
c   Ri, the inverse normalizing radius, is used in evaluating the
c   polynomials.
      ri=1d0/rn
c   Read the first line.
      read(in,100,end=15,err=17)tag,w
 100  format(1x,a3,5f9.4)
c   It is an error if the first line is not a "layer" line.
      if(tag.ne.' ') go to 13
      i=0
c   Loop over all nlyr segments.  A zero upper segment boundary
c   terminates the input.
 1    if(w(1).le.0d0) go to 14
c   Too many layers is an error.
      if(i.ge.max) go to 16
c   Set the segment upper radius.
      i=i+1
      rd(i)=w(1)
c   Initialize found flags to false and all polynomials to zeros.
      do 2 j=1,8
      fnd(j)=.false.
      do 2 k=1,5
 2    a(k,j,i)=0d0
c   Read the next line.
 3    read(in,100,end=14,err=17)tag,w
c   If it is a "layer" line, go initialize the next segment.
      if(tag.eq.' ') go to 1
c   Determine which model parameter is involved.
      do 4 j=1,8
      if(tag.eq.ids(j)) go to 5
 4    continue
c   It is an error if the tag is unknown.
      write(ntrm,101)tag
 101  format(1x,'Rdpoly:  polynomial tag ',a,' is not defined.')
      call exit
c   It is an error if this tag has already been found in this segment.
 5    if(fnd(j)) go to 18
c   Set the found flag and the polynomial.
      fnd(j)=.true.
      do 6 k=1,5
 6    a(k,j,i)=w(k)
      go to 3
c   All done.  See if all segments were indeed found.
 14   lyr=i
      if(lyr.eq.nlyr) go to 7
c   Issue a warning if the number of segments found is different from
c   the number on the model id line.
      write(ntrm,102)lyr,nlyr
 102  format(1x,'Rdpoly:',i4,' layers found.',i4,' expected.')
 7    if(rd(lyr).eq.rn) return
c   Issue a warning if the normalizing radius is different from the free
c   surface radius.
      write(ntrm,103)rd(lyr),rn
 103  format(1x,'Rdpoly:  free surface radius (',f6.1,
     1 ') is inconsistent with normalizing'/10x,'radius (',f6.1,').')
      return
c   Issue various error messages and exit.
 13   if(tag.ne.'  1') write(ntrm,104)
 104  format(1x,'Rdpoly:  second line must be a layer radius.')
      if(tag.eq.'  1') write(ntrm,110)
 110  format(1x,'Rdpoly:  unexpected model type.')
      call exit
 15   write(ntrm,105)
 105  format(1x,'Rdpoly:  no layers have been defined.')
      call exit
 16   write(ntrm,106)max
 106  format(1x,'Rdpoly:  attempt to define more than',i3,' layers.')
      call exit
 17   write(ntrm,107)i
 107  format(1x,'Rdpoly:  bad format in layer',i3,'.')
      call exit
 18   write(ntrm,108)ids(j),i
 108  format(1x,'Rdpoly:  attempt to redefine ',a,' in layer',i3,'.')
      call exit
      end
      subroutine regrid(fmax,ppw,ftru)
c
c $$$$$ calls drspln and mkstep $$$$$
c
c   Subroutine regrid breaks a model up into layers bounded by first
c   order discontinuities (flagged by adjacent, identical radii).  After
c   mkstep determines the optimal node spacing, regrid copies it into
c   array x.  On return, the n points of x will be the optimized grid
c   for a new version of the model.  Regrid ensures that at all
c   frequencies smaller than fmax (mhz), there will be at least ppw
c   nodes per wavelength radially.  Ftru is the highest frequency (mhz)
c   that actually has the requisite sampling everywhere.
c
      implicit double precision (a-h,o-z)
      dimension w(3,1000)
      common/model/r(1000),rho(1000),vpv(1000),vph(1000),vsv(1000),
     1 vsh(1000),eta(1000),qim(1000),qik(1000),npts
      common/isomod/x(1000),vp(1000),vs(1000),n
      common/fakvel/vv(1000),qv(3,1000),z(1000)
      data max,eps/1000,5d-3/
      rfac=(r(npts)-r(1))/r(npts)
      ftry=fmax/rfac
c   Vv holds the minimum (non-zero) local velocity.
      do 2 j=1,npts
      vv(j)=vs(j)
      if(vv(j).le.0d0) vv(j)=vp(j)
 2    continue
      call drspln(1,npts,r,vv,qv,w)
c   Set up for a pass through the model.
 7    n=0
      k=1
      n1=npts-1
      ftru=1d38
c   Step through the model looking for discontinuities.
      do 1 j=1,n1
      if(dabs(r(j+1)-r(j)).gt.eps.and.j.lt.n1) go to 1
      i=j
      if(j.eq.n1) i=npts
c   A discontinuity has been found, make up the node spacing and return
c   it in array z.
      call mkstep(max,ftry,ppw,i,k,m,fgot)
c   Fgot is the highest frequency actually achieved.
      ftru=dmin1(ftru,fgot)
c   If there are too many nodes don't transfer them to x.
      if(n+m.le.max) go to 6
      n=n+m
      go to 5
c   Transfer points z to x.
 6    nn=m+1
      do 3 l=1,m
      n=n+1
 3    x(n)=z(nn-l)
 5    k=i+1
 1    continue
      ftru=ftru*rfac
c   Return the maximum frequency that has the desired number of nodes
c   per wavelength everywhere.
      if(n.le.max) return
c   If there are too many nodes readjust the frequency limit.
      ftry=ftry*(max*.99d0/n+.01d0)
      go to 7
      end
      double precision function drsple(i1,i2,x,y,q,s)
c
c C$$$$ calls only library routines $$$$$
c
c   Drsple returns the value of the function y(x) evaluated at point s
c   using the cubic spline coefficients computed by drspln and saved in
c   q.  If s is outside the interval (x(i1),x(i2)) drsple extrapolates
c   using the first or last interpolation polynomial.  The arrays must
c   be dimensioned at least - x(i2), y(i2), and q(3,i2).
c
      implicit double precision (a-h,o-z)
      dimension x(1),y(1),q(3,1)
      data i/1/
      ii=i2-1
c   Guarantee i within bounds.
      i=max0(i,i1)
      i=min0(i,ii)
c   See if x is increasing or decreasing.
      if(x(i2)-x(i1))1,2,2
c   X is decreasing.  Change i as necessary.
 1    if(s-x(i))3,3,4
 4    i=i-1
      if(i-i1)11,6,1
 3    if(s-x(i+1))5,6,6
 5    i=i+1
      if(i-ii)3,6,7
c   X is increasing.  Change i as necessary.
 2    if(s-x(i+1))8,8,9
 9    i=i+1
      if(i-ii)2,6,7
 8    if(s-x(i))10,6,6
 10   i=i-1
      if(i-i1)11,6,8
 7    i=ii
      go to 6
 11   i=i1
c   Calculate drsple using spline coefficients in y and q.
 6    h=s-x(i)
      drsple=y(i)+h*(q(1,i)+h*(q(2,i)+h*q(3,i)))
      return
      end
      subroutine drspln(i1,i2,x,y,q,f)
c   Rspln computes cubic spline interpolation coefficients
c   for y(x) between grid points i1 and i2 saving them in q.  The
c   interpolation is continuous with continuous first and second
c   derivatives.  It agrees exactly with y at grid points and with the
c   three point first derivatives at both end points (i1 and i2).
c   X must be monotonic but if two successive values of x are equal
c   a discontinuity is assumed and seperate interpolation is done on
c   each strictly monotonic segment.  The arrays must be dimensioned at
c   least - x(i2), y(i2), q(3,i2), and f(3,i2).  F is working storage
c   for rspln.
c                                                     -rpb
      implicit double precision (a-h,o-z)
      dimension x(1),y(1),q(3,1),f(3,1),yy(3)
      equivalence (yy(1),y0)
      data yy/3*0.d0/
      j1=i1+1
      y0=0.d0
c   Bail out if there are less than two points total.
      if(i2-i1)13,17,8
 8    a0=x(j1-1)
c   Search for discontinuities.
      do 3 i=j1,i2
      b0=a0
      a0=x(i)
      if(a0-b0)3,4,3
 3    continue
 17   j1=j1-1
      j2=i2-2
      go to 5
 4    j1=j1-1
      j2=i-3
c   See if there are enough points to interpolate (at least three).
 5    if(j2+1-j1)9,10,11
c   Only two points.  Use linear interpolation.
 10   j2=j2+2
      y0=(y(j2)-y(j1))/(x(j2)-x(j1))
      do 15 j=1,3
      q(j,j1)=yy(j)
 15   q(j,j2)=yy(j)
      go to 12
c   More than two points.  Do spline interpolation.
 11   a0=0.d0
      h=x(j1+1)-x(j1)
      h2=x(j1+2)-x(j1)
      y0=h*h2*(h2-h)
      h=h*h
      h2=h2*h2
c   Calculate derivative at near end.
      b0=(y(j1)*(h-h2)+y(j1+1)*h2-y(j1+2)*h)/y0
      b1=b0
c   Explicitly reduce banded matrix to an upper banded matrix.
      do 1 i=j1,j2
      h=x(i+1)-x(i)
      y0=y(i+1)-y(i)
      h2=h*h
      ha=h-a0
      h2a=h-2.d0*a0
      h3a=2.d0*h-3.d0*a0
      h2b=h2*b0
      q(1,i)=h2/ha
      q(2,i)=-ha/(h2a*h2)
      q(3,i)=-h*h2a/h3a
      f(1,i)=(y0-h*b0)/(h*ha)
      f(2,i)=(h2b-y0*(2.d0*h-a0))/(h*h2*h2a)
      f(3,i)=-(h2b-3.d0*y0*ha)/(h*h3a)
      a0=q(3,i)
 1    b0=f(3,i)
c   Take care of last two rows.
      i=j2+1
      h=x(i+1)-x(i)
      y0=y(i+1)-y(i)
      h2=h*h
      ha=h-a0
      h2a=h*ha
      h2b=h2*b0-y0*(2.d0*h-a0)
      q(1,i)=h2/ha
      f(1,i)=(y0-h*b0)/h2a
      ha=x(j2)-x(i+1)
      y0=-h*ha*(ha+h)
      ha=ha*ha
c   Calculate derivative at far end.
      y0=(y(i+1)*(h2-ha)+y(i)*ha-y(j2)*h2)/y0
      q(3,i)=(y0*h2a+h2b)/(h*h2*(h-2.d0*a0))
      q(2,i)=f(1,i)-q(1,i)*q(3,i)
c   Solve upper banded matrix by reverse iteration.
      do 2 j=j1,j2
      k=i-1
      q(1,i)=f(3,k)-q(3,k)*q(2,i)
      q(3,k)=f(2,k)-q(2,k)*q(1,i)
      q(2,k)=f(1,k)-q(1,k)*q(3,k)
 2    i=k
      q(1,i)=b1
c   Fill in the last point with a linear extrapolation.
 9    j2=j2+2
      do 14 j=1,3
 14   q(j,j2)=yy(j)
c   See if this discontinuity is the last.
 12   if(j2-i2)6,13,13
c   No.  Go back for more.
 6    j1=j2+2
      if(j1-i2)8,8,7
c   There is only one point left after the latest discontinuity.
 7    do 16 j=1,3
 16   q(j,i2)=yy(j)
c   Fini.
 13   return
      end
