c  program to calculate the reflection and transmission coefficents
c  of elastic waves across boundaries or at free surfaces.
c  made by V. Maupin for G385, with a subroutine from Cormier and
c  Richards.

      parameter (nnrt=6,nnhs=100)
c... description of the model, and waves
      complex r1,a1,b1,r2,a2,b2,hs
      complex vs(8),vw(8)
      integer m(8)
      common/coin/vs,vw,m
c... reflection and transmission coefficients
      complex slow
      integer icoef(nnrt)
      complex  sd(nnrt)
      integer msd(nnrt)
      common/coout/sd,msd
      complex value
      real amp(nnrt,nnhs),fas(nnrt,nnhs)
c...
      character*16 filename
      character*12 intftyp
      character*8 wavetype
      character*4 answer   

c.. specify and open output file (OP, 1/11-06)
      write (*,*) 'Enter filename for output: '
      read (*,'(a)') filename
      open (UNIT=2,FILE=filename,STATUS='UNKNOWN')

c.. definition of the model
      write (*,*) 
     &  'Type of interface? 1:solid-solid,2:free,3:liquid-solid'
      read(*,*) interf

c.. (OP, 1/11-06)
      if (interf.eq.1) intftyp(1:12)='solid-solid '
      if (interf.eq.2) intftyp(1:12)='free surface'
      if (interf.eq.3) intftyp(1:12)='liquid-solid'

      write (*,*) 'Default model? (y/n)'
      read(*,'(a)') answer  
      if (answer(1:1).eq.'y') then
        rho1=3.4
        alp1=8.1
        bet1=4.5
        rho2=2.9
        alp2=6.8
        bet2=3.9
      else
        if (interf.eq.1) 
     &   write (*,*) 'enter density, alpha, beta in lower layer'
        if (interf.eq.2) 
     &   write (*,*) 'enter density, alpha, beta in half-space'
        if (interf.eq.3) 
     &   write (*,*) 'enter density, alpha, beta in solid layer'
        read(*,*) rho1,alp1,bet1
        if (interf.eq.1)  then
         write (*,*) 'enter density, alpha, beta in upper layer'
         read(*,*) rho2,alp2,bet2
        endif
        if (interf.eq.3) then
         write (*,*) 'enter density, alpha, in liquid layer'
         read(*,*) rho2,alp2
         bet2=0.
        endif
      endif

      r1=cmplx(rho1,0.)
      a1=cmplx(alp1,0.)
      b1=cmplx(bet1,0.)
      r2=cmplx(rho2,0.)
      a2=cmplx(alp2,0.)
      b2=cmplx(bet2,0.)
c
c.. choice of incident wave
      write (*,*) 'Incident wave type? example: "SH down"'
      read(*,'(a)') wavetype

c.. (OP, 2/11-06)
      write (2,1002) intftyp,wavetype
 1002 format (' Interface: ',a12,'; incident wave: ',a8)
 1003 format (' No. of waves generated: ',i1)

c.. determine interface type for subroutine coeff
      if (wavetype(1:2).eq.'SH'.and.interf.eq.1) kint=6
      if (wavetype(1:2).eq.'SH'.and.interf.eq.2) kint=8
      if (wavetype(1:2).eq.'SH'.and.interf.eq.3) kint=7
      if (wavetype(1:2).eq.'SV'.and.interf.eq.1) kint=4
      if (wavetype(1:2).eq.'SV'.and.interf.eq.2) kint=5
      if (wavetype(1:2).eq.'SV'.and.interf.eq.3) kint=2
      if (wavetype(1:2).eq.'P '.and.interf.eq.1) kint=4
      if (wavetype(1:2).eq.'P '.and.interf.eq.2) kint=5
      if (wavetype(1:2).eq.'P '.and.interf.eq.3) kint=2

c.. determine waves coming out
      if (kint.eq.6) then
        ntype=2
        write (2,1003) ntype
        if (wavetype(4:5).eq.'up') then
          icoef(1)=23
          icoef(2)=24
          write (*,*) 'sh upward transmitted'
          write (*,*) 'sh downward reflected'
          write (2,'(a)') 'SH transm up'
          write (2,'(a)') 'SH refl down'
        else
          icoef(1)=21
          icoef(2)=22
          write (*,*) 'sh upward reflected  '
          write (*,*) 'sh downward transmitted'
          write (2,'(a)') 'SH refl up'
          write (2,'(a)') 'SH transm down'
        endif
      endif

      if (kint.eq.8.or.kint.eq.7) then
        ntype=1
        write (2,1003) ntype
        if (wavetype(4:5).eq.'up') then
          icoef(1)=24
          write (*,*) 'sh downward reflected'
          write (2,'(a)') 'SH refl down'
        else
          stop 'no downgoing SH wave in this case'
        endif
      endif

      if (kint.eq.4) then
        ntype=4
        write (2,1003) ntype
        if (wavetype(1:5).eq.'SV up') then
          icoef(1)=13
          icoef(2)=14
          icoef(3)=15
          icoef(4)=16
          write (*,*) 'p upward transmitted'
          write (*,*) 'sv upward transmitted'
          write (*,*) 'p downward reflected'
          write (*,*) 'sv downward reflected'
          write (2,'(a)') 'P transm up'
          write (2,'(a)') 'SV transm up'
          write (2,'(a)') 'P refl down'
          write (2,'(a)') 'SV refl down'
        endif
        if (wavetype(1:7).eq.'SV down') then
          icoef(1)=5 
          icoef(2)=6 
          icoef(3)=7 
          icoef(4)=8 
          write (*,*) 'p upward reflected  '
          write (*,*) 'sv upward reflected  '
          write (*,*) 'p downward transmitted'
          write (*,*) 'sv downward transmitted'
          write (2,'(a)') 'P refl up'
          write (2,'(a)') 'SV refl up'
          write (2,'(a)') 'P transm down'
          write (2,'(a)') 'SV transm down'
        endif
        if (wavetype(1:4).eq.'P up') then
          icoef(1)=9 
          icoef(2)=10
          icoef(3)=11
          icoef(4)=12
          write (*,*) 'p upward transmitted'
          write (*,*) 'sv upward transmitted'
          write (*,*) 'p downward reflected'
          write (*,*) 'sv downward reflected'
          write (2,'(a)') 'P transm up'
          write (2,'(a)') 'SV transm up'
          write (2,'(a)') 'P refl down'
          write (2,'(a)') 'SV refl down'
        endif
        if (wavetype(1:6).eq.'P down') then
          icoef(1)=1
          icoef(2)=2
          icoef(3)=3
          icoef(4)=4
          write (*,*) 'p upward reflected  '
          write (*,*) 'sv upward reflected ' 
          write (*,*) 'p downward transmitted'
          write (*,*) 'sv downward transmitted'
          write (2,'(a)') 'P refl up'
          write (2,'(a)') 'SV refl up' 
          write (2,'(a)') 'P transm down'
          write (2,'(a)') 'SV transm down'
        endif
      endif

      if (kint.eq.5) then
        ntype=2
        write (2,1003) ntype
        if (wavetype(1:5).eq.'SV up') then
          icoef(1)=15
          icoef(2)=16
          write (*,*) 'p downward reflected'
          write (*,*) 'sv downward reflected'
          write (2,'(a)') 'P refl down'
          write (2,'(a)') 'SV refl down'
        endif
        if (wavetype(1:4).eq.'P up') then
          icoef(1)=11
          icoef(2)=12
          write (*,*) 'p downward reflected'
          write (*,*) 'sv downward reflected'
          write (2,'(a)') 'P refl down'
          write (2,'(a)') 'SV refl down'
        endif
        if (wavetype(4:7).eq.'down'.or.wavetype(3:6).eq.'down') 
     &            stop 'no downgoing incident wave in this case'
      endif

      if (kint.eq.2) then
        ntype=3
        write (2,1003) ntype
        if (wavetype(1:5).eq.'SV up') then
          icoef(1)=13
          icoef(2)=15
          icoef(3)=16
          write (*,*) 'p upward transmitted'
          write (*,*) 'p downward reflected'
          write (*,*) 'sv downward reflected'
          write (2,'(a)') 'P transm up'
          write (2,'(a)') 'P refl down'
          write (2,'(a)') 'SV refl down'
        endif
        if (wavetype(1:7).eq.'SV down') 
     &       stop 'no downgoing incident S wave in this case'
        if (wavetype(1:4).eq.'P up') then
          icoef(1)=9 
          icoef(2)=11
          icoef(3)=12
          write (*,*) 'p upward transmitted'
          write (*,*) 'p downward reflected'
          write (*,*) 'sv downward reflected'
          write (2,'(a)') 'P transm up'
          write (2,'(a)') 'P refl down'
          write (2,'(a)') 'SV refl down'
        endif
        if (wavetype(1:6).eq.'P down') then
          icoef(1)=1
          icoef(2)=3
          icoef(3)=4
          write (*,*) 'p upward reflected  '
          write (*,*) 'p downward transmitted'
          write (*,*) 'sv downward transmitted'
          write (2,'(a)') 'P refl up'
          write (2,'(a)') 'P transm down'
          write (2,'(a)') 'SV transm down'
        endif
      endif

c.. definition of the slowness interval
c   we choose it such that the incoming wave is not inhomogeneous
      hsmin=0.
      if (wavetype(1:6).eq.'P down')  hsmax=1./alp2
      if (wavetype(1:4).eq.'P up')  hsmax=1./alp1
      if (wavetype(1:7).eq.'SV down')  hsmax=1./bet2
      if (wavetype(1:5).eq.'SV up')  hsmax=1./bet1
      if (wavetype(1:7).eq.'SH down')  hsmax=1./bet2
      if (wavetype(1:5).eq.'SH up')  hsmax=1./bet1-0.00001
      nhs=41.
      dhs=(hsmax-hsmin)/(nhs-1)

c.. wave amplitudes (equal to 1 in our case)
      do i=1,8
        vw(i)=1.
        m(i)=1
      enddo

c.. calculate the coefficients as a function of slowness

      do ihs=1,nhs
        hs=(ihs-1)*dhs
        vs(1)=slow(hs,1./alp1)
        vs(3)=vs(1)
        vs(2)=slow(hs,1./bet1)
        vs(4)=vs(2)
        if (interf.eq.2) goto 10
        vs(5)=slow(hs,1./alp2)
        vs(7)=vs(5)
        if (interf.eq.3) goto 10
        vs(6)=slow(hs,1./bet2)
        vs(8)=vs(6)
 10     continue

        call coeff(r1,a1,b1,r2,a2,b2,hs,kint,icoef,ntype)

        do itype=1,ntype
          value=sd(itype)*exp(float(msd(itype)))
          amp(itype,ihs)=cabs(value)
          if (value.eq.(0.,0.)) then
            fas(itype,ihs)=0.
          else
            fas(itype,ihs)=atan2(imag(value),real(value))
          endif
        enddo
      enddo

c..(OP, 1/11-06)
c      write (*,*) 'enter filename for output'
c      read (*,'(a)') filename
c      open (UNIT=2,FILE=filename,STATUS='UNKNOWN')

      write (2,*) nhs
      do ihs=1,nhs
        hs=(ihs-1)*dhs
        write (2,1001) 
     &       real(hs),(amp(itype,ihs),fas(itype,ihs),itype=1,ntype)
      enddo
 1001 format (9e12.4)
      close (2)
      
      stop
      end
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      complex function slow(hs,slowave)
c calculate the vertical slowness, given the horizontal slowness hs
c  and the wave slowness slowave
c
      if (hs.lt.slowave) then
        slow =cmplx(sqrt(slowave*slowave-hs*hs),0.)
      else
        slow =cmplx(0.,sqrt(hs*hs-slowave*slowave))
      endif
      return
      end
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
      subroutine coeff(r1,a1,b1,r2,a2,b2,hs,kint,icoef,ntype)
c        this subroutine has been written by Cormier and Richards,
c        and is extracted from their program for spectral synthesis
c        of body waves, described in the book "Seismolgical
c        algorithms" ed by Doornbos. (additional comment by VM)
c
c        this subroutine selects and evaluates from among 25 useful
c     reflection/conversion/transmission/correction coefficients for a
c     spherically symmetric earth.
c
c        r1,a1,b1 are the density and p- and s-wave speeds just below
c     the interface, and r2,a2,b2 are these elastic constants just
c     above the interface. note that this is a natural ordering
c     in terms of radius, and is opposite to an ordering with depth.
c
c        hs is the horizontal slowness (in units which are reciprocal
c     to a1,...,b2). for a spherical earth, use of ray parameter (sec/
c     radian) is via hs = ray parameter/radius of interface. in fact,
c     this subroutine can be used unchanged for flat earth problems.
c
c        kint describes the kind of interface problem,
c            =1 for solid over fluid (p-sv)
c            =2 for fluid over solid (p-sv)
c            =3 for fluid over fluid (p-sv)
c            =4 for solid over solid (p-sv)
c            =5 for a stress-free upper surface (p/sv)
c            =6 for solid over solid (sh)
c            =7 for solid over fluid (sh)
c            =8 for a shear-free upper surface (sh).
c
c        icoef is an array of up to ntype integers, ntype.le.20
c    (or, from dimensions in the main program, ntype.le.(e.g.)6),
c     and each integer in the array will cause a specific coefficient
c     to be evaluated. integers stored in icoef must either be in the
c     range 1,...,20 (which will generate p/sv coefficients) or within
c     21,...,25 (to generate sh coefficients). thus, p/sv and sh
c     problems cannot be studied together in the same call to this
c     subroutine.
c
c        the answers generated are placed in the arrays sd(complex)
c     and msd(integer), which are in common/coout/. these arrays
c     have length 25. since only a few of these locations may be
c     filled by the subroutine, do initialising elsewhere.
c
c        some input to the subroutine is provided by common/coin/,
c     which contains up to eight generalized vertical slownesses, vs,
c     and up to eight vertical wavefunctions, vw*e**m. these too should
c     have components which are initialized to zero if they are not
c     applicable to the coefficients which this subroutine is called
c     to evaluate. to get the notation, use the following table
c
c     xidown1  --  vs(1)          pidown1  --  vw(1)*e**m(1)
c     etadown1  -  vs(2)          sigdown1  -  vw(2)*e**m(2)
c     xiup1  ----  vs(3)          piup1  ----  vw(3)*e**m(3)
c     etaup1  ---  vs(4)          sigup1  ---  vw(4)*e**m(4)
c     xidown2  --  vs(5)          pidown2  --  vw(5)*e**m(5)
c     etadown2  -  vs(6)          sigdown2  -  vw(6)*e**m(6)
c     xiup2  ----  vs(7)          piup2  ----  vw(7)*e**m(7)
c     etaup2  ---  vs(8)          sigup2  ---  vw(8)*e**m(8).
c
c        in spirit, this subroutine rests on work of nafe(bssa,
c     pages 205-219, july 1957), who figured out a neat way to generate
c     all coefficients relevant to the same interface, and scholte
c     (koninkl. med. meteorol. inst. publ. 65, pages 1-55, . good
c     luck with finding this one), who saw the advantages of taking the
c     turning points, defined for each of the two layers on either side
c     of the interface, as the reference level.
c
c        i use sign conventions as in chapters 5 and 9 of
c     aki and richards, *methods of quantitative seismology*.
c
c        in summary, the subroutine yields
c
c         sd(k), msd(k), k=1,ntype
c
c     where the coefficients in sd are defined by the following table
c
c***********************************************************************
c
c  p incident downwards, to  p upward reflected  ----  sd( 1)*e**msd( 1)
c  p incident downwards, to sv upward reflected  ----  sd( 2)*e**msd( 2)
c  p incident downwards, to  p downward transmitted    sd( 3)*e**msd( 3)
c  p incident downwards, to sv downward transmitted    sd( 4)*e**msd( 4)
c
c sv incident downwards, to  p upward reflected  ----  sd (5)*e**msd( 5)
c sv incident downwards, to sv upward reflected  ----  sd( 6)*e**msd( 6)
c sv incident downwards, to  p downward transmitted    sd( 7)*e**msd( 7)
c sv incident downwards, to sv downward transmitted    sd( 8)*e**msd( 8)
c
c  p incident upwards, to  p upward transmitted  ----  sd( 9)*e**msd( 9)
c  p incident upwards, to sv upward transmitted  ----  sd(10)*e**msd(10)
c  p incident upwards, to  p downward reflected  ----  sd(11)*e**msd(11)
c  p incident upwards, to sv downward reflected  ----  sd(12)*e**msd(12)
c
c sv incident upwards, to  p upward transmitted  ----  sd(13)*e**msd(13)
c sv incident upwards, to sv upward transmitted  ----  sd(14)*e**msd(14)
c sv incident upwards, to  p downward reflected  ----  sd(15)*e**msd(15)
c sv incident upwards, to sv downward reflected  ----  sd(16)*e**msd(16)
c
c
c  p incident upwards, to horizontal displacement
c                      of free surface  -------------  sd(17)*e**msd(17)
c  p incident upwards, to vertical(upward) displacement
c                      of free surface  -------------  sd(18)*e**msd(18)
c
c sv incident upwards, to horizontal displacement
c                      of free surface  -------------  sd(19)*e**msd(19)
c sv incident upwards, to vertical(upward) displacement
c                      of free surface  -------------  sd(20)*e**msd(20)
c
c***********************************************************************
c
c     now follow the sh coefficients
c
c sh incident downwards, to sh upward reflected  ----  sd(21)*e**msd(21)
c sh incident downwards, to sh downward transmitted    sd(22)*e**msd(22)
c
c sh incident upwards, to sh upward transmitted  ----  sd(23)*e**msd(23)
c sh incident upwards, to sh downward reflected  ----  sd(24)*e**msd(24)
c
c sh incident upwards, to horizontal displacement
c                      of free surface  -------------  sd(25)*e**msd(25)
c
c
      parameter 
     &   (nnrt=6)
      integer icoef(nnrt)
      complex  sd(nnrt), vs(8),vw(8)
      integer m(8),msd(nnrt)
      complex r1,a1,b1,r2,a2,b2,c2,hs,hssq,a,b,c,d,e,f,g,h,delta
     1  ,rig1,rig2
      common/coin/vs,vw,m
      common/coout/sd,msd
      data c2/(2.,0.)/
      hssq=hs**2
      rig1=r1*b1**2
      rig2=r2*b2**2
      d=c2*(rig2-rig1)
      b=r2-d*hssq
      c=r1+d*hssq
      a=b-r1
c.....
      go to (101,102,103,104,105,106,107,108),kint
c.....
  101 continue
c....................     solid over fluid
      e=b*vs(1)+c*vs(7)
      f=b
      g=a-d*vs(1)*vs(8)
      h=-d*vs(7)
      delta=e*f+g*h*hssq
      do 100 k=1,ntype
      i=icoef(k)
      go to (1,2,3,20,5,6,7,20,9,10,11,20,20,20,20,20),i
    1 sd(k)=-(((b*vs(1)-c*vs(5))*f+d*vs(5)*g*hssq)/delta)
     1   *vw(5)/vw(7)
      msd(k)=m(5)-m(7)
      go to 99
    2 sd(k)=-((vs(7)+vs(5))*b*d*vs(1)*hs*a2/b2/delta)*vw(5)/vw(8)
      msd(k)=m(5)-m(8)
      go to 99
    3 sd(k)=(r2*(vs(7)+vs(5))*f*a2/a1/delta)*vw(5)/vw(1)
      msd(k)=m(5)-m(1)
      go to 99
    5 sd(k)=((vs(8)+vs(6))*b*d*vs(1)*hs*b2/a2/delta)*vw(6)/vw(7)
      msd(k)=m(6)-m(7)
      go to 99
    6 sd(k)=-((b*e+(a+d*vs(1)*vs(6))*h*hssq)/delta)*vw(6)/vw(8)
      msd(k)=m(6)-m(8)
      go to 99
    7 sd(k)=(r2*(vs(8)+vs(6))*h*hs*b2/a1/delta)*vw(6)/vw(1)
      msd(k)=m(6)-m(1)
      go to 99
    9 sd(k)=(r1*(vs(1)+vs(3))*f*a1/a2/delta)*vw(3)/vw(7)
      msd(k)=m(3)-m(7)
      go to 99
   10 sd(k)=-(r1*(vs(1)+vs(3))*h*hs*a1/b2/delta)*vw(3)/vw(8)
      msd(k)=m(3)-m(8)
      go to 99
   11 sd(k)=(((b*vs(3)-c*vs(7))*f-(a+d*vs(3)*vs(8))*h*hssq)/delta)
     1   *vw(3)/vw(1)
      msd(k)=m(3)-m(1)
   99 continue 
      go to 100
   20 write(* ,2001)  i
 2001 format('you asked for  coefficient number ',i3,' for fluid'
     1 ,' below')
  100 continue
      return
c.....
  102 continue
c....................     fluid over solid
      e=b*vs(1)+c*vs(7)
      f=c
      g=-d*vs(1)
      h=a-d*vs(7)*vs(2)
      delta=e*f+g*h*hssq
      do 200 k=1,ntype
      i=icoef(k)
      go to (21,40,23,24,40,40,40,40,29,40,31,32,33,40,35,36),i
   21 sd(k)=-(((b*vs(1)-c*vs(5))*f+(a+d*vs(5)*vs(2))*g*hssq)/delta)
     1   *vw(5)/vw(7)
      msd(k)=m(5)-m(7)
      go to 199
   23 sd(k)=(r2*(vs(7)+vs(5))*f*a2/a1/delta)*vw(5)/vw(1)
      msd(k)=m(5)-m(1)
      go to 199
   24 sd(k)=-(r2*(vs(7)+vs(5))*g*hs*a2/b1/delta)*vw(5)/vw(2)
      msd(k)=m(5)-m(2)
      go to 199
   29 sd(k)=(r1*(vs(1)+vs(3))*f*a1/a2/delta)*vw(3)/vw(7)
      msd(k)=m(3)-m(7)
      go to 199
   31 sd(k)=(((b*vs(3)-c*vs(7))*f-d*vs(3)*h*hssq)/delta)*vw(3)/vw(1)
      msd(k)=m(3)-m(1)
      go to 199
   32 sd(k)=-((vs(1)+vs(3))*c*d*vs(7)*hs*a1/b1/delta)*vw(3)/vw(2)
      msd(k)=m(3)-m(2)
      go to 199
   33 sd(k)=(r1*(vs(2)+vs(4))*g*hs*b1/a2/delta)*vw(4)/vw(7)
      msd(k)=m(4)-m(7)
      go to 199
   35 sd(k)= ((vs(2)+vs(4))*c*d*vs(7)*hs*b1/a1/delta)*vw(4)/vw(1)
      msd(k)=m(4)-m(1)
      go to 199
   36 sd(k)=-((c*e+(a+d*vs(7)*vs(4))*g*hssq)/delta)*vw(4)/vw(2)
      msd(k)=m(4)-m(2)
  199 continue
      go to 200
   40 write(* ,2002)  i
 2002 format('you asked for coefficient number ',i3,' for fluid'
     1 ,' above')
  200 continue
      return
c.....
  103 continue
c....................     fluid over fluid
      e=b*vs(1)+c*vs(7)
      do 300 k=1,ntype
      i=icoef(k)
      go to (41,60,43,60,60,60,60,60,49,60,51,60,60,60,60,60),i
   41 sd(k)=-(b*vs(1)-c*vs(5))/e*vw(5)/vw(7)
      msd(k)=m(5)-m(7)
      go to 299
   43 sd(k)=r2*(vs(7)+vs(5))*a2/a1/e*vw(5)/vw(1)
      msd(k)=m(5)-m(1)
      go to 299
   49 sd(k)=r1*(vs(1)+vs(3))*a1/a2/e*vw(3)/vw(7)
      msd(k)=m(3)-m(7)
      go to 299
   51 sd(k)=(b*vs(3)-c*vs(7))/e*vw(3)/vw(1)
      msd(k)=m(3)-m(1)
  299 continue
      go to 300
   60 write(* ,2003)  i
 2003 format('you asked for coefficient number ',i3,' with fluid'
     1  ,' over fluid')
  300 continue
      return
c.....
  104 continue
c....................     solid over solid
      e=b*vs(1)+c*vs(7)
      f=b*vs(2)+c*vs(8)
      g=a-d*vs(1)*vs(8)
      h=a-d*vs(7)*vs(2)
      delta=e*f+g*h*hssq
      do 400 k=1,ntype
      i=icoef(k)
      go to (61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76),i
   61 sd(k)=-(((b*vs(1)-c*vs(5))*f+(a+d*vs(5)*vs(2))*g*hssq)/delta)
     1   *vw(5)/vw(7)
      msd(k)=m(5)-m(7)
      go to 399
   62 sd(k)=-((vs(7)+vs(5))*(a*c+b*d*vs(1)*vs(2))*hs*a2/b2/delta)
     1   *vw(5)/vw(8)
      msd(k)=m(5)-m(8)
      go to 399
   63 sd(k)=(r2*(vs(7)+vs(5))*f*a2/a1/delta)
     1   *vw(5)/vw(1)
      msd(k)=m(5)-m(1)
      go to 399
   64 sd(k)=-(r2*(vs(7)+vs(5))*g*hs*a2/b1/delta)
     1   *vw(5)/vw(2)
      msd(k)=m(5)-m(2)
      go to 399
   65 sd(k)=((vs(8)+vs(6))*(a*c+b*d*vs(1)*vs(2))*hs*b2/a2/delta)
     1   *vw(6)/vw(7)
      msd(k)=m(6)-m(7)
      go to 399
   66 sd(k)=-(((b*vs(2)-c*vs(6))*e+(a+d*vs(1)*vs(6))*h*hssq)/delta)
     1   *vw(6)/vw(8)
      msd(k)=m(6)-m(8)
      go to 399
   67 sd(k)=(r2*(vs(8)+vs(6))*h*hs*b2/a1/delta)
     1   *vw(6)/vw(1)
      msd(k)=m(6)-m(1)
      go to 399
   68 sd(k)=(r2*(vs(8)+vs(6))*e*b2/b1/delta)
     1   *vw(6)/vw(2)
      msd(k)=m(6)-m(2)
      go to 399
   69 sd(k)=(r1*(vs(1)+vs(3))*f*a1/a2/delta)
     1   *vw(3)/vw(7)
      msd(k)=m(3)-m(7)
      go to 399
   70 sd(k)=-(r1*(vs(1)+vs(3))*h*hs*a1/b2/delta)
     1   *vw(3)/vw(8)
      msd(k)=m(3)-m(8)
      go to 399
   71 sd(k)=(((b*vs(3)-c*vs(7))*f-(a+d*vs(3)*vs(8))*h*hssq)/delta)
     1   *vw(3)/vw(1)
      msd(k)=m(3)-m(1)
      go to 399
   72 sd(k)=-((vs(1)+vs(3))*(a*b+c*d*vs(7)*vs(8))*hs*a1/b1/delta)
     1   *vw(3)/vw(2)
      msd(k)=m(3)-m(2)
      go to 399
   73 sd(k)= (r1*(vs(2)+vs(4))*g*hs*b1/a2/delta)
     1   *vw(4)/vw(7)
      msd(k)=m(4)-m(7)
      go to 399
   74 sd(k)= (r1*(vs(2)+vs(4))*e*b1/b2/delta)
     1   *vw(4)/vw(8)
      msd(k)=m(4)-m(8)
      go to 399
   75 sd(k)= ((vs(2)+vs(4))*(a*b+c*d*vs(7)*vs(8))*hs*b1/a1/delta)
     1   *vw(4)/vw(1)
      msd(k)=m(4)-m(1)
      go to 399
   76 sd(k)= (((b*vs(4)-c*vs(8))*e-(a+d*vs(7)*vs(4))*g*hssq)/delta)
     1   *vw(4)/vw(2)
      msd(k)=m(4)-m(2)
  399 continue
  400 continue
      return
c.....
  105 continue
c....................     free surface
      delta=c**2+d**2*vs(1)*vs(2)*hssq
      do 500 k=1,ntype
      i=icoef(k)
      go to (80,80,80,80,80,80,80,80,80,80,81,82,80,80,85,86,
     1  87,88,89,90),i
   81 sd(k)=(-c**2+d**2*vs(3)*vs(2)*hssq)/delta*vw(3)/vw(1)
      msd(k)=m(3)-m(1)
      go to 499
   82 sd(k)=-(vs(1)+vs(3))*c*d*hs*a1/b1/delta*vw(3)/vw(2)
      msd(k)=m(3)-m(2)
      go to 499
   85 sd(k)=-(vs(2)+vs(4))*c*d*hs*b1/a1/delta*vw(4)/vw(1)
      msd(k)=m(4)-m(1)
      go to 499
   86 sd(k)=(c**2-d**2*vs(1)*vs(4)*hssq)/delta*vw(4)/vw(2)
      msd(k)=m(4)-m(2)
      go to 499
   87 sd(k)=-r1*a1*hs*d*(vs(1)+vs(3))*vs(2)/delta
      go to 499
   88 sd(k)=r1*a1*c*(vs(1)+vs(3))/delta
      go to 499
   89 sd(k)=r1*b1*c*(vs(2)+vs(4))/delta
      go to 499
   90 sd(k)=r1*b1*hs*d*(vs(2)+vs(4))*vs(1)/delta
  499 continue
      go to 500
   80 write(* ,2005)  i
 2005 format('you asked for coefficient number ',i3,' for a free'
     1  ,' surface')
  500 continue
      return
c.....
  106 continue
c....................     general sh problem
      delta=rig1*vs(2)+rig2*vs(8)
      do 600 k=1,ntype
      i=icoef(k)
      im=i-20
      go to (91,92,93,94),im
   91 sd(k)=(rig2*vs(6)-rig1*vs(2))/delta*vw(6)/vw(8)
      msd(k)=m(6)-m(8)
      go to 599
   92 sd(k)=rig2*(vs(6)+vs(8))/delta*vw(6)/vw(2)
      msd(k)=m(6)-m(2)
      go to 599
   93 sd(k)=rig1*(vs(2)+vs(4))/delta*vw(4)/vw(8)
      msd(k)=m(4)-m(8)
      go to 599
   94 sd(k)=(rig1*vs(4)-rig2*vs(8))/delta*vw(4)/vw(2)
      msd(k)=m(4)-m(2)
  599 continue
  600 continue
      return
c.....
  107 continue
c....................     solid over fluid (sh)
      do 700 k=1,ntype
      i=icoef(k)
      im=i-20
      go to (701,702,702,702),im
  701 sd(k)=vs(6)/vs(8)*vw(6)/vw(8)
      msd(k)=m(6)-m(8)
      go to 700
  702 write( *,2007) i
 2007 format('you asked for coefficient number ',i3,' for a',
     1  ' solid(sh) over a fluid')
  700 continue
      return
c.....
  108 continue
c....................     shear-free surface above (sh)
      do 800 k=1,ntype
      i=icoef(k)
      im=i-20
      go to (810,810,810,804,805),im
  804 sd(k)=vs(4)/vs(2)*vw(4)/vw(2)
      msd(k)=m(4)-m(2)
      go to 799
  805 sd(k)=(vs(4)+vs(2))/vs(2)
  799 continue
      go to 800
  810 write(* ,2008)  i
 2008 format('you asked for coefficient number ',i3,' for a'
     1  ,' shear-free surface (sh)')
  800 continue
      return
      end