C 
      subroutine input(lu,nlay,ifree)

      implicit real*8(a-h,o-z)

c-LK  complex*16 px,py,omega

c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c Rmatrix and associated programs may be used freely by University
c researchers and teachers. Commercial parties interested in using
c the program should contact C. J. Thomson.
c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c Copyright (c) 1996 C. J. Thomson.
c All rights reserved by the author.
c
c-LK: Input logical unit number 10 changed to LU.
c-LK: Dummy arguments PX,PY,OMEGA of subroutine INPUT discarded and
c-LK: replaced with new dummy argument LU.
c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
      parameter (nlaymx=30)
c
c Earth model
c
      common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
     .              rho(nlaymx)

      complex*16 evals,evecs
      common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

      complex*16 q11i,q12i,q21i,q22i
      complex*16 tui,rui,tdi,rdi
      common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
     .               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
     .               tui(3,3,nlaymx),rui(3,3,nlaymx),
     .               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

      complex*16 tusi,rusi,tdsi,rdsi
      common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
     .                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
      complex*16 diag1,j2cst1,j2cstn
      common /symts/diag1(6,nlaymx)
      common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c


      read(lu,*)nlay,ifree
c-LK  read(lu,*)omega,px,py

      do 99 ilay=1,nlay

      read(lu,*)isoflg(ilay),thickn(ilay),rho(ilay)

      if(isoflg(ilay).eq.0)then

        do 100 l=1,3
         do 101 k=1,3

          read(lu,*)((a(i,j,k,l,ilay),i=1,3),j=1,3)

  101    continue
  100   continue

      else

        read(lu,*)a(3,3,3,3,ilay),a(2,3,2,3,ilay)

        iflat=0
        if(iflat.ne.0)then

c-LK     print *,' '
c-LK     print *,'Earth flattening transformation being called.'
c-LK     print *,' '

         call flatten(ilay,thickn(ilay),
     .          a(3,3,3,3,ilay),a(2,3,2,3,ilay),rho(ilay))

        else

c-LK     print *,' '
c-LK     print *,'Earth flattening transformation not being called.'
c-LK     print *,' '

        endif

c next lines for conversion from velocity to elasticity
        a(3,3,3,3,ilay)=(a(3,3,3,3,ilay)**2)
        a(2,3,2,3,ilay)=(a(2,3,2,3,ilay)**2)

      endif

   99 continue

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine flatten(ilay,thickn,vp,vs,rho)

      implicit real*8 (a-h,o-z)

c Copyright (c) 1996 C. J. Thomson.
c All rights reserved by the author.

      save zlast,rad,rad2,a
      data a/6371.d3/
c
      if(ilay.eq.1)then
          rad=a
          rad2=a+thickn/2.d0
          thickl=thickn
          zlast=0.d0
      endif

          rad=rad-thickn
          rad2=rad2-(thickn+thickl)/2.d0
          thickl=thickn

          z=a*dlog(a/rad)
          thickn=z-zlast
          zlast=z
          vp=vp*a/rad2
          vs=vs*a/rad2
          rho=rho*(rad2/a)**5

c          write(11,*)ilay,thickn,vp,vs,rho

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine rmatrix(px,py,omega,nlay,ifree,tus,rus,tds,rds,iter)
c
c Program to compute r/t matrices for a stack of homogeneous
c layers. Uses reflection matrix algorithm of B. L. N. Kennett.
c General anisotropy permitted.
c
c C. J. Thomson, August--November 1993.
c
c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c Rmatrix and associated programs may be used freely by University
c researchers and teachers. Commercial parties interested in using
c the program should contact C. J. Thomson.
c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c Copyright (c) 1996 C. J. Thomson.
c All rights reserved by the author.

      implicit real*8(a-h,o-z)

      complex*16 px,py,omega

c eigenvalue, eigenvector and inverse eigenvector arrays

      complex*16 eval1(6), evec1(6,6), eveci1(6,6)
      complex*16 eval2(6), evec2(6,6), eveci2(6,6)


c partitions of wave propagator

      complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3)

c r/t coefficient matrices - single interface and stack of layers

      complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)
      complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3)

c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
      parameter (nlaymx=30)
c
c Earth model
c
      common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
     .              rho(nlaymx)
c

      complex*16 evals,evecs
      common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

      complex*16 q11i,q12i,q21i,q22i
      complex*16 tui,rui,tdi,rdi
      common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
     .               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
     .               tui(3,3,nlaymx),rui(3,3,nlaymx),
     .               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

      complex*16 tusi,rusi,tdsi,rdsi
      common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
     .                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
      complex*16 diag1,j2cst1,j2cstn
      common /symts/diag1(6,nlaymx)
      common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c

c start at bottom layer
      ilay=nlay

c only compute e-vecs, etc, if on first pass for this slowness
c as they are frequency independent

      if(iter.eq.1)then

       call evset(px,py,omega,ilay,eval1,evec1)

       call estore(eval1,evec1,evals(1,ilay),evecs(1,1,ilay))

      endif

c invert eigenvalue matrix
       call xeveci(nlay,ilay,evecs(1,1,ilay),eveci1)
c      call write2(ilay,evecs(1,1,ilay),eveci1)

c set j2 symmetry constant for bottom half space
       call j2sym(nlay,ilay,evecs(1,1,ilay),j2cstn)

c next layer and loop point
  999 ilay=ilay-1

c check if on first pass for this slowness

      if(iter.eq.1)then

       call evset(px,py,omega,ilay,eval2,evec2)

       call estore(eval2,evec2,evals(1,ilay),evecs(1,1,ilay))

c wave propagator for interface

       call qmatrix(eveci1,evec2,q11,q12,q21,q22)

       call qstore(q11,q12,q21,q22,q11i(1,1,ilay),q12i(1,1,ilay),
     .             q21i(1,1,ilay),q22i(1,1,ilay))

c r/t matrices for interface from wave propagator

       call retra(q11,q12,q21,q22,tu,ru,td,rd)
c      call write3(ilay,tu,ru,td,rd,evec1,evec2)

       call rtstore(tu,ru,td,rd,tui(1,1,ilay),rui(1,1,ilay),
     .              tdi(1,1,ilay),rdi(1,1,ilay))

      else

       call rtget(tu,ru,td,rd,tui(1,1,ilay),rui(1,1,ilay),
     .              tdi(1,1,ilay),rdi(1,1,ilay))

      endif

c addition rule for previous and present r/t matrices
      call addit(nlay,ilay,tus,rus,tds,rds,tu,ru,td,rd)

       call rtstore(tus,rus,tds,rds,tusi(1,1,ilay),rusi(1,1,ilay),
     .              tdsi(1,1,ilay),rdsi(1,1,ilay))

      if(ilay.gt.1)then

c layer up to next interface
c-LK    print *,'entering layer, ilay=',ilay
        call layer(ilay,omega,thickn(ilay),evals(1,ilay),
     .              tus,rus,tds,rds)

        if(iter.eq.1)then

c transfer data
          call trnsf6(eval1,eval2,evec1,evec2)

c invert eigenvalue matrix
          call xeveci(nlay,ilay,evec1,eveci1)
c         call write2(ilay,evec1,eveci1)

        endif

c loop
         go to 999

      endif

c top layer of stack -- if no free surface above just write out results
c-LK    call write4(tus,rus,tds,rds)

c now test some symmetries of stack results... note requires that the
c simple inverse constant of proportionality for the upper half-
c space be found first...

        call xeveci(nlay,ilay,evecs(1,1,ilay),eveci2)
c-LK    call j1symtst(nlay,ilay,tus,rus,tds,rds)

c set j2 symmetry constant for top half space
        call j2sym(nlay,ilay,evecs(1,1,ilay),j2cst1)
c and test j2 symmetry
c-LK    call j2symtst(nlay,ilay,tus,rus,tds,rds)

      if(ifree.eq.1)then

c free surface at top of first layer -- so include phase shift from
c first interface up to free surface

c-LK    print *,'entering layer, ilay=',ilay
        call layer(ilay,omega,thickn(ilay),evals(1,ilay),
     .              tus,rus,tds,rds)

      endif

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
C