C 
C Subroutine file 'rm.for' to support program 'greenss.for'.
C
C Version: 5.20
C Date: 1998, September 6
C
C=======================================================================
C
      SUBROUTINE RM(NLAY,IFREE,ITER,OMEGA,TR,P1,P2,P3,
     *                  AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3)
      INTEGER NLAY,IFREE,ITER
      REAL OMEGA,TR,P1,P2,P3
      REAL AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3
C
C Subroutine to calculate the response of the transmission through the
C stack of fine horizontal layers at the receiver site.  The parameters
C of the layer stack should have already been stored in the memory by
C the invocation of subroutine INPUT.  The frequency-dependent
C transmission coefficiens are calculated by subroutines of package
C RMATRIX
C written by C. J. Thomson (Copyright (c) 1996 C. J. Thomson).
C
C Input:
C     NLAY... Number of layers, including the bottom and top halfspaces.
C     IFREE...IFREE=1: The thickness THICKN(1) of the first (top) layer
C               is used for the calculation of the travel time
C               correction at the geophone situated in the top halfspace
C               at the elevation THICKN(1) above the uppermost
C               interface.
C             IFREE=0: Has the same effect as THICKN(1)=0.
C     ITER... ITER=1: Polarization vectors in the layers are calculated.
C             ITER.GT.1: Polarization vectors from the previous
C               invocation are used.
C     OMEGA...Angular frequency.
C     TR...   Travel time at the receiver.
C     P1,P2,P3... Slowness vector at the receiver.
C     AR1,AI1,AR2,AI2,AR3,AI3... Complex-valued displacement at the
C             bottom of the layer stack.
C
C Output:
C     TR...   Travel time at the bottom of the layer stack.
C     BR1,BI1,BR2,BI2,BR3,BI3... Complex-valued displacement at the
C             geophone (IFREE=1) or at the top of the layer stack
C             (IFREE=0).
C
C Date: 1998, September 4
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Common blocks for RMATRIX:
      INTEGER NLAYMX
      PARAMETER (NLAYMX=30)
      INTEGER ISOFLG
      DOUBLE PRECISION THICKN,A,RHO
      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)
C
C     Input and output arguments for RMATRIX:
      COMPLEX*16 PX,PY,COMEGA
      COMPLEX*16 TUS(3,3),RUS(3,3),TDS(3,3),RDS(3,3)
C
C     Temporary storage locations:
      COMPLEX*16 E11,E21,E31,E12,E22,E32,E13,E23,E33
      COMPLEX    F11,F21,F31,F12,F22,F32,F13,F23,F33,E,A1,A2,A3,B1,B2,B3
C
C-----------------------------------------------------------------------
C
C
C     Travel-time correction corresponding to the elevation of the
C     bottom of the layer stack with respect to the receiver:
      TR=TR+THICKN(NLAY)*ABS(P3)
C
C     Calculation of the transmission matrices:
      COMEGA=DCMPLX(DBLE(OMEGA),0.0D0)
      PX=DCMPLX(DBLE(P1),0.0D0)
      PY=DCMPLX(DBLE(P2),0.0D0)
      CALL RMATRIX(PX,PY,COMEGA,NLAY,IFREE,TUS,RUS,TDS,RDS,ITER)
C
C     Complex-valued amplitude at the bottom of the layer stack
      A1=CMPLX(AR1,AI1)
      A2=CMPLX(AR2,AI2)
      A3=CMPLX(AR3,AI3)
C
C     Transforming the amplitude from Cartesian coordinates to
C     the eigenvectors Eij
      E11=EVECS(1,4,NLAY)
      E21=EVECS(2,4,NLAY)
      E31=EVECS(3,4,NLAY)
      E12=EVECS(1,5,NLAY)
      E22=EVECS(2,5,NLAY)
      E32=EVECS(3,5,NLAY)
      E13=EVECS(1,6,NLAY)
      E23=EVECS(2,6,NLAY)
      E33=EVECS(3,6,NLAY)
C     Fij=inv(Eij)*det(E)
      F11=E22*E33-E23*E32
      F12=E32*E13-E33*E12
      F13=E12*E23-E13*E22
      F21=E23*E31-E21*E33
      F22=E33*E11-E31*E13
      F23=E13*E21-E11*E23
      F31=E21*E32-E22*E31
      F32=E31*E12-E32*E11
      F33=E11*E22-E12*E21
C     E=det(E)
      E=E11*F11+E12*F21+E13*F31
C     Bi=inv(Eij)*Ai
      B1=(F11*A1+F12*A2+F13*A3)/E
      B2=(F21*A1+F22*A2+F23*A3)/E
      B3=(F31*A1+F32*A2+F33*A3)/E
C
C     Transition through the layer stack
      A1=TUS(1,1)*B1+TUS(1,2)*B2+TUS(1,3)*B3
      A2=TUS(2,1)*B1+TUS(2,2)*B2+TUS(2,3)*B3
      A3=TUS(3,1)*B1+TUS(3,2)*B2+TUS(3,3)*B3
C
C     Transforming the amplitude from the eigenvectors to the
C     Cartesian coordinates
      B1=EVECS(1,4,1)*A1+EVECS(1,5,1)*A2+EVECS(1,6,1)*A3
      B2=EVECS(2,4,1)*A1+EVECS(2,5,1)*A2+EVECS(2,6,1)*A3
      B3=EVECS(3,4,1)*A1+EVECS(3,5,1)*A2+EVECS(3,6,1)*A3
C
C     Complex-valued amplitude at the top of the layer stack
      BR1=REAL (B1)
      BI1=AIMAG(B1)
      BR2=REAL (B2)
      BI2=AIMAG(B2)
      BR3=REAL (B3)
      BI3=AIMAG(B3)
C
      RETURN
      END
C
C=======================================================================
C