! program to interactively explore mantle V(Z) models

program playmod
   implicit none
   
   integer, parameter :: nvz0 = 1000, np0 = 5000
   
   integer :: i, nvz, nvzref, iw, np, iscr, ii, iocean, icrustop, nboot, nseis, &
              i10, n1200, ioff, igeocor, ip, iref1, iref2, ithou, iz, n_syn, &
              nplot, i410, i560, menu, np2
   real :: erad, pi, ecircum, kmdeg, degrad, dfact, dfactinv, r, vp, vs, &
           xw1a, xw1b, xw2a, xw2b, xw3a, xw3b, yw1a, yw1b, yw2a, yw2b, yw3a, yw3b, &
           xp1a, xp1b, xp2a, xp2b, xp3a, xp3b, yp1a, yp1b, yp2a, yp2b, yp3a, yp3b, &
           xtic1, xtic2, ytic1, ytic2, x, y, rfirstmin, rfirstmax, dist, distmin, &
           dz, dx, xmin, xmax, delavg, delfix, qdepmin, qdepmax, stnmin, stnmax, wexp, &
           z10, zmax, pmin, pmax, ttest0, vv1, vv2, refl, tmoho1, tmoho2, &
           x_ss1, t_ss1, x_ss2, t_ss2, frac, dp, dxdp, pavg, u0, s1, s2, t_ss0, &
           tref, tref2, tref1, zz1, zz2, den1, den2, zavg, ravg, psph, the1, the2, &
           f1, f2, refl2, geocor, standev, r_earth, x_ss, t_ss, tthou, xthou, &
           velden, xx, yy, x1, x2, tpred
   
   real, dimension(nvz0) :: zref, z, den, vref, v 
   real, dimension(np0) :: deltab, tttab
   real, dimension(100) :: tfirstref, tfirst, rfirst
   real, dimension(-600:300) :: time, amp, ampmin, ampmax, syn1, syn2, syn3, &
                                xsyn, psyn, gsyn, zsyn
   real, dimension(900) :: z_syn, v_syn, p_syn, t_tab(900), p_tab(900), &
                           z_examp, v_examp, g_tab, den_syn
   real, dimension(900,50) :: x_syn, t_syn, x_comb, t_comb, tstar, x_iasp, t_iasp 
   real, dimension(50) :: x_syn0, t_syn0, tstar_tot, x_iasp0, t_iasp0, p_iasp   
   
   character (len = 100) :: vzfile, ssfile
   
   erad = 6371.
   pi = 3.1415927
   ecircum = 2.*pi*erad
   kmdeg = ecircum/360.
   degrad = 180./pi

   r_earth = erad
   igeocor = 1
!   standev = 1.5               !2.5 used in Shearer 1996
   nplot = 1

   print *, 'Enter reference V(Z) file name'
   read (*, '(a)') vzfile
   
   print *, 'Enter wavetype: (1) P   or   (2) S'
   read *, iw
   
!   open (13, file='out.debug')

! read reference velocity model file   
   open (11, file=vzfile, status='old')   
   do i = 1, nvz0
      read (11, *, end=10) zref(i), r, vp, vs, den(i)
      if (vs == 0.) go to 10
      if (iw == 1) then
         vref(i) = vp
      else
         vref(i) = vs
      endif
      z(i) = zref(i)      
      v(i) = vref(i)
   enddo
   i = nvz0 + 1
10 close (11)
   nvzref = i - 1
   nvz = nvzref
   print *, 'nvzref = ', nvzref

! ray trace model to get predicted travel times (surface-to-surface)      
   np = 500   
   call GETTX(zref, vref(1), nvzref, np, deltab, tttab, tfirstref)
   

   print *, 'Enter SS precursor stack file'
   read (*, '(a)') ssfile
   open (11, file=ssfile, status='old')
   
   read (11,13) xmin,xmax,delavg,delfix,qdepmin,      &
     qdepmax,stnmin,stnmax,wexp,iocean,icrustop,nboot,nseis
13 format (20x,8f6.1,f4.1,2i2,i4,i5)
   print *, 'average range from stack = ', delavg
   print *, 'fixed range for reduction = ', delfix   
   
   do i=-600,300
      read (11,21) time(i), amp(i), ampmin(i), ampmax(i)
21    format (2f10.5,30x,2f10.5)      
      if (i < -60) then
         amp(i) = amp(i)/10.               !remove stacker ten boost
         ampmin(i) = ampmin(i)/10.
         ampmax(i) = ampmax(i)/10.
      end if
   enddo   
   close (11)
   
   
   do i = 1, nvzref
      if (zref(i).eq.410.) i410=i
      if (zref(i).eq.560.) i560=i
      if (zref(i).eq.1205.5) ithou=i    !this was 1000, hence the name
      if (zref(i).eq.6371.) zref(i)=6370.   !don't go to center         
   enddo
   print *,'i410, i560, ithou = ',i410, i560, ithou
   
30    print *,'Enter pmin,pmax,np2 (np2.le.50)'
      print *,'Example: .096 .109 20'
      read *,pmin,pmax,np2
      zmax=2889.

! ray trace reference model to get SS times
      call RAYTRACE(zref,vref,nvzref,pmin,pmax,np2,zmax,             &
            p_iasp,x_iasp,t_iasp,tstar,x_iasp0,t_iasp0,tstar_tot)
      do ip=1,np2
         x_ss=2.*x_iasp0(ip)
         t_ss=2.*t_iasp0(ip)
         xthou=x_ss-x_iasp(ithou,ip)
         tthou=t_ss-t_iasp(ithou,ip)         
         print 34, p_iasp(ip),x_ss,t_ss,xthou,tthou, &
                   tstar_tot(ip),tstar(i410,ip),tstar(i560,ip)
34       format (f8.5,4f8.3,2x,3f6.3)
      enddo
      print *,'(1) continue  or (2) redo'
      read *,menu
      if (menu.eq.2) go to 30
      
      
   print *, 'Enter S.D. for Gaussian spread function (s)(e.g., 2.5 s)'
   read *, standev          

   print *, 'Enter display size scale factor (e.g., 1 for 15 inch)'
   read *, dfact
   dfactinv = 1./dfact   
   
   call EZXINIT(nint(1200.*dfact),nint(800.*dfact))   
   call EZXLOADFONT(i,'9x15')    
    
   xw1a = 50.*dfact
   xw1b = 600.*dfact
   yw1a = 60.*dfact
   yw1b = 750.*dfact
   
   xw2a = 670.*dfact
   xw2b = 1150.*dfact
   yw2a = 450.*dfact
   yw2b = 750.*dfact
   
   xw3a = xw2a
   xw3b = xw2b
   yw3a = 60.*dfact
   yw3b = 400.*dfact
   
   xp1a = vref(1) - 0.2              !starting values, user can change
   xp1b = vref(nvzref) + 0.4
   yp1a = zref(nvzref) + 50.
   yp1b = zref(1) -50.   
   
            
50 call GETTX(z, v(1), nvz, np, deltab, tttab, tfirst)
   rfirstmin = 999.
   rfirstmax = -999.
   do i = 11, 99
      rfirst(i) = tfirst(i) - tfirstref(i)
      if (abs(rfirst(i)) < 1000.) then
         if (rfirst(i) < rfirstmin) rfirstmin = rfirst(i)
         if (rfirst(i) > rfirstmax) rfirstmax = rfirst(i)
      endif
   enddo
!   print *, 'rfirstmin, rfirstmax = ', rfirstmin, rfirstmax
     
   
! ------------------------------------------display V(Z)
   call EZXWIND(xw1a, xw1b, yw1a, yw1b, xp1a, xp1b, yp1a, yp1b)
   call EZXNICE(xp1a, xp1b, xtic1, xtic2)
   call EZXNICE(yp1a, yp1b, ytic1, ytic2)
   call EZXAXES(xtic1, xtic2, 5, 'f7.1', ytic1, ytic2, 5, 'i8')
   call EZXLAX('Velocity (km/s)', 20, 'Depth (km)', 30)
   
   do i = 1, nvz
      if (i == 1) then
         call EZXMOVE(vref(i),zref(i))
      else
         call EZXDRAW(vref(i),zref(i))
      endif
   enddo   
   call EZXCOL(2)
   do i = 1, nvz
      if (i == 1) then
         call EZXMOVE(v(i),z(i))
      else
         call EZXDRAW(v(i),z(i))
      endif
      call EZXSYMB(1, 4)
   enddo
   call EZXCOL(1)   
   call EZXMOVE(xp1b, yp1b)
   call EZXTIC(-5, -5, 0)
   call EZXLORG(9)
   if (iw == 1) then
      call EZXLAB('P-wave velocity')
   else
      call EZXLAB('S-wave velocity')   
   endif
   
! ----------------------------------------------plot residuals
   xp2a = 9.5
   xp2b = 100.
   yp2a = rfirstmin - 1.7
   yp2b = rfirstmax + 1.7
   call EZXWIND(xw2a, xw2b, yw2a, yw2b, xp2a, xp2b, yp2a, yp2b)
   call EZXNICE(xp2a, xp2b, xtic1, xtic2)
   call EZXNICE(yp2a, yp2b, ytic1, ytic2)
   call EZXAXES(xtic1, xtic2, 5, 'i8', ytic1, ytic2, 5, 'f9.1')
   call EZXLAX('Distance (degrees)', 20, 'Residual (s)', 40)
   
      
!   do i = 11, 99
!      x = real(i) - 0.5
!      y = rfirst(i)
!      if (abs(y) < 500) then
!         if (i.eq.11) then
!            call EZXMOVE(x,y)
!         else
!            call EZXDRAW(x,y)
!         endif
!      endif
!   enddo
   
   do i = 1, np
      x = deltab(i)
      if (x > xp2a .and. x < xp2b) then
         ii = int(x - 0.5) + 1
         x1 = real(ii)-0.5
         x2 = x1 + 1.
         frac = (x-x1)/(x2-x1)
         tpred = tfirstref(ii)+frac*(tfirstref(ii+1)-tfirstref(ii))
         y = tttab(i) - tpred
         call EZXMOVE(x,y)
         call EZXSYMB(5, 2)
      endif
   enddo
               
   call EZXLORG(9)
   call EZXMOVE(xp2b, yp2b)
   call EZXTIC(-5, -5, 0)
   call EZXLAB('ak135 residuals')   
   
   
! now compute SS precursor synthetics (section adapted from old modvelden code)

! compute kluged densities that will always align with velocities
   velden = 1.39
   do i = 1, nvz
      den(i)=v(i)/velden
   enddo
   
!--------------------------------SET UP 10-km sampled V_SYN(Z_SYN) from V(Z)
      v_syn(1) = v(1)
      den_syn(1) = den(1)
      i = 2
      i10 = 1
      z10 = 10.
      n1200 = 0
    
60    if (z10.ge.1200..and.n1200.le.2) then
         z10 = 1200.
         n1200 = n1200+1
      else if (z10.ge.1210.) then
         go to 65
      end if
!      print *,i,i10,z10,z(i)
      i10 = i10 + 1
      if (z10.lt.z(i)-1.) then
         z_syn(i10) = z10
         call GETVEL(z10, z, v, nvz, v_syn(i10))
         call GETVEL(z10, z, den, nvz, den_syn(i10))
         z10 = z10 + 10.
         go to 60
      else
         z_syn(i10) = z(i)
         v_syn(i10) = v(i)
         den_syn(i10) = den(i)
         if (z(i).gt.1200.) go to 65
         i = i + 1
         if (z10 .lt. z(i)-1) z10 = z10 + 10.
         go to 60
      end if
65    n_syn = i10
      print *,'nvz,n_syn = ',nvz,n_syn
      print *,'z_syn(1),v_syn(1),den_syn(i) = ', &
               z_syn(1),v_syn(1),den_syn(i)
      print *,'z_syn(n_syn),v_syn(n_syn),den_syn(n_syn) = ', &
               z_syn(n_syn),v_syn(n_syn),den_syn(n_syn)
!      do 67 i=n_syn-5,n_syn
!67       print *,i,z_syn(i),v_syn(i),den_syn(i)


!   call EZXCOL(3)
!   call EZXWIND(xw1a, xw1b, yw1a, yw1b, xp1a, xp1b, yp1a, yp1b)
!   do i = 1, n_syn
!      if (i.eq.1) then
!         call EZXMOVE(v_syn(i),z_syn(i))
!      else
!         call EZXDRAW(v_syn(i),z_syn(i))
!      endif
!      call EZXSYMB(1,3)
!   enddo
!   call EZXCOL(1)

!
!-----------------------------------RAY TRACE UPPER 1200 km of model
      zmax = 1200.  
!      print *,'Entering raytrace: ',pmin,pmax,np2,zmax
      call RAYTRACE(z_syn,v_syn,n_syn,pmin,pmax,np2,zmax,     &
           p_syn,x_syn,t_syn,tstar,x_syn0,t_syn0,tstar_tot)
!      print *,'Ray trace results follow: '
!      do 68 ip=1,np2
!         print 69,p_iasp(ip),p_syn(ip),x_iasp(ithou,ip),  &
!                  t_iasp(ithou,ip),x_syn0(ip),t_syn0(ip)
!68    continue
!69    format (2f8.5,4f9.3)
      ttest0 = t_syn0(1)
      do 74 iz = 1, n_syn
         if (iz.gt.1) then              !check for Moho
            vv1=v_syn(iz-1)
            vv2=v_syn(iz)
            if (nplot.eq.1.and.vv1.lt.4.0.and.vv2.gt.4.0) then
               refl=(vv2**2-vv1**2)/(vv1**2+vv2**2)
               tmoho1=t_syn(iz,1)*60.
               tmoho2=t_syn(iz,np2)*60.
               print *,'Moho 2dt,refl = ',tmoho1,tmoho2,refl
            end if
         end if
         do ip=1,np2
            x_ss1=2.*x_iasp0(ip)
            t_ss1=2.*t_iasp0(ip)
            x_ss2=x_ss1-x_iasp(ithou,ip)+x_syn0(ip)
            t_ss2=t_ss1-t_iasp(ithou,ip)+t_syn0(ip)            
            x_comb(iz,ip)=x_ss2-x_syn(iz,ip)
            
!            print *, iz, ip, x_iasp0(ip), x_iasp(ithou,ip), x_syn0(ip), x_comb(iz,ip)
            
            t_comb(iz,ip)=t_ss2-t_syn(iz,ip)
         enddo
         do ip=2,np2
            if (x_comb(iz,ip).le.delfix.and.     &
                x_comb(iz,ip-1).gt.delfix) then
               frac=(delfix-x_comb(iz,ip-1))/   &
                    (x_comb(iz,ip)-x_comb(iz,ip-1))
               t_tab(iz)=t_comb(iz,ip-1)+               &
                         frac*(t_comb(iz,ip)-t_comb(iz,ip-1))
               p_tab(iz)=p_syn(ip-1)+frac*(p_syn(ip)-p_syn(ip-1))
               dx=x_comb(iz,ip)-x_comb(iz,ip-1)
               dp=p_syn(ip)-p_syn(ip-1)
               dxdp=abs(dx/dp)
               pavg=(p_syn(ip)+p_syn(ip-1))/2.
               u0=1./4.                           !assumed surface slowness
               s1=sqrt(pavg/(u0**2-pavg**2))
               s2=sqrt(1./dxdp)
               g_tab(iz)=s1*s2
!               print *,iz,dx,dp,s1,s2,g_tab(iz)
               go to 74
            end if
         enddo
         print *, '***WARNING: could not get delfix = ', delfix
         print *, 'iz,z_syn(iz) = ',iz,z_syn(iz)
         stop
74    continue
      print *,'finished raytracing upper 1200 km'

!75    print *,'Enter iz (0 to cont)'
!      read *,iz
!      if (iz.ne.0) then
!         s1=g_tab(iz)/g_tab(1)
!         print *,z_syn(iz),t_tab(iz),p_tab(iz),g_tab(iz),s1 
!         go to 75
!      end if

!-----------------------------------COMPUTE SYNTHETIC REFLECTIVITY TRACE
      do i=-600,300
         syn1(i)=0.
         syn2(i)=0.
      enddo

      t_ss0=t_tab(1)
      do 90 i=-600,-29
         tref=(float(i)-0.5)/60.  !pre-SS in minutes, offset 1/2 sample
         do iz=2,n_syn
            tref2=t_tab(iz)-t_ss0 
            if (tref2.le.tref) then
               tref1=t_tab(iz-1)-t_ss0
               if (tref2-tref1.ne.0.) then
                  frac=(tref-tref1)/(tref2-tref1)
                  zsyn(i)=z_syn(iz-1)+frac*(z_syn(iz)-z_syn(iz-1))
                  psyn(i)=p_tab(iz-1)+frac*(p_tab(iz)-p_tab(iz-1))
                  gsyn(i)=g_tab(iz-1)+frac*(g_tab(iz)-g_tab(iz-1))
                  go to 90
               else
                  print *,'***ERROR: tref1=tref2, iz = ',iz 
                  stop
               end if
            end if
          enddo
!         print *,'***depth not found in z_syn for i = ',i
         zsyn(i)=-999.
90    continue
!      do 95 i=-200,-150
!95       print *,i,zsyn(i)


!      write (13,*) nvz, n_syn

      do 100 i=-600,-30
         zz2=zsyn(i)
         zz1=zsyn(i+1)
         if (zz1.eq.-999..or.zz2.eq.-999.) go to 100
         call GETVEL(zz1,z,v,nvz,vv1)
         call GETVEL(zz1,z,den,nvz,den1)     !****now do density
         call GETVEL(zz2,z,v,nvz,vv2)
         call GETVEL(zz2,z,den,nvz,den2)
         if (vv1.eq.0..or.vv2.eq.0.) then
            print *,'Error in GETVEL. zz1,zz2,vv1,vv2 = ', zz1,zz2,vv1,vv2
            stop
         end if

         pavg=(psyn(i)+psyn(i+1))/2.    !this is flat earth p
         zavg=(zsyn(i)+zsyn(i+1))/2.
         ravg=r_earth-zavg
         psph=pavg*r_earth
         the1=asin(psph*vv1/ravg)
         the2=asin(psph*vv2/ravg)
!****changed to handle density jumps                  
         f1=cos(the2)*vv2*den2-cos(the1)*vv1*den1
         f2=cos(the2)*vv2*den2+cos(the1)*vv1*den1
         refl2=f1/f2
         geocor=gsyn(i)/g_tab(1)
         if (igeocor.eq.1) refl2=refl2*geocor
!          if (mod(i,50).eq.0) print *,i,geocor,igeocor,refl2
         syn1(i)=refl2
         if (zz1 < 660. .and. zz2 > 660.) then
!            write (13,*) zz1,zz2,vv1,vv2,den1,den2,refl2,geocor,pavg, the1, the2
         endif
!         if (zz1.gt.1000.) then
!         print *,zz1,zz2,vv1,vv2,pavg,ravg,psph,the1,the2,refl2,refl1
!         end if
100   continue          
         

!-----------------------------------------CONVOLVE with REFERENCE PULSE
      iref1 = -30
      iref2 = 60
      do 120 i=-600,300
         if (syn1(i).eq.0.) go to 120
         do 110 ioff=iref1,iref2
            ii=i+ioff
            if (ii.lt.-600.or.ii.gt.300) go to 110
            syn2(ii)=syn2(ii)+syn1(i)*amp(ioff)
110      continue
120   continue

! -----------------------------------------CONVOLVE WITH GAUSS FUNC
      if (standev.ne.0.) then
         call GAUSSOP(syn2,syn3,901,1.,standev)
         do i=-600,300
            syn2(i)=syn3(i)
         enddo
      end if
            
   
! plot SS precursors from stack
   xp3a = -360.
   xp3b = -60.
   yp3a = -0.04
   yp3b = 0.06
   call EZXWIND(xw3a, xw3b, yw3a, yw3b, xp3a, xp3b, yp3a, yp3b)
   call EZXNICE(xp3a, xp3b, xtic1, xtic2)
   call EZXNICE(yp3a, yp3b, ytic1, ytic2)
   call EZXAXES(xtic1, xtic2, 5, 'i8', ytic1, ytic2, 5, 'f9.2')
   call EZXLAX('Time (s)', 20, 'Relative amplitude', 40)
      
   do i = -360, -60
      x = time(i)*60.
      y = amp(i)
      if (i == -360) then
         call EZXMOVE(x,y)
      else
         call EZXDRAW(x,y)
      endif
   enddo
   
! plot spike functions
!   call EZXCOL(3)
!   do i = -360, -60
!      x = time(i)*60.
!      y = syn1(i)
!      if (i == -360) then
!         call EZXMOVE(x,y)
!      else
!         call EZXDRAW(x,y)
!      endif
!   enddo   

! plot synthetic SS precursors
   call EZXCOL(2)
   do i = -360, -60
      x = time(i)*60.
      y = syn2(i)
      if (i == -360) then
         call EZXMOVE(x,y)
      else
         call EZXDRAW(x,y)
      endif
   enddo
   call EZXCOL(1)
   call EZXLORG(9) 
   call EZXMOVE(xp3b, yp3b)
   call EZXTIC(-5, -5, 0)
   call EZXLAB('SS precursors')       

      
150 call EZXMESS('(1) change V(Z) point, (2) remove V(Z) point, (3) add V(Z) point, &
   (4) rescale V(Z), (5) output model, (6) quit')

   call EZXWIND(xw1a, xw1b, yw1a, yw1b, xp1a, xp1b, yp1a, yp1b)   
   call EZXDIG(x, y, iscr)
!   print *, 'x,y,iscr = ',x,y,iscr
   
   
   if (iscr == 1 .or. iscr == 2) then
      distmin = 9999.
      do i = 1, nvz
         dz = abs( (y-z(i))*(yw1b-yw1a)/(yp1b-yp1a) )
         dx = abs( (x-v(i))*(xw1b-xw1a)/(xp1b-xp1a) )
         dist = sqrt(dx**2 + dz**2)
         if (dist < distmin) then
            distmin = dist 
            ii = i
         endif
      enddo
      if (iscr == 1) then
         z(ii) = y
         v(ii) = x
         do i = ii+1, nvz
            if (z(i) < z(ii)) z(i) = z(ii)
         enddo
         do i = ii-1, 1, -1
            if (z(i) > z(ii)) z(i) = z(ii)
         enddo
         
      else
         do i = ii, nvz - 1
            z(i) = z(i+1)
            v(i) = v(i+1)
         enddo
         nvz = nvz -1
      endif
      call EZXCLR
      go to 50         
   else if (iscr == 3) then
      do i = 2, nvz
         if (z(i) > y) then
            ii = i 
            exit
         endif
      enddo
      do i = nvz, ii, -1
         z(i+1) = z(i)
         v(i+1) = v(i)
      enddo
      z(ii) = y
      v(ii) = x
      nvz = nvz+1
      call EZXCLR
      go to 50
   else if (iscr == 4) then
      call EZXMESS('(1) enter lower right point (3) zoom out')
      call EZXDIG(xx, yy, iscr)
      if (iscr == 3) then
         xp1a = vref(1) - 0.2
         xp1b = vref(nvzref) + 0.4
         yp1a = zref(nvzref) + 50.
         yp1b = zref(1) -50.
         call EZXCLR
         go to 50         
      else
         xp1a = min(x, xx)
         xp1b = max(x, xx)
         yp1a = max(y, yy)
         yp1b = min(y, yy)
         call EZXCLR
         go to 50      
      endif
   else if (iscr == 5) then
      open (14, file = 'out.vzmod')
      do i = 1, nvz
         write (14, 177) z(i), v(i), den(i)
177      format (3f10.4)
      enddo
      close (14)
      print *,'Model written to out.vzmod, npts = ', nvz 
      go to 150
   endif 

   call EZXQUIT
!   close (13)


end program playmod




subroutine GETTX(z, v, nvz, np, deltab, tttab, tfirst)
   implicit none
   integer, parameter :: nvz0 = 1000, np0 = 5000
   integer :: nvz, np, i, ip, imth, irtr
   real :: erad, pi, ecircum, kmdeg, degrad, pmin, pmax, dp, p, x, t, h, dx, dt, &
           tmin, frac
   real, dimension(nvz0) :: z, zf, v, vf, slow
   real, dimension(np0) :: deltab, tttab
   real, dimension(100) :: tfirst
   
   erad = 6371.
   pi = 3.1415927
   ecircum = 2.*pi*erad
   kmdeg = ecircum/360.
   degrad = 180./pi     

   do i = 1, nvz 
      call FLATTEN(z(i), v(i), zf(i), vf(i))
      slow(i) = 1./vf(i)
   enddo

   pmin = slow(nvz)   
   do i = 1, nvz
      if (slow(nvz) < pmin) pmin = slow(nvz)
   enddo
      
   pmin = pmin + 0.0001
   pmax = slow(1) - 0.0001
   
   dp = (pmax - pmin)/real(np - 1)
   
   open (12, file = 'tx.out')
   
   do ip = 1, np
      p = pmin + dp*real(ip - 1)
      x = 0.
      t = 0.
      do i = 1, nvz - 1
         h = zf(i+1) - zf(i)
         imth = 3
         call LAYERTRACE(p, h, slow(i), slow(i+1), imth, dx, dt, irtr)
         x = x + dx
         t = t + dt         
         if (irtr == 0 .or. irtr == 2) exit
      enddo
      x = 2.*x
      t = 2.*t
      deltab(ip) = x/kmdeg
      tttab(ip) = t
      write (12, *) p, deltab(ip), tttab(ip), tttab(ip) - 15.*deltab(ip)
   enddo
   close (12)
   
   do i = 1, 10
      tfirst(i) = 0.
   enddo

   open (12, file ='tx2.out')   
   do i = 11, 100
      x = real(i) - 0.5
      tmin = 99999.
      do ip = 1, np-1
         if ((deltab(ip) <= x .and. deltab(ip+1) >= x) .or. &
             (deltab(ip+1) <= x .and. deltab(ip) >= x) ) then
            frac = (x - deltab(ip))/(deltab(ip+1)-deltab(ip))
            t = tttab(ip) + (tttab(ip+1)-tttab(ip))*frac
            if (t < tmin) tmin = t
         endif
      enddo
      tfirst(i) = tmin
      write (12,*) x, tmin
   enddo
   close (12)
   
   return
end subroutine GETTX





! FLATTEN calculates flat Earth transformation
subroutine FLATTEN(z_s, vel_s, z_f, vel_f)
   implicit none
   real :: z_s, vel_s, z_f, vel_f, erad, r   
   erad = 6371.
   r = erad - z_s
   z_f = -erad*alog(r/erad)
   vel_f = vel_s*(erad/r)
   return
end subroutine FLATTEN



