C==COULOMB MATRIX ELEMENTS==(m*f)====
      SUBROUTINE MAEL4(ALFA,BET,NE5,H,NK,IT,F,NM,M,L,
     *        LL,LF,LM,LK,LV,
     *        R1,PR,ME,ME1,ME2,ME3,R8,R2,R5,R6,R7)

       INTEGER NE5,L,NK,NM,M,F,LL,LM,LF,LK,LV,
     *         NE,I,L1,L2,J1,J
       REAL*8 ALFA,BET,H,SUM,A1,D1,B,
     *      ME(M,*),ME1(M,*),ME2(M,*),ME3(M,*),
     *      R1(*),R2(*),R5(*),R6(*),R7(*),
     *      R8(NE5,*),PR(*)
       NE=NE5-15
C       PRINT 5,NK,IT,F,NM,M,L,LL,LF,LM,LK,LV
C 5     FORMAT(' NK,IT,F,NM,M,L,LL,LF,LM,LK,LV'/11I3)
       CALL MATR1(ALFA,BET,NE5,H,F,M,NK,IT,L,NM,LV,
     *          ME,R1,R8,R2,R5,R6,R7)
       CALL IOT31(LK,LV,LF,0,0,0,A1)
       CALL IOT31(LL,LV,LM,0,0,0,D1)
       A1=A1*D1*2*DSQRT(DFLOAT((2*LK+1)*(2*LF+1)*
     *           (2*LL+1)*(2*LM+1)))
      DO 10 I=1, F
      DO 10 J=1, M
       ME1(J,I)=ME(J,I)*A1
10    CONTINUE
      DO 20 I=1, F
      DO 20 J=1, M
       ME2(J,I)=0.00
20    CONTINUE
      IF ( IABS(LK-LM).LE.IABS(LF-LL) ) THEN
         L1=IABS(LF-LL)
      ELSE
         L1=IABS(LK-LM)
      END IF
      IF (IABS(LK+LM).LE.IABS(LF+LL)) THEN
         L2=IABS(LK+LM)
      ELSE
         L2=IABS(LF+LL)
      END IF
       SUM=0
C       PRINT 25,L1,L2,F,M
C25   FORMAT('L1,L2',2I3,'F,M',2I3)
      DO 30 J=L1, L2, 2
      CALL MATR1(ALFA,BET,NE5,H,F,M,L,IT,NK,NM,J,
     *           ME,R1,R8,R2,R5,R6,R7)
      CALL IOT31(LK,J,LM,0,0,0,A1)
      CALL IOT31(LF,J,LL,0,0,0,D1)
      CALL IOT61(LK,LV,LF,LL,J,LM,B)
       PR(J)=A1*D1*B*(-1)**(LV+J)*(2*LV+1)*2
       A1=DSQRT(DFLOAT((2*LK+1)*(2*LF+1)*
     *            (2*LL+1)*(2*LM+1)))
C       PRINT 32,A1
C32   FORMAT('A1',F10.5)
      DO 40 I=1, F
      DO 50 J1=1,M
           ME2(J1,I)=ME2(J1,I)+PR(J)*ME(J1,I)*A1
50    CONTINUE
40    CONTINUE
30    CONTINUE
      DO 60 I=1, F
      DO 70 J=1, M
         ME3(J,I)=2*ME1(J,I)-ME2(J,I)
70    CONTINUE
60    CONTINUE
      RETURN
      END
