      subroutine emdlv(r,vp,vs)
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(13)
      common/emdlc/np,rd(13)
      data np,rd/13,1010.,1221.5,1621.5,3480.,3631.,
     1 5611.,5711.,5961.,6161.,
     1 6251.,6336.,6351.,6371./,rn,vn/1.5696123e-4,6.8501006/
      data modnam/'iasp93s'/
c
      call emias5(rn*r,rho,vp,vs)
      vp=vn*vp
      vs=vn*vs
      return
c
      entry emdld(n,cpr,name)
      n=np
      do 1 i=1,np
 1    cpr(i)=rd(i)
      name=modnam
      return
      end
c
      subroutine emias5(x0,ro,vp,vs)
c
c $$$$$ calls no other routine $$$$$
c
c   Emias4 returns model parameters for the IASPEI 93a model (June 
c   1993).  The upper mantle S structure of the IASPEI 91a model has 
c   been modified to make regional S phases faster by about 1 s.  Also, 
c   the SP6 core has been adopted for this version making it a hybred 
c   of the two models.
c
c   Given non-dimensionalized radius x0, emias4 returns
c   non-dimensionalized density, ro, compressional velocity, vp, and
c   shear velocity, vs.  Non-dimensionalization is according to the
c   scheme of Gilbert in program EOS:  x0 by a (the radius of the
c   Earth), ro by robar (the mean density of the Earth), and velocity
c   by a*sqrt(pi*G*robar) (where G is the universal gravitational
c   constant.
c
c
      save
      dimension r(16),p(15,4),s(15,4)
      data r/0. ,1010.,1221.5, 1621.5 ,3480.0  ,3631.,  
     1 5611.   ,5711.   ,
     1 5961.   ,6161.   ,6251.   ,6336.   ,6351.    ,6371.    ,
     2 6371.,6371./
      data p/11.2622, 11.3041, 4.0354, 11.0487,
     1 14.49470,25.1486 ,25.969838,29.38896,
     1 30.78765,25.41389, 8.785412, 6.5   , 5.8   ,2*0.,
     2        0.   ,-1.2730, 82.008, -4.0362,
     2 -1.47089,-41.1538,-16.934118,-21.40656,
     2-23.25415,-17.69722,-0.7495294, 4*0.,
     3  -6.3640, 0., -347.769, 4.8023,
     3  0.0  ,51.9932,9*0.,
     4     0.,   0.     ,468.786, -13.5732,
     4 0.     ,-26.6083,9*0./
      data s/ 3.66780, 3.66780, 0.  , 0.      , 
     1 8.16616,12.9303 ,20.768902,17.70732,
     1 15.24213,5.065646, 5.24541, 3.85   , 3.46   ,2*0.,
     2        0.     , 0.      , 0. , 0.,
     2 -1.58206,-21.2590,-16.531471,-13.50652,
     2-11.08553, -0.566312,-0.749528, 4*0.,
     3       -4.44749, -4.44749, 0.  ,0.    ,
     3  0.0    ,27.8988 ,9*0.,
     4        0.     , 0.      ,  0. , 0. ,
     4 0.     ,-14.1080,9*0./
      data xn,rn,vn/6371.,.18125793,.14598326/,i/1/
c
      x=amax1(x0,0.)
      x1=xn*x
 2    if(x1.ge.r(i)) go to 1
      i=i-1
      go to 2
 1    if(x1.le.r(i+1).or.i.ge.11) go to 3
      i=i+1
      if(i.lt.11) go to 1
 3    vp=vn*(p(i,1)+x*(p(i,2)+x*(p(i,3)+x*p(i,4))))
      vs=vn*(s(i,1)+x*(s(i,2)+x*(s(i,3)+x*s(i,4))))
      return
      end



