! these subroutines are from old modvelden.f program

      subroutine GETVEL(zz,z,v,nvz,vv)
      real z(*),v(*)
      do 10 i=2,nvz
         if (z(i).gt.zz) go to 20
10    continue
      i=nvz
20    dz=(z(i)-z(i-1))
      if (dz.eq.0.) then
         print *,'***Error in GETVEL, zz,i = ',zz,i
         stop
      end if
      frac=(zz-z(i-1))/dz
      vv=v(i-1)+frac*(v(i)-v(i-1))
      return
      end
 
c RAYTRACE does raytracing for v(s) model
c  Inputs:  z_s   =  array with depth points (dimension 700)
c           v_s   =  array with velocities (dimension 700)
c           nvz   =  number of z(v) points
c           pmin  =  minimum ray parameter (flat earth)
c           pmax  =  maximum ray parameter
c           np    =  number of ray parameters (max 50)
c           zmax  =  maximum depth
c  Returns: p     =  array with ray parameters
c           x     =  700x50 array with range offset (km) two-way to surface 
c           t     =  700x50 array with time (minutes) two-way to surface
c           tstar =  700x50 array with tstar values two-way to surface
c           x_tot =  x50 array with total range (km)
c           t_tot =  x50 array with total time (minutes)
c           tstar_tot = x50 array with total tstar
c  NOTE: both up and downgoing legs are counted.  tstar values are from PREM
c        Q model assuming S-waves.  
      subroutine RAYTRACE(z_s,v_s,nvz,pmin,pmax,np,zmax,
     &              p,x,t,tstar,x_tot,t_tot,tstar_tot)
      parameter (nz0=900)  
      real z_s(nz0),v_s(nz0),z(nz0),v(nz0),slow(nz0)
      real p(50),t(nz0,50),x(nz0,50),tstar(nz0,50)
      real t_tot(50),x_tot(50),tstar_tot(50),q(nz0)
      real kmdeg

      if (nvz.gt.nz0) then
         print *,'***ERROR in RAYTRACE: nvz,nz0 = ',nvz,nz0
      end if

      erad=6371.
      pi=3.1415927
      ecircum=2.*pi*erad
      kmdeg=ecircum/360.
      degrad=180./pi

      do 10 i=1,nvz
         call FLATTEN(z_s(i),v_s(i),z(i),v(i))
         if (v(i).ne.0.) slow(i)=1./v(i)
         if (i.eq.1) go to 10
         zmid=(z_s(i)+z_s(i-1))/2.         
         if (zmid.lt.80.) then
            q(i)=600.
         else if (zmid.lt.220.) then
            q(i)=80.
         else if (zmid.lt.670.) then
            q(i)=143.
         else
            q(i)=312.
         end if
c         if (i.gt.nvz-5) print *,i,z_s(i),z(i)
10    continue

20    dp=(pmax-pmin)/float(np-1)
      do 50 ip=1,np
         p(ip)=pmin+dp*float(ip-1)
         x_tot(ip)=0.
         t_tot(ip)=0.
         tstar_tot(ip)=0.
         x(1,ip)=0.
         t(1,ip)=0.
         imth=3
         do 40 i=1,nvz-1
            if (z_s(i).ge.zmax) go to 50
            h=z(i+1)-z(i)
            call LAYERTRACE(p(ip),h,slow(i),slow(i+1),imth,dx,dt,irtr)
            dx=dx/kmdeg
            dt=dt/60.
            x_tot(ip)=x_tot(ip)+2.*dx
            t_tot(ip)=t_tot(ip)+2.*dt
            x(i+1,ip)=x_tot(ip)
            t(i+1,ip)=t_tot(ip)
            if (q(i+1).ne.0.) then
               tstar_tot(ip)=tstar_tot(ip)+120.*dt/q(i+1)
               tstar(i+1,ip)=tstar_tot(ip)
            end if
c            if (ip.eq.1.and.i.gt.nvz-5) then
c               print *,i,h,slow(i),slow(i+1),dx,dt,irtr
c            end if
            if (irtr.eq.0.or.irtr.eq.2) go to 50
40       continue
50    continue

100   return
      end
c
c
c CRUSTOP convolves seismogram with crustal operator for surface reflection
c
c  Inputs:  a      =  input trace
c           n      =  number of points in a
c           dt     =  sample interval (s)
c           dtmoho =  two-way surface-to-Moho travel time
c           rmoho  =  Moho reflection coefficient
c  Returns: b  =  convolved trace
c
      subroutine CRUSTOP(a,b,n,dt,dtmoho,rmoho)
      real a(n),b(n)
      tmoho=(1.-rmoho**2)
      ioff=nint(dtmoho/dt)
c      print *,'ioff = ',ioff
      do 100 i=1,n
         b(i)=a(i)*tmoho**2
         i1=i-ioff
         if (i1.ge.1) b(i)=b(i)-a(i1)*rmoho
         i2=i+ioff
         if (i2.le.n) b(i)=b(i)+a(i2)*rmoho
100   continue
      return
      end

c
c
c GAUSSOP convolves seismogram with normal distribution function
c  Inputs:  a      =  input trace
c           n      =  number of points in a
c           dt     =  sample interval (s)
c           sd     =  standard deviation of normal distribution
c  Returns: b      =  convolved trace
c
      subroutine GAUSSOP(a,b,n,dt,sd)
      real a(n),b(n),btot(10000)
c      real a(*)
      pi=3.1415927
      fact=1./sqrt(2.*pi)
      do 10 i=1,n
         b(i)=0.
         btot(i)=0.
10    continue
      nlim=nint(5.*sd/dt)     !truncate at 5 sd
      do 50 ii=-nlim,nlim
         t=float(ii)*dt
         x=t/sd
         f=fact*exp(-0.5*x**2)
         do 40 i=1,n
            j=i+ii
            if (j.lt.1.or.j.gt.n) go to 40
            b(j)=b(j)+a(i)*f
            btot(j)=btot(j)+f
40       continue
50    continue
      do 60 i=1,n
         if (btot(i).eq.0.) go to 60
         b(i)=b(i)/btot(i)
60    continue

      return
      end
