      SUBROUTINE CRSCUPD(H1,H2,M,EGVL,EGVR,CLR,CRL,MX,
     1                  BMATL,BMATR,W4,ALPHAR,BETAR,FNORMR,
     2                  ALPHAL,BETAL,FNORML,FKWL,FKWR,ROHBL,ROHBR,
     3                  NPBL,DEPBL,RHOBL,GRHOBL,LHBL,
     4                  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 IS GENERALIZED TO ALLOW
C     THE BOTTOM DENSITY TO BE DIFFERENT ON EACH SIDE OF THE VERTICAL
C     INTERFACE AND EXIBIT EXPONENTIAL VARIATIONS IN LAYERS. 
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION DEPBL(MX),RHOBL(MX),GRHOBL(MX)
      DIMENSION DEPBR(MX),RHOBR(MX),GRHOBR(MX)
      COMPLEX*16  CI,ONE
      COMPLEX*16  CDSRT,CIW,CEXIB,CEXIWB
      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
      LHBLM1=LHBL-1
      LHBRM1=LHBR-1
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,ROHW,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,ROHW,ROHBR)
C
   50 CONTINUE
C     
      IF(H1 .GT. H2) THEN
C
C     THIS IS THE UPSLOPE CASE
C
C     FIND THE INDEX LRBH1 SUCH THAT H1 IS BETWEEN DEPBR(LRBH1)
C     AND DEPBR(LRBH1+1)
C
      LRBH1=0
      NPBRM1=NPBR-1
      DO 100 L=1,NPBRM1
      IF(DEPBR(L+1) .GE. H1) GO TO 200
  100 CONTINUE
C     
C      H1 IS .GT. DEPBR(NPBR)
C
      LRBH1=NPBR
      GO TO 250
C
  200 CONTINUE
      LRBH1=L
  250 CONTINUE
      LRBH1M1=LRBH1-1
      LRBH1P1=LRBH1+1
C
C      WRITE(9,*) 'UPSLOPE CASE'
C      WRITE(9,*) 'H2= ',H2
C      IF(LRBH1M1 .GE. 1) THEN
C      WRITE(9,*) 'LRBH1M1, DEPBR(LRBH1M1)= ',LRBH1M1, DEPBR(LRBH1M1)
C      END IF
C      WRITE(9,*) 'LRBH1, DEPBR(LRBH1)= ',LRBH1, DEPBR(LRBH1)
C      WRITE(9,*) 'H1= ',H1
C      IF(LRBH1P1 .LE. NPBR) THEN
C      WRITE(9,*) 'LRBH1P1, DEPBR(LRBH1P1)= ',LRBH1P1, DEPBR(LRBH1P1)
C      END IF
C      WRITE(9,*) 'LHBR, DEPBR(LHBR)= ',LHBR, DEPBR(LHBR)
C      WRITE(9,*) 'NPBR,HB= ',NPBR,HB
C
      DO 550 I=1,M
      DO 500 II=1,M
C
      CLR(II,I)=CIW(H1,H2,ALPHAL(I),ALPHAR(II))
      CRL(II,I)=CLR(II,I)
C
      GRHO0=0.0D0
      CLR(II,I)=CLR(II,I)+CEXIWB(H1,H2,H2,H1,H2,ROHW,GRHO0,
     1          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C     
      IF(LHBLM1 .LT. 1 ) GO TO 320
      DO 300 L=1,LHBLM1
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,DEPBL(L),DEPBL(L+1),
     1          DEPBL(L),RHOBL(L),GRHOBL(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
  300 CONTINUE
  320 CONTINUE
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,DEPBL(LHBL),HB,
     1          DEPBL(LHBL),RHOBL(LHBL),GRHOBL(LHBL),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      IF(LRBH1M1 .LT. 1 ) GO TO 420
      DO 400 L=1,LRBH1M1
      CRL(II,I)=CRL(II,I)+CEXIWB(H1,H2,DEPBR(L),DEPBR(L+1),
     1          DEPBR(L),RHOBR(L),GRHOBR(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
  400 CONTINUE
  420 CONTINUE
C
      CRL(II,I)=CRL(II,I)+CEXIWB(H1,H2,DEPBR(LRBH1),H1,
     1          DEPBR(LRBH1),RHOBR(LRBH1),GRHOBR(LRBH1),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      IF(LRBH1 .EQ. LHBR) THEN
C
C     H1 IS IN THE LAYER WITH HB
C
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,H1,HB,
     1          DEPBR(LHBR),RHOBR(LHBR),GRHOBR(LHBR),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      ELSE
C      
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,H1,DEPBR(LRBH1P1),
     1          DEPBR(LRBH1),RHOBR(LRBH1),GRHOBR(LRBH1),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      IF(LHBRM1 .LT. LRBH1P1) GO TO 460
      DO 450 L=LRBH1P1,LHBRM1
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,DEPBR(L),DEPBR(L+1),
     1          DEPBR(L),RHOBR(L),GRHOBR(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
  450 CONTINUE
  460 CONTINUE
C
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,DEPBR(LHBR),HB,
     1          DEPBR(LHBR),RHOBR(LHBR),GRHOBR(LHBR),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      END IF
C
      V=CDSRT(FNORML(I)*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 LLBH2 SUCH THAT H2 IS BETWEEN DEPBL(LLBH2)
C     AND DEPBL(LLBH2+1)
C
      LLBH2=0
      NPBLM1=NPBL-1
      DO 600 L=1,NPBLM1
      IF(DEPBL(L+1) .GE. H2) GO TO 700
  600 CONTINUE
C     
C      H2 IS .GT. DEPBL(NPBL)
C
      LLBH2=NPBL
      GO TO 750
C
  700 CONTINUE
      LLBH2=L
  750 CONTINUE
      LLBH2M1=LLBH2-1
      LLBH2P1=LLBH2+1
C
C      WRITE(9,*) 'DOWN SLOPE CASE'
C      WRITE(9,*) 'H1= ',H1
C      IF(LLBH2M1 .GE. 1) THEN
C      WRITE(9,*) 'LLBH2M1, DEPBL(LLBH2M1)= ',LLBH2M1, DEPBL(LLBH2M1)
C      END IF
C      WRITE(9,*) 'LLBH2, DEPBL(LLBH2)= ',LLBH2, DEPBL(LLBH2)
C      WRITE(9,*) 'H2= ',H2
C      IF(LLBH2P1 .LE. NPBL) THEN
C      WRITE(9,*) 'LLBH2P1, DEPBL(LLBH2P1)= ',LLBH2P1, DEPBL(LLBH2P1)
C      END IF
C      WRITE(9,*) 'LHBL, DEPBL(LHBL)= ',LHBL, DEPBL(LHBL)
C      WRITE(9,*) 'NPBL,HB= ',NPBL,HB
C
      DO 1050 I=1,M
      DO 1000 II=1,M
C
      CLR(II,I)=CIW(H1,H2,ALPHAL(I),ALPHAR(II))
      CRL(II,I)=CLR(II,I)
C
      GRHO0=0.0D0
      CRL(II,I)=CRL(II,I)+CEXIWB(H2,H1,H1,H2,H1,ROHW,GRHO0,
     1          ALPHAR(II),ALPHAL(I),BETAR(II),BETAL(I))
C     
      IF(LLBH2M1 .LT. 1 ) GO TO 820
      DO 800 L=1,LLBH2M1
      CLR(II,I)=CLR(II,I)+CEXIWB(H2,H1,DEPBL(L),DEPBL(L+1),
     1          DEPBL(L),RHOBL(L),GRHOBL(L),
     2          ALPHAR(II),ALPHAL(I),BETAR(II),BETAL(I))
  800 CONTINUE
  820 CONTINUE
      CLR(II,I)=CLR(II,I)+CEXIWB(H2,H1,DEPBL(LLBH2),H2,
     1          DEPBL(LLBH2),RHOBL(LLBH2),GRHOBL(LLBH2),
     2          ALPHAR(II),ALPHAL(I),BETAR(II),BETAL(I))
C
      IF(LLBH2 .EQ. LHBL) THEN
C
C     H2 IS IN THE LAYER WITH HB
C
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,H2,HB,
     1          DEPBL(LHBL),RHOBL(LHBL),GRHOBL(LHBL),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      ELSE
C
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,H2,DEPBL(LLBH2P1),
     1          DEPBL(LLBH2),RHOBL(LLBH2),GRHOBL(LLBH2),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      IF(LHBLM1 .LT. LLBH2P1) GO TO 920
      DO 900 L=LLBH2P1,LHBLM1
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,DEPBL(L),DEPBL(L+1),
     1          DEPBL(L),RHOBL(L),GRHOBL(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
  900 CONTINUE
  920 CONTINUE
C
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,DEPBL(LHBL),HB,
     1          DEPBL(LHBL),RHOBL(LHBL),GRHOBL(LHBL),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      END IF
C
      IF(LHBRM1 .LT. 1 ) GO TO 960
      DO 950 L=1,LHBRM1
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,DEPBR(L),DEPBR(L+1),
     1          DEPBR(L),RHOBR(L),GRHOBR(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
  950 CONTINUE
  960 CONTINUE
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,DEPBR(LHBR),HB,
     1          DEPBR(LHBR),RHOBR(LHBR),GRHOBR(LHBR),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
C
      V=CDSRT(FNORML(I)*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
C      WRITE(9,*) 'NO SLOPE CASE'
C      WRITE(9,*) 'H1= ',H1
C      WRITE(9,*) 'H2= ',H2
C      WRITE(9,*) 'LHBL, DEPBL(LHBL)= ',LHBL, DEPBL(LHBL)
C      WRITE(9,*) 'NPBL,HB= ',NPBL,HB
C      WRITE(9,*) 'LHBR, DEPBR(LHBR)= ',LHBR, DEPBL(LHBR)
C      WRITE(9,*) 'NPBR,HB= ',NPBR,HB
C
      DO 1550 I=1,M
      DO 1500 II=1,M
C
      CLR(II,I)=CIW(H1,H2,ALPHAL(I),ALPHAR(II))
      CRL(II,I)=CLR(II,I)
C
      IF(LHBLM1 .LT. 1 ) GO TO 1220
      DO 1200 L=1,LHBLM1
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,DEPBL(L),DEPBL(L+1),
     1          DEPBL(L),RHOBL(L),GRHOBL(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
 1200 CONTINUE
 1220 CONTINUE
      CLR(II,I)=CLR(II,I)+CEXIB(H1,H2,DEPBL(LHBL),HB,
     1          DEPBL(LHBL),RHOBL(LHBL),GRHOBL(LHBL),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      IF(LHBRM1 .LT. 1) GO TO 1360
      DO 1350 L=1,LHBRM1
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,DEPBR(L),DEPBR(L+1),
     1          DEPBR(L),RHOBR(L),GRHOBR(L),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
 1350 CONTINUE
 1360 CONTINUE
      CRL(II,I)=CRL(II,I)+CEXIB(H1,H2,DEPBR(LHBR),HB,
     1          DEPBR(LHBR),RHOBR(LHBR),GRHOBR(LHBR),
     2          ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
      V=CDSRT(FNORML(I)*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
