       SUBROUTINE ellip()
C==========================================================================
C                                                                         
C    Ellipticity correction for any given phase using
C    Dziewonski & Gilbert representation                                                               C                                                                         
C      The ellipticity corrections are found by linear interpolation       
C    in terms of values calculated for the ak135 model for a wide 
C    range of phases to match the output of the iasp software 
C
C     first call:  ellref(ecolat) 
C                        - to set up source dependent constants
C     2nd call  :  ellcor(phase,edist,depth,ecolat,azim,tcor,abrt) 
C                        - to calculate correction for a station                                                                                     C                                                                         
C    Parematers: 
C    character  
C          phase : a  string specifying the PHASE,   -e.g P, ScP etc.  
C                                                        
C    real 
C          edist  :  epicentral distance to station (in degrees)     
C          edepth :  depth of event         
C          ecolat :  epicentral co-latitude of source (in radians) 
C          azim   :  azimuth from source to station (in radians)
C                                
C          tcor   :  time correction for path to allow for ellipticity
C 
C    logical 
C          abrt   :  a logical variable -usally set to .FALSE.  
C                    which is set to .TRUE. if a phase for      
C                    which no data is available is chosen       
C                                                                         
C==========================================================================
C   B.L.N. Kennett RSES,ANU        May 1995                 
C   (based on earlier routine by D.J. Brown)                                                                    C=========================================================================
      character *(*) phase
      character*8 phcod(55)
      integer phind(55),phspn(55),phnch(55)
      real edist,edepth,ecolat,azim,
     ^     sc0,sc1,sc2,s3,tcor,
     ^     tau0, a0,b0,h0,d0,e0,f0,g0,
     ^     tau1, a1,b1,h1,d1,e1,f1,g1,
     ^     tau2, a2,b2,h2,d2,e2,f2,g2
      real dpth(6),delta(50)
      real t0(50,6),t1(50,6),t2(50,6)
      integer Ne,Nd
      logical abrt
      data phcod/
     & "Pup   ","P     ","PKPab ","PKPbc ","PKPdf ","PKiKP ",
     & "pP    ","pPKPab","pPKPbc","pPKPdf","pPKiKP","sP    ",
     & "sPKPab","sPKPbc","sPKPdf","sPKiKP","PcP   ","ScP   ",
     & "SKPab ","SKPbc ","SKPdf ","SKiKP ","PKKPab","PKKPbc",
     & "PKKPdf","SKKPab","SKKPbc","SKKPdf","PP    ","P'P'  ",
     & "Sup   ","S     ","SKSac ","SKSdf ","pS    ","pSKSac",
     & "pSKSdf","sS    ","sSKSac","sSKSdf","ScS   ","PcS   ",
     & "PKSab ","PKSbc ","PKSdf ","PKKSab","PKKSbc","PKKSdf",
     & "SKKSac","SKKSdf","SS    ","S'S'  ","SP    ","PS    ",
     & "PnS   "/
      data phind/
     &        1,      14,      95,     124,     137,     194,
     &      323,     392,     421,     434,     491,     620,
     &      701,     730,     743,     800,     929,    1006,
     &     1059,    1072,    1093,    1154,    1275,    1296,
     &     1341,    1466,    1475,    1532,    1661,    1786,
     &     1891,    1904,    1981,    2046,    2111,    2148,
     &     2209,    2270,    2339,    2404,    2465,    2542,
     &     2595,    2608,    2625,    2686,    2695,    2752,
     &     2881,    3054,    3187,    3312,    3501,    3570,
     &     3611/
      data phspn/
     &        3,      20,       7,       3,      14,      32,
     &       17,       7,       3,      14,      32,      20,
     &        7,       3,      14,      32,      19,      13,
     &        3,       5,      15,      30,       5,      11,
     &       31,       2,      14,      32,      31,      26,
     &        3,      19,      16,      16,       9,      15,
     &       15,      17,      16,      15,      19,      13,
     &        3,       4,      15,       2,      14,      32,
     &       43,      33,      31,      47,      17,      10,
     &        6/
      data phnch/
     &        3,       1,       5,       5,       5,       5,
     &        2,       6,       6,       6,       6,       2,
     &        6,       6,       6,       6,       3,       3,
     &        5,       5,       5,       5,       6,       6,
     &        6,       6,       6,       6,       2,       4,
     &        3,       1,       5,       5,       2,       6,
     &        6,       2,       6,       6,       3,       3,
     &        5,       5,       5,       6,       6,       6,
     &        6,       6,       2,       4,       2,       2,
     &        3/
      data dpth/ 0.0, 100.0, 200.0, 300.0, 500.0, 700.0 /
      save sc0,sc1,sc2
c...
c<sc>
c	                     initial call to set up source dependent constants
      entry ellref(ecolat)
c                                            
      s3 = sqrt(3.0)/2.0
      sc0 = 0.25*(1.0+3.0*cos(2.0*ecolat))
      sc1 = s3*sin(2.0*ecolat)
      sc2 = s3*sin(ecolat)*sin(ecolat)
      return
c<sc>
c<ec>                                           phase identification
      entry ellcor(phase,edist,edepth,ecolat,azim,tcor,abrt)
*      write(6,*) "phase,edist,edepth,ecolat,azim"
*      write(6,*)  phase,edist,edepth,ecolat,azim
      Nd = 6
      deldst = 5.0
      abrt = .FALSE.
c                                             select phase
      ip = -1
      do 10 i=1,54
        nc = phnch(i)
        nc1 = nc+1
        if (phase(1:nc) .eq. phcod(i)(1:nc)) then
          if(phase(nc1:nc1).eq."c".or.phase(nc1:nc1).eq."K") go to 10
          ip = i
          Ne = phspn(ip)
          go to 11
        endif
 10   continue
 11   continue
*      write(6,*) "ip:",ip
c                                              phase not found
      if(ip.lt.0) then
        write(6,*)' Selected phase is not available'
        abrt = .true.
        return
      endif
c                                               special case of upgoing waves
*
c
c                                                acquire phase information
       nr = phind(ip)
*       write(6,*) "nrec:",nr
       read(21,61,rec=nr) phcod(ip),np,d1,d2
*       write(6,*) "phcode,np,d1,d2: ", phcod(ip),np,d1,d2
       nr = nr+1
       if(np.ne.Ne) write(6,*) "HELP! - index wrong"
       do 15 i=1,np
         read(21,62,rec=nr) delta(i)
         nr = nr+1
         read(21,63,rec=nr) (t0(i,m),m=1,6)
         nr = nr+1
         read(21,63,rec=nr) (t1(i,m),m=1,6)
         nr = nr+1
         read(21,63,rec=nr) (t2(i,m),m=1,6)
         nr = nr+1
 15    continue         
 61    format(a8,i10,2f10.0)
 62    format(f10.0)
 63    format(6f10.4)
c                                  distance index
       idist = 1 + int( (edist-d1)/ deldst )
       if(edist.lt.d1) idist =1
       if(edist.gt.d2) idist= np-1
c                                  depth index
       do 25 j = 1,Nd-1
         if ((dpth(j).le.edepth).and.(dpth(j+1).ge.edepth))then
            jdepth = j
            goto 26
         endif
 25    continue
 26    continue
*       write(6,*) "idist, jdepth;",idist,jdepth
c
*                      need to allow for zero entries (where phase
*                      description strongly depth dependent)
c tau0
         a0 = t0(idist,jdepth)
         b0 = t0(idist,jdepth+1)
         h0 = t0(idist+1,jdepth+1)
         d0 = t0(idist+1,jdepth)
         e0 = a0 + 
     ^       (d0-a0)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         f0 = b0 + 
     ^       (h0-b0)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         g0 = e0 + 
     ^       (f0-e0)*(edepth-dpth(jdepth))/(dpth(jdepth+1)-dpth(jdepth))
         tau0 = g0
c tau1
         a1 = t1(idist,jdepth)
         b1 = t1(idist,jdepth+1)
         h1 = t1(idist+1,jdepth+1)
         d1 = t1(idist+1,jdepth)
         e1 = a1 + 
     ^       (d1-a1)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         f1 = b1 + 
     ^       (h1-b1)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         g1 = e1 + 
     ^       (f1-e1)*(edepth-dpth(jdepth))/(dpth(jdepth+1)-dpth(jdepth))
         tau1 = g1
c tau2
         a2 = t2(idist,jdepth)
         b2 = t2(idist,jdepth+1)
         h2 = t2(idist+1,jdepth+1)
         d2 = t2(idist+1,jdepth)
         e2 = a2 + 
     ^       (d2-a2)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         f2 = b2 + 
     ^       (h2-b2)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         g2 = e2 + 
     ^       (f2-e2)*(edepth-dpth(jdepth))/(dpth(jdepth+1)-dpth(jdepth))
         tau2 = g2
c
*         write(6,*) "tau0,tau1,tau2:",tau0,tau1,tau2
         caz = cos(azim)
         cbz = cos(2.0*azim)
*         write(6,*) "azim,caz,cbz",azim,caz,cbz    
c
         tcor = sc0*tau0 + sc1*cos(azim)*tau1 + sc2*cos(2.0*azim)*tau2
c
      return
c<ec>
      end
