      subroutine delth23(del,dep,depq,th,t,dtdr,dtddel,kps)
      common /cvel/al(2,2),be(2,2),r,rm,ntype_velmodel
cc      write(*,*) 'dep,depq=',dep,depq,((al(k,kq),be(k,kq),k=1,2),kq=1,2)
      r1=r-dep
      r2=r-depq
      if(r2.lt.r1) then
      r2=r-dep
      r1=r-depq 
      endif
cc      write(*,*) 'r1,r2=',r1,r2,rm,r
      va1=al(1,kps)
      va2=al(2,kps)
      be1=be(1,kps)
      be2=be(2,kps)
      if(r1.gt.rm) then
      v1=va1*(r1/r)**be1
      else
      v1=va2*(r1/rm)**be2
      endif
      if(r2.gt.rm) then
      v2=va1*(r2/r)**be1
      else
      v2=va2*(r2/rm)**be2
      endif
      vd=va1*(rm/r)**be1
      q1=r1/v1
      qm=rm/va2
      qd=rm/vd
      qs=r2/v2
      a1=1.0-be1
      a2=1.0-be2
      delt=del/r
cc      write(*,*) 'qs=',qs
      call hadel(qm,a1,qs,qd,crdel)
      if(r1.le.rm.and.r2.le.rm) then
      qs1=q1-qs
      s=sin(a2*delt/2.0)
      t=sqrt2(qs1*qs1+4.0*qs*q1*s*s)/a2
      th=atan3(sin(a2*delt),q1/qs-cos(a2*delt))
      return
      endif
      if(r1.lt.rm) go to 1
      qs1=q1-qs
      qsd=qs-qd
      s=sin(a1*delt/2.0)
      tg=sqrt2(qs1*qs1+4.0*qs*q1*s*s)/a1
      call hadel(qm,a1,q1,qd,d1)
      d2=delt-d1-crdel
cc      write(*,*) 'd2=',d2,' tg=',tg
      if(d2.lt.0.0) go to 2
ccc	travel time of Pn	***********************************
      do 3 k=1,30
      dth=d2/2.0*a2
      hp=cos(dth)*qm
      call hadel(hp,a1,q1,qd,d11)
      call hadel(hp,a1,qs,qd,d12)
      dsa=abs(delt-d2-d11-d12)*r
      d2=delt-d11-d12
      if(dsa.lt.0.1) go to 4
    3 continue
      go to 999
    4 continue
      s1=sin(d11*a1/2.0)
      s2=sin(d12*a1/2.0)
      tq=(sqrt2((q1-qd)*(q1-qd)+4.0*q1*qd*s1*s1)+sqrt2(qsd*qsd+
     1 4.0*qs*qd*s2*s2))/a1+2.0*qm*sin(d2*a2/2.0)/a2
cc      write(*,*) 'tq=',tq,tg
      if(tq.gt.tg) go to 2
    7 t=tq
      s=hp/q1
      c=sqrt2(1.0-s*s)
      th=atan3(s,c)

c     write(*,*) 'tq=',tq

      go to 998
    2 th=atan3(sin(a1*delt),q1/qs-cos(a1*delt))
      sg=del/(abs(r2-r1)+0.01)
      if(sg.lt.0.01.and.th.lt.1.5) th=3.14159265
      sg=sin(th)-qd/q1
cc      write(*,*) 'after do 2 th=',th
      if(th.lt.1.57.and.sg.lt.0.0) then
	if(d2.gt.0.0) go to 7		!r1>radius
	endif
      delmax=3.14/a1
c      write(*,*) 'delmax=',delmax,delt,'tg=',tg
      if(delt.gt.delmax) go to 7
      t=tg
      dtddel=sin(th)/v1
      dtdr=cos(th)/v1
c      write(*,*) 't=',t
      if(depq.gt.dep) then
         thr=atan3(sin(a1*delt),qs/q1-cos(a1*delt))
      th=thr
      dtddel=sin(th)/v2
      dtdr=cos(th)/v2
         endif
         return
    1 continue
      qds2=(qd-qs)*(qd-qs)
      qm2=(q1-qm)*(q1-qm)
      d1=delt-crdel
      if(d1.gt.0.0) go to 11
	if(r.eq.r1) go to 11
      d1=delt*(rm-r1)/(r-r1)
      dh=rm-r1
      if(dh.gt.0.2) go to 11
      d2=delt
      thd=atan3(sin(a1*delt),qd/qs-cos(a1*delt))
      s=sin(d2*a1/2.0)
      t=sqrt2((qs-qd)*(qs-qd)+4.0*qs*qd*s*s)/a1
      s=v1/vd*sin(thd)*rm/r1
      c=-sqrt2(abs(1.0-s*s))
      th=atan3(s,c)
      t=t-dh*c/v1
      go to 998
   11 continue
      do 5 k=1,30
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
           hp=sin(th)*q1
      hp2=hp*hp
      call hadel(hp,a1,qs,qd,d2)
      dsa=abs(delt-d1-d2)*r
      sc=sin(a1*d2/2.0)
      sm=sin(a2*d1/2.0)
      tc=sqrt2(qds2+4.0*qd*qs*sc*sc)/a1
      tm=sqrt2(qm2+4.0*q1*qm*sm*sm)/a2
      if(dsa.lt.0.1) go to 6
	qe=(qs*qs-hp2)*(qd*qd-hp2)
	if(qe.gt.-0.00001.and.qe.le.0.00001) qe=0.00001
	if(tm.gt.-0.00001.and.tm.le.0.00001) tm=0.00001	
      ddd1=1.0+sqrt2((q1*q1-hp2)*(qm*qm-hp2)/qe)*tc/tm
      db=d1
	if(ddd1.gt.-0.00001.and.ddd1.le.0.00001) ddd1=0.00001		
      d1=d1+(delt-d1-d2)/ddd1
      if(d1.lt.0.000001) d1=db/2.0
      if(d1.gt.delt) d1=delt
    5 continue
      go to 999
    6 t=tc+tm
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
  998 continue
      if(th.lt.0.01) th=3.14159265-th
      dtddel=sin(th)/v1
      dtdr=cos(th)/v1
      return
  999 t=9999.0
      return
      end
        subroutine delhh(heitk,dep,del,th,tc,dtdr,dtddel,kps)
      parameter(vpsf=2.5,vssf=1.0)
	include 'common_stncod.inc'    
cc        ntype_velmodel=1:Tohoku,2:JMA vel model,3:Layered       
	if(ntype_velmodel.eq.1) then
       depq=-heitk
cc       call delth3(del,dep,depq,th,tc,kps)
cc       r1=radius-dep
cc       call  delth2(r1,del,th,tc,dtdr,dtddel,kps)
cc       tc=tc+heitk/al(1,kps)
       call delth23(del,dep,depq,th,tc,dtdr,dtddel,kps)
       go to 88
       endif
	if(ntype_velmodel.eq.2) then
ccc   JMA velocity structure
       call delth_jma2(heitk,dep,del,th,tc,dtdr,dtddel,kps)
       go to 88
       endif
       if(ntype_velmodel.gt.6.or.ntype_velmodel.le.0) then
          write(*,*) ntype_velmodel
       write(*,*) 'velocity type in read_vel_model must be 1,2,3!!!'
       stop
       endif
       call delhh_lay(heitk,dep,del,th,tc,dtdr,dtddel,kps)
       return
 88    continue
c      if(kps.eq.1) tc=tc+heitk/vpsf
c      if(kps.eq.2) tc=tc+heitk/vssf
       return
       end
       subroutine del_trv_hei(tp,dep,del,hei)
      del=0.0
      kps=1
       call delhh(hei,dep,del,th,tc,dtdr,dtddel,kps)
cc       write(*,*) 'hei,dep,del,th,tc,tp=',hei,dep,del,th,tc,tp
       if(tc.gt.tp) then
          del=-0.01
          return
          endif
          del1=0.0
          dd=20.0
          do kkk=1,10
             if(dd.lt.0.01) go to 5
          do k=1,1000
             del=float(k)*dd+del1
                 call delhh(hei,dep,del,th,tc,dtdr,dtddel,kps)
                 sa=tp-tc
cc                 write(*,*) dd,hei,dep,del,tc,' sa=',sa
                 if(tc.gt.tp) go to 1
                 enddo
 1               continue
                 del1=del-dd
                 dd=dd/5.0
            enddo
 5          continue
           
            return
            end
        subroutine del_trv(tp,dep,del)
	include 'common_stncod.inc'    
ccc   calculate approximate epicentral distance by putting Tp and depth,
ccc   input tp:P wave travel time, dep;depth
ccc   output del:epicentral distance 
        if(al(1,1).le.0.1) then
           al(1,1)=5.4
           al(2,1)=7.75
           be(1,1)=-59.0
           be(2,1)=-2.3
           endif
           kps=1
           delsa=100.0
           delk=0.0
           do m=1,4
           do k=1,30
              del=delk+float(k)*delsa
           call delth3(del,dep,0.0,th,tpq,kps)
cc           write(*,*) del,dep,tp,tpq
           if(tpq.gt.tp) go to 1
           enddo
 1         continue
           delk=del-delsa
           delsa=delsa/10.0
           enddo
           call delth3(del,dep,dep,th,tpq,kps)
           sa=tp-tpq
cc           write(*,*)  'del_trv',tp,dep,del,sa
           return
           end    
  
      subroutine delth3(del,dep,depq,th,t,kps)
      common /cvel/al(2,2),be(2,2),r,rm,ntype_velmodel
cc      write(*,*) 'dep,depq=',dep,depq,((al(k,kq),be(k,kq),k=1,2),kq=1,2)
      r1=r-dep
      r2=r-depq
      if(r2.lt.r1) then
      r2=r-dep
      r1=r-depq 
      endif
cc      write(*,*) 'r1,r2=',r1,r2,rm,r
      va1=al(1,kps)
      va2=al(2,kps)
      be1=be(1,kps)
      be2=be(2,kps)
      if(r1.gt.rm) then
      v1=va1*(r1/r)**be1
      else
      v1=va2*(r1/rm)**be2
      endif
      if(r2.gt.rm) then
      v2=va1*(r2/r)**be1
      else
      v2=va2*(r2/rm)**be2
      endif
      vd=va1*(rm/r)**be1
      q1=r1/v1
      qm=rm/va2
      qd=rm/vd
      qs=r2/v2
      a1=1.0-be1
      a2=1.0-be2
      delt=del/r
cc      write(*,*) 'qs=',qs
      call hadel(qm,a1,qs,qd,crdel)
      if(r1.le.rm.and.r2.le.rm) then
         qs1=q1-qs
      s=sin(a2*delt/2.0)
      t=sqrt2(qs1*qs1+4.0*qs*q1*s*s)/a2
      th=atan3(sin(a2*delt),q1/qs-cos(a2*delt))
      return
      endif
      if(r1.lt.rm) go to 1
      qs1=q1-qs
      qsd=qs-qd
      s=sin(a1*delt/2.0)
      tg=sqrt2(qs1*qs1+4.0*qs*q1*s*s)/a1
      call hadel(qm,a1,q1,qd,d1)
      d2=delt-d1-crdel
      if(d2.lt.0.0) go to 2
ccc	travel time of Pn	***********************************
      do 3 k=1,30
      dth=d2/2.0*a2
      hp=cos(dth)*qm
      call hadel(hp,a1,q1,qd,d11)
      call hadel(hp,a1,qs,qd,d12)
      dsa=abs(delt-d2-d11-d12)*r
      d2=delt-d11-d12
      if(dsa.lt.0.1) go to 4
    3 continue
      go to 999
    4 continue
      s1=sin(d11*a1/2.0)
      s2=sin(d12*a1/2.0)
      tq=(sqrt2((q1-qd)*(q1-qd)+4.0*q1*qd*s1*s1)+sqrt2(qsd*qsd+
     1 4.0*qs*qd*s2*s2))/a1+2.0*qm*sin(d2*a2/2.0)/a2
cc      write(*,*) 'tq=',tq,tg
      if(tq.gt.tg) go to 2
    7 t=tq
      s=hp/q1
      c=sqrt2(1.0-s*s)
      th=atan3(s,c)
      go to 998
    2 th=atan3(sin(a1*delt),q1/qs-cos(a1*delt))
      sg=del/(abs(r2-r1)+0.01)
      if(sg.lt.0.01.and.th.lt.1.5) th=3.14159265
      sg=sin(th)-qd/q1
      if(th.lt.1.57.and.sg.lt.0.0) then
	if(d2.gt.0.0) go to 7		!r1>radius
	endif
      delmax=3.14/a1
      if(delt.gt.delmax) go to 7
      t=tg
      go to 998
    1 continue
      qds2=(qd-qs)*(qd-qs)
      qm2=(q1-qm)*(q1-qm)
      d1=delt-crdel
      if(d1.gt.0.0) go to 11
	if(r.eq.r1) go to 11
      d1=delt*(rm-r1)/(r-r1)
      dh=rm-r1
      if(dh.gt.0.2) go to 11
      d2=delt
      thd=atan3(sin(a1*delt),qd/qs-cos(a1*delt))
      s=sin(d2*a1/2.0)
      t=sqrt2((qs-qd)*(qs-qd)+4.0*qs*qd*s*s)/a1
      s=v1/vd*sin(thd)*rm/r1
      c=-sqrt2(abs(1.0-s*s))
      th=atan3(s,c)
      t=t-dh*c/v1
      go to 998
   11 continue
      do 5 k=1,30
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
           hp=sin(th)*q1
      hp2=hp*hp
      call hadel(hp,a1,qs,qd,d2)
      dsa=abs(delt-d1-d2)*r
      sc=sin(a1*d2/2.0)
      sm=sin(a2*d1/2.0)
      tc=sqrt2(qds2+4.0*qd*qs*sc*sc)/a1
      tm=sqrt2(qm2+4.0*q1*qm*sm*sm)/a2
      if(dsa.lt.0.1) go to 6
	qe=(qs*qs-hp2)*(qd*qd-hp2)
	if(qe.gt.-0.00001.and.qe.le.0.00001) qe=0.00001
	if(tm.gt.-0.00001.and.tm.le.0.00001) tm=0.00001	
      ddd1=1.0+sqrt2((q1*q1-hp2)*(qm*qm-hp2)/qe)*tc/tm
      db=d1
	if(ddd1.gt.-0.00001.and.ddd1.le.0.00001) ddd1=0.00001		
      d1=d1+(delt-d1-d2)/ddd1
      if(d1.lt.0.000001) d1=db/2.0
      if(d1.gt.delt) d1=delt
    5 continue
      go to 999
    6 t=tc+tm
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
  998 continue
      if(th.lt.0.01) th=3.14159265-th
      dtddel=sin(th)/v1
      dtdr=cos(th)/v1
      return
  999 t=9999.0
      return
      end
         subroutine trv_head(k1,k2,heit,dep,del,tn,thn,kps)
      common  /c_vel_hh/alnso(50,2),hdepa(50),nso
cc         write(*,*) 'head wave ',heit,dep,del
         tmin=999999.9
         do kk=k2,nso-1
            tz=del/alnso(kk+1,kps)
            tz=0.0
            DHH=0.0
            dd=0.0
           v2=alnso(kk+1,kps)      
            do kd=k2,kk
            if(kd.eq.k2) then
               h=hdepa(kd+1)-dep
               else
                  h=hdepa(kd+1)-hdepa(kd)
                  endif
              v1=alnso(kd,kps)
              tz=tz+sqrt(v2*v2-v1*v1)/v2/v1*h
              DHH=DHH+H
              s=v1/v2
              c=sqrt(1.0-s*s)
              dd=dd+h*s/c

cc             write(*,*) 'down ',kk,kd,h,tz,V1,V2
              enddo
ccc   going up
              do ku=k1,kk
                 if(ku.eq.k1) then
                    h=hdepa(ku+1)-heit
                    else
                       h=hdepa(ku+1)-hdepa(ku)
                       endif
              v1=alnso(ku,kps)
              tz=tz+sqrt(v2*v2-v1*v1)/v2/v1*h
              DHH=DHH+H
              s=v1/v2
              c=sqrt(1.0-s*s)
              dd=dd+h*s/c
cc              write(*,*) 'up   ',kk,ku,h,tz,V1,V2
  
              enddo
              tz1=tz 
           tz=tz+del/alnso(kk+1,kps)
           if(dd.gt.del) tz=999.9
cc           WRITE(*,108) k1,k2,KK,dep,DHH,alnso(kk+1,kps),TZ1,TZ,tmin
 108       format(3i3,6f11.3)
cc           write(*,*) (hdepa(j),j=1,nso)

             if(tz.lt.tmin) then
                tmin=tz
                smin=v1/v2
                endif
                enddo
                tn=tmin
               c=sqrt(1.0-smin*smin)
                thn=atan2(s,c)
                return
                end

      subroutine delhh_lay(heitk,depk,del,th,t,dtdr,dtddel,kps)
      common  /c_vel_hh/alnso(50,2),hdepa(50),nso
         pideg=3.14159265/180.0
          pi2=3.14158265/2.0 
          heit=-heitk    !SIGN change heit become positeve if deeper
          dep=depk
            t=0.0
          th=0.0
cc          write(*,*) 'nso=',nso
ccc   search velocity at focus. 
      do k=2,nso
         k2=k-1
         if(dep.lt.hdepa(k)) go to 1
         enddo
 1       continue
             do k=2,nso
         k1=k-1
         if(heit.lt.hdepa(k)) go to 2
         enddo  
 2       continue
         kstndep=0
         if(k1.gt.k2) then
            kk=k1
            k1=k2
            k2=kk
            dep=-heitk
            heit=depk
            kstndep=1
            endif
cc            write(*,*) 'k1,k2,dep,heit=',k1,k2,dep,heit
            if(k2.le.0) k2=1
        h=dep-hdepa(k2)
        kk2=k2
        if(h.le.0.01) kk2=k2-1
        if(kk2.le.1) kk2=2
ccc     travel time of head wave.
        call trv_head(k1,kk2,heit,dep,del,tpn,thn,kps)
cc        write(*,*) 'tpn=',tpn  
ccc   direct wave *******************************************
       
       z=dep-heit
        dis=sqrt(del*del+z*z)    
         v1=alnso(k2,kps)     
        if(k1.eq.k2) then
           t=dis/v1
           th=3.14159
           if(dis.gt.0.01) th=atan2(del,z)
           go to 4
           endif
ccc   k2.ge.2

         th1=atan2(del,z)
         thw=th1/pideg
        call trv_up(th1,k1,kk2,heit,dep,dd1,t,kps)
cc        write(*,*) 'begin th d,t=',kk2,thw,dd1,t
        if(th1.le.0.01) then
           th=3.14159-th1
           go to 888
           endif
        pi2s=pi2*0.95
           do m=1,30
              th2=th1+0.05*float(m)
cc              write(*,*) 'th1,th2=',th1,th2
              if(th2.gt.pi2s) th2=pi2s
        call trv_up(th2,k1,kk2,heit,dep,dd2,t2,kps)
cc        write(*,*) ' del,th2,dd2,t=',del,th2,dd2,t2
        if(dd2.gt.del) go to 23
        if(th2.ge.pi2s) go to 24
        enddo
 24     continue
ccc   long distance.
        dep2=hdepa(kk2-1)
        vu=alnso(kk2-1,kps) 
        s=vu/v1
        c=sqrt(1.0-s*s)
        thh=atan2(s,c)
        call trv_up(thh,k1,kk2-1,heit,dep2,dd2,t2,kps)
        th=pi2
        h=dep-hdepa(kk2-1)
        t=t2+sqrt(h*h+(del-dd2)*(del-dd2))/v1
        dsa=del-dd2
cc        write(*,*) 'long dis ',thh,dep2,heit,dd2,t2,t
cc        write(*,*) 'dep,h,dsa=',dep,h,dsa
cc        write(*,*) 'k1,k2,v1,v2=',k1,k2,v1,vu
        if(dsa.le.0.01) go to 4


        do mm=1,10
        th=atan2(dsa,h)
        s=vu/v1*sin(th)
        c=sqrt(1.0-s*s)
        thh=atan2(s,c)
        call trv_up(thh,k1,kk2-1,heit,dep2,dd2,t2,kps)
        dd=dsa+dd2
        er=abs(del-dd)
        th=3.14159265-thh
cc        write(*,*) del,dd,dsa,dd2,thh,'  er=',er
        t=t2+sqrt(h*h+(del-dd2)*(del-dd2))/v1
        if(er.le.0.001) go to 4
        dsa=del-dd2
        enddo
        go to 4
 23     continue
cc         write(*,*) 'after 23 th1,d1,2=',th1,th2,dd1,dd2,del



ccc   wave going up.
         do k=1,100
            th=th1+(del-dd1)/(dd2-dd1)*(th2-th1)*1.0
        call trv_up(th,k1,kk2,heit,dep,dd,t,kps)
        dsa=del-dd
cc        write(*,*) k,' th,d,t=',th,dd,t,' dsa=',dsa
 100    format(i4,8f8.4)
        dsa=del-dd
        if(dsa.gt.-0.01.and.dsa.le.0.01) go to 4
        if(dsa.gt.0) then
          dd1=dd
          th1=th
          else
             dd2=dd
             th2=th
             endif
             enddo
 4           continue
cc             write(*,*) 'direct ',t,th

             th=3.14159265-th
cc             write(*,*) 't,tpn=',t,tpn , ' th,thn=',th,thn
             if(kstndep.eq.1) th=3.14159265-th
             if(t.gt.tpn) then
                t=tpn
                th=thn
                endif
 888            continue

      dtddel=sin(th)/v1
      dtdr=cos(th)/v1

        return
        end 






         subroutine trv_up(thg,k1,k2,heit,dep,dd,t,kps)
      common  /c_vel_hh/alnso(50,2),hdepa(50),nso
         dd=0.0
         tt=0.0
         v1=alnso(k2,kps)
         s1=sin(thg)
                  thw=thg*180.0/3.14159
cc         write(*,*) 'v1,k1,k2=',v1,k1,k2
cc         write(*,*) 'thw=',thw, 'dep=',dep
           if(k1.le.0) k1=1
         do k=k1,k2
         v=alnso(k,kps)
         if(k.eq.k2) then
         th=thg
         h=dep-hdepa(k)
         else
         v=alnso(k,kps)
         s=v*s1/v1
cc         write(*,*) 's=',s
         c=sqrt(abs(1.0-s*s))
         th=atan2(s,c)
         h=hdepa(k+1)-hdepa(k)
         endif
         if(k.eq.k1) h=hdepa(k+1)-heit
         d=h*tan(th)
         dd=dd+d
         dis=sqrt(d*d+h*h)
         tt=tt+dis/v
         ts=dis/v
         thw=th*180.0/3.14159
cc         write(*,*) k,hdepa(k),v,' h,d,ts=',h,d,ts,tt
         enddo
         t=tt
cc         write(*,*) tt,dd
cc         read(*,*) iok
         return
         end


      subroutine delth_jma2(heit,dep,del,th,tc,dtdr,dtddel,kps)
      common /c_trav_jma/idelg(236),idepg(106),travp(106,236),
     1 travs(106,236) 
      data ibgn /99/
         if(ibgn.eq.99) then
         ibgn=0
         call read_jma_trv
               endif
                   if(del.lt.2.0.and.dep.lt.2.0) then
                     v=4.8
                     if(kps.eq.2) v=2.84
                     tc=sqrt(del*del+dep*dep)/v
                     if(del.lt.0.01) del=0.01
                     th=atan2(del,-dep)
                     dtddel=v*sin(th)
                     dtdr=v*cos(th)
                     return
                     endif 
                 
                     if(del.gt.1999.0.or.dep.gt.699.0) then
                        tc=del/10.0+dep/10.0
                        return
                        endif

                  call delth_jma(del,dep,tc,kps)

cc                  write(*,*)  'del,dep,tc,kps',del,dep,tc,kps

                  call delth_jma(del+1.0,dep,t1,kps)
                  call delth_jma(del,dep+1.0,t2,kps)
                  dtddel=t1-tc
                  dtdr=-(t2-tc)
cc               write(*,*) kps,del,dep,tc,dtddel,dtdr
               return
               end


        subroutine cal_idel_idep(del,dep,idel,idep,ddel,ddep)
        if(del.lt.50.0) then
         idel=ifix(del/2.0)+1
         ddel=2.0
         go to 1
         endif
         if(del.lt.200.0) then
            idel=ifix((del-50.0)/5.0)+26
            ddel=5.0
            go to 1
            endif
            idel=ifix((del-200.0)/10.0)+56
            ddel=10.0
 1          continue

cc            write(*,*) 'dep=',dep
      if(dep.lt.50.0) then
         idep=ifix(dep/2.0)+1
         ddep=2.0
         go to 2
         endif
         if(dep.lt.200.0) then
            idep=ifix((dep-50.0)/5.0)+26
            ddep=5.0
            go to 2
            endif
            idep=ifix((dep-200.0)/10.0)+56
            ddep=10.0
 2          continue
            if(idep.ge.106) idep=106
            if(idel.gt.236) idel=236
            if(idel.le.0) idel=1
            if(idep.le.0) idep=1
            return
            end
      subroutine delth_jma(del,dep,t,kps)
      common /c_trav_jma/idelg(236),idepg(106),travp(106,236),
     1 travs(106,236) 
      
      call cal_idel_idep(del,dep,idel,idep,ddel,ddep)

c            WRITE(*,*) 'delth_jma, del,dep=',kps,del,dep,idel,idep,
c     1 'delg=',idelg(idel),idepg(idep)
            depsa=abs(dep-float(idepg(idep)))
            if(depsa.gt.ddep) then
cc               write(*,*) 'dep,idep,sa=',dep,idepg(idep),depsa
cc               read(*,*) iok
               endif
            if(kps.eq.1) then

               if(idel.le.0)  then
                  write(*,*) 'delthjma error ',del,dep,kps,idel
                  
                  stop
                  endif
 


        tra1=travp(idep,idel)
        tra2=travp(idep,idel+1)
        tra3=travp(idep+1,idel)
        tra4=travp(idep+1,idel+1)
        else
        tra1=travs(idep,idel)
        tra2=travs(idep,idel+1)
        tra3=travs(idep+1,idel)
        tra4=travs(idep+1,idel+1)
        
        endif
        q1=(del-float(idelg(idel)))/ddel
        q2=(dep-float(idepg(idep)))/ddep
        t1=tra1+q1*(tra2-tra1)
        t2=tra3+q1*(tra4-tra3)
        t=t1+(q2*(t2-t1))
cc       write(*,'(11f8.2)') tra1,tra2,tra3,tra4,q1,q2,ddel,ddep,t1,t2,t
cc        read(*,*) iok
        return
        end



      subroutine read_jma_trv
      common /c_trav_jma/idelg(236),idepg(106),travp(106,236),
     1 travs(106,236) 
        open(7,file='../etc/jma_trav_del_dep.dat')
      ndel=236
      ndep=106
      do kdep=1,ndep
      do kdel=1,ndel
         read(7,100) travp(kdep,kdel),travs(kdep,kdel),
     1 idepg(kdep),idelg(kdel)
 100     format(1x,f9.3,2x,f9.3,i4,i7)
c         if(kdel.eq.1) write(*,100) travp(kdep,kdel),travs(kdep,kdel),
c     1 idepg(kdep),kdep
         enddo
         enddo
         close(7)
c         write(*,'(10i5)') (idelg(j),j=1,ndel)
c         write(*,'(10i7)') (idepg(j),j=1,ndep)
         return
         end
