      SUBROUTINE CRSCUPB(H1,H2,M,EGVL,EGVR,CLR,CRL,MX,
     1                  BMATL,BMATR,W4,ALPHAR,BETAR,FNORMR,
     2                  ALPHAL,BETAL,FNORML,FKWL,FKWR,
     3                  ROHWL,NPWL,DEPWL,RHOWL,GRHOWL,LHWL,
     4                  ROHWR,NPWR,DEPWR,RHOWR,GRHOWR,LHWR,
     5                  ROHBL,NPBL,DEPBL,RHOBL,GRHOBL,LHBL,
     6                  ROHBR,NPBR,DEPBR,RHOBR,GRHOBR,LHBR)
C
C     THE FOLLOWING SUBROUTINE COMPUTES THE CROSS COUPLING BETWEEN THE
C     EIGEN FUNCTIONS CORRESPONDING TO EGVL AND THE EIGENFUNCTIONS
C     CORRESPONDING TO EGVR. THIS VERSION INCORPORATES A DENSITY 
C     PROFILE ON EACH SIDE OF THE VERTICAL INTERFACE THAT EXIBITS 
C     PIECEWISE EXPONENTIAL VARIATION IN BOTH THE WATER AND BOTTOM 
C     SEDIMENT LAYERS. 
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION DEPWL(MX),RHOWL(MX),GRHOWL(MX)
      DIMENSION DEPWR(MX),RHOWR(MX),GRHOWR(MX)
      DIMENSION DEPBL(MX),RHOBL(MX),GRHOBL(MX)
      DIMENSION DEPBR(MX),RHOBR(MX),GRHOBR(MX)
      COMPLEX*16  CI,ONE
      COMPLEX*16  CDSRT,CUPW,CUPB,CUPWB,CUPBW
      COMPLEX*16  BMATL(MX,MX),BMATR(MX,MX),W4(MX,MX)
      COMPLEX*16  EGVL(MX),EGVR(MX),CLR(MX,MX),CRL(MX,MX)
      COMPLEX*16  GAMA00,V
      COMPLEX*16  ALPHAR(MX), BETAR(MX), FNORMR(MX)
      COMPLEX*16  ALPHAL(MX), BETAL(MX), FNORML(MX)
C
      COMMON /BLKEVN/ HB,CW,CB,FKW,FKB,ROHW,ROHB,ATEN
C
      ONE=DCMPLX(1.0,0.0)
      CI=DCMPLX(0.0,1.0)
      FKWL2=FKWL**2
      FKWR2=FKWR**2
      FKB2=FKB**2
C
      DO 50 I=1,M
C
      ALPHAL(I)=CDSRT(FKWL2+FKB2*(EGVL(I)-ONE))
      BETAL(I)=CDSRT(FKB2*EGVL(I))
      FNORML(I)=GAMA00(H1,ALPHAL(I),BETAL(I),HB,FKB,ROHWL,ROHBL)
C
      ALPHAR(I)=CDSRT(FKWR2+FKB2*(EGVR(I)-ONE))
      BETAR(I)=CDSRT(FKB2*EGVR(I))
      FNORMR(I)=GAMA00(H2,ALPHAR(I),BETAR(I),HB,FKB,ROHWR,ROHBR)
C
   50 CONTINUE
C     
      IF(H1 .GT. H2) THEN
C
C     THIS IS THE UPSLOPE CASE
C
C     FIND THE INDEX LLWH2 SUCH THAT H2 IS BETWEEN DEPWL(LLWH2)
C     AND DEPWL(LLWH2+1)
C
      LLWH2=INDEP(H2,DEPWL,NPWL,MX)
C
C     FIND THE INDEX LRBH1 SUCH THAT H1 IS BETWEEN DEPBR(LRBH1)
C     AND DEPBR(LRBH1+1)
C
      LRBH1=INDEP(H1,DEPBR,NPBR,MX)
C
      DO 550 I=1,M
      DO 500 II=1,M
C
C     UPSLOPE CASE
C
      CLR(II,I)=CUPWB(H1,H2,ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II),
     1          DEPWL,RHOWL,GRHOWL,LHWL,LLWH2,MX)
      CLR(II,I)=CLR(II,I)+CUPB(H1,H2,ALPHAL(I),ALPHAR(II),
     1          BETAL(I),BETAR(II),DEPBL,RHOBL,GRHOBL,LHBL,MX)
C
      CRL(II,I)=CUPW(H2,ALPHAL(I),ALPHAR(II),
     1               DEPWR,RHOWR,GRHOWR,LHWR,MX)
      CRL(II,I)=CRL(II,I)+CUPBW(H1,H2,ALPHAL(I),ALPHAR(II),
     1          BETAL(I),BETAR(II),DEPBR,RHOBR,GRHOBR,LHBR,LRBH1,MX)
C
C     CHANGED ORDER SQUARE ROOT AND MULTIPLICATION 6/8/96
C     FROM V=CDSRT(FNORML(I)*FNORMR(II)) TO:
C
      V=CDSRT(FNORML(I))*CDSRT(FNORMR(II))
      CLR(II,I)=CLR(II,I)/V
      CRL(II,I)=CRL(II,I)/V
C
  500 CONTINUE
  550 CONTINUE
C
      ELSE IF(H1 .LT. H2) THEN
C
C     THIS IS THE DOWNSLOPE CASE (H1 .LT. H2).
C
C     FIND THE INDEX LRWH1 SUCH THAT H1 IS BETWEEN DEPWR(LRWH1)
C     AND DEPWR(LRWH1+1)
C
      LRWH1=INDEP(H1,DEPWR,NPWR,MX)
C
C     FIND THE INDEX LLBH2 SUCH THAT H2 IS BETWEEN DEPBL(LLBH2)
C     AND DEPBL(LLBH2+1)
C
      LLBH2=INDEP(H2,DEPBL,NPBL,MX)
C
      DO 1050 I=1,M
      DO 1000 II=1,M
C
C     DOWN SLOPE CASE
C
      CLR(II,I)=CUPW(H1,ALPHAL(I),ALPHAR(II),
     1               DEPWL,RHOWL,GRHOWL,LHWL,MX)
      CLR(II,I)=CLR(II,I)+CUPBW(H2,H1,ALPHAR(II),ALPHAL(I),
     1          BETAR(II),BETAL(I),DEPBL,RHOBL,GRHOBL,LHBL,LLBH2,MX)
C
      CRL(II,I)=CUPWB(H2,H1,ALPHAR(II),ALPHAL(I),BETAR(II),BETAL(I),
     1          DEPWR,RHOWR,GRHOWR,LHWR,LRWH1,MX)
      CRL(II,I)=CRL(II,I)+CUPB(H1,H2,ALPHAL(I),ALPHAR(II),
     1          BETAL(I),BETAR(II),DEPBR,RHOBR,GRHOBR,LHBR,MX)
C
      V=CDSRT(FNORML(I))*CDSRT(FNORMR(II))
      CLR(II,I)=CLR(II,I)/V
      CRL(II,I)=CRL(II,I)/V
C
 1000 CONTINUE
 1050 CONTINUE
C
      ELSE
C
C     THIS IS THE NO SLOPE CASE (H1 .EQ. H2)
C
      DO 1550 I=1,M
      DO 1500 II=1,M
C
C     NO SLOPE CASE
C
      CLR(II,I)=CUPW(H1,ALPHAL(I),ALPHAR(II),
     1               DEPWL,RHOWL,GRHOWL,LHWL,MX)
      CLR(II,I)=CLR(II,I)+CUPB(H1,H1,ALPHAL(I),ALPHAR(II),
     1          BETAL(I),BETAR(II),DEPBL,RHOBL,GRHOBL,LHBL,MX)
      CRL(II,I)=CUPW(H1,ALPHAL(I),ALPHAR(II),
     1               DEPWR,RHOWR,GRHOWR,LHWR,MX)
      CRL(II,I)=CRL(II,I)+CUPB(H1,H1,ALPHAL(I),ALPHAR(II),
     1          BETAL(I),BETAR(II),DEPBR,RHOBR,GRHOBR,LHBR,MX)
C
      V=CDSRT(FNORML(I))*CDSRT(FNORMR(II))
      CLR(II,I)=CLR(II,I)/V
      CRL(II,I)=CRL(II,I)/V
C
 1500 CONTINUE
 1550 CONTINUE
C
      END IF
C
C     THIS IS THE MAIN ADVANTAGE OF THE GALERKIN PROCEDURE IN
C     A COUPLED MODE CALCULATION. THE COUPLING CAN BE DONE 
C     BETWEEN THE BASIC EIGENFUNCTIONS WHERE THE INTEGRALS 
C     ARE DONE EXACTLY, THEN THE TRUE COUPLING MATRIX IS 
C     OBTAINED BY MATRIX MULTIPLICATON
C
      CALL TMMULT(M,BMATR,CLR,BMATL,MX,W4)
      CALL TMMULT(M,BMATR,CRL,BMATL,MX,W4)
C
C
      RETURN
      END
