C====CALCULATION OF COULOMB MATRIX ELEMENTS==
C ==== (DIRECT,EXCHANGE,COMBINED)==============
      SUBROUTINE MAEL5(AL,BET,NE5,H,L,IT,NK,NM,M,
     *LL,LI,LK,LM,LV,R1,PR,ME,ME4,ME5,ME6,ME7,R8,R2,R5,R6,R7)
      INTEGER NE5,L,IT,NK,NM,M,LL,LI,LK,LM,LV
      REAL*8 C,AL,BET,H,A1,D1,B,PR(*),
     *ME(M,*),ME4(*),ME5(*),ME6(*),ME7(*),
     *R1(*),R2(*),R5(*),R6(*),R7(*),R8(NE5,*)
C      PRINT 2,R8(1,L),R8(1,IT),R8(1,NK),R8(1,NM)
C  2   FORMAT(4E15.8)
      CALL MATR1(AL,BET,NE5,H,1,M,NK,IT,L,NM,LV,ME,
     *R1,R8,R2,R5,R6,R7)
      CALL IOT31(LK,LV,LI,0,0,0,A1)
      CALL IOT31(LM,LV,LL,0,0,0,D1)
      A1=A1*D1*2.0*DSQRT(DFLOAT((2*LK+1)*(2*LI+1)*(2*LM+1)*(2*LL+1)))
      DO 10 I=1,M
 10   ME4(I)=ME(I,1)*A1
C     PRINT 15,LK,LL,LM,LI,LV
C15   FORMAT( ' LK,LL,LM,LI,LV',5I5)
      IF(IABS(LK-LL).GE.IABS(LM-LI)) THEN
      L1=IABS(LL-LK)
      ELSE
      L1=IABS(LI-LM)
      END IF
      IF(IABS(LK+LL).GE.IABS(LM+LI)) THEN
      L2=LM+LI
      ELSE
      L2=LK+LL
      END IF
      C=DSQRT(DFLOAT((2*LK+1)*(2*LL+1)*(2*LM+1)*(2*LI+1)))
      DO 20 I=1,M
      ME5(I)=0.0D 00
 20   CONTINUE
      DO 30 J=L1,L2,2
      J100=J
C     PRINT 76,R8(1,L),R8(1,IT),R8(1,NK),R8(1,NM)
C76   FORMAT(4E15.8)
      CALL MATR1(AL,BET,NE5,H,1,M,NK,L,IT,NM,J100,ME,
     *R1,R8,R2,R5,R6,R7)
      CALL IOT31(LK,J100,LL,0,0,0,A1)
      CALL IOT31(LM,J100,LI,0,0,0,D1)
      CALL IOT61(LK,LV,LI,LM,J100,LL,B)
      PR(J)=A1*D1*B*(-1)**(LV+J)*(2*LV+1)*2
C     PRINT 77,A1,D1,B,PR(J100),L1,L2
C77   FORMAT('A1D1BPRL1L2 ',4E10.4,2I3)
      DO 40 I=1,M
      ME5(I)=ME5(I)+ME(I,1)*PR(J)*C
 40   CONTINUE
 30   CONTINUE
C     PRINT 130,(ME4(I),I=1,M)
C130  FORMAT('DIR  M.E. ',8(E15.8))
C     PRINT 135,(ME5(I),I=1,M)
C135  FORMAT('EXC  M.E. ',8(E15.8))
      DO 50 I=1,M
      ME6(I)=2*ME4(I)-ME5(I)
 50   CONTINUE
C     PRINT 140,(ME6(I),I=1,M)
C140  FORMAT(' COMB.M E      ',8(E15.8))
      DO 60 I=1,M
 60   ME7(I)=ME6(I)**2/DFLOAT(2*LV+1)/DFLOAT(2*LK+1)
      RETURN
      END
