      SUBROUTINE CRSCUP(H1,H2,M,EGVL,EGVR,CLR,CRL,MX,
     1                  IMAT,BMATL,BMATR,W4,ALPHAR,BETAR,FNORMR,
     2                  ALPHAL,BETAL,FNORML,FKWL,FKWR,ROHBL,ROHBR)
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. THE FUNCTIONS CIW0 AND CIB0 USED HERE HAVE NOT BEEN 
C     DIVIDED BY THE CORESPONDING DENSITIES ROHW, ROHBL AND ROHBR.
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16  CI,ONE
      COMPLEX*16  CDSRT,CIW0,CIB0,CIWB
      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  CHW0,CHWB,CHB0,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   400   I=1,M
      ALPHAR(I)=CDSRT(FKWR2+FKB2*(EGVR(I)-ONE))
C
      IF(IMAT .GE. 1) THEN
      BETAR(I)=CDSRT(FKB2*EGVR(I))
      ELSE
      BETAR(I)=CDSRT(FKB2*(EGVR(I)+CI*ATEN))
      END IF
C
      FNORMR(I)=GAMA00(H2,ALPHAR(I),BETAR(I),HB,FKB,ROHW,ROHBR)
      ALPHAL(I)=CDSRT(FKWL2+FKB2*(EGVL(I)-ONE))
C
      IF(IMAT .GE. 1) THEN
      BETAL(I)=CDSRT(FKB2*EGVL(I))
      ELSE
      BETAL(I)=CDSRT(FKB2*(EGVL(I)+CI*ATEN))
      END IF
C
      FNORML(I)=GAMA00(H1,ALPHAL(I),BETAL(I),HB,FKB,ROHW,ROHBL)
  400 CONTINUE
C
      DO 700 II=1,M
C
      DO 600 I=1,M
C
      CHW0=CIW0(H1,H2,ALPHAL(I),ALPHAR(II))
      CHB0=CIB0(H1,H2,ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
C
C     GHANGED ORDER OF THE SQUARE ROOT AND PRODUCT 6/8/96
C     FROM V=CDSRT(FNORML(I))*FNORMR(II)) TO:
C
      V=CDSRT(FNORML(I))*CDSRT(FNORMR(II))
C
      IF(H2 .GE. H1) GO TO 550
C
C     THIS IS THE UPSLOPE CASE.
C
      CHWB=CIWB(H1,H2,ALPHAL(I),ALPHAR(II),BETAL(I),BETAR(II))
      CLR(II,I)=(CHW0/ROHW+CHWB/ROHW +CHB0/ROHBL)/V
      CRL(II,I)=(CHW0/ROHW+CHWB/ROHBR+CHB0/ROHBR)/V
C
      GO TO 600
C
  550 CONTINUE
C
C     THIS IS THE DOWNSLOPE CASE.
C
      CHWB=CIWB(H2,H1,ALPHAR(II),ALPHAL(I),BETAR(II),BETAL(I))
      CLR(II,I)=(CHW0/ROHW+CHWB/ROHBL+CHB0/ROHBL)/V
      CRL(II,I)=(CHW0/ROHW+CHWB/ROHW +CHB0/ROHBR)/V
C
  600 CONTINUE
  700 CONTINUE
C
      IF(IMAT .LT. 1) GO TO 1000
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
 1000 CONTINUE
C
      RETURN
      END
