DIF Radix-2 FFT Algorithm

Below is the Fortran code for a Decimation-in-Frequency, Radix-2, three butterfly Cooley-Tukey FFT followed by a bit-reversing unscrambler.

C   A COOLEY-TUKEY RADIX 2, DIF  FFT PROGRAM
C   THREE-BF, MULT BY  1  AND  J  ARE REMOVED
C   COMPLEX INPUT DATA IN ARRAYS X AND Y
C   TABLE LOOK-UP OF W VALUES
C      C. S. BURRUS, RICE UNIVERSITY, SEPT 1983
C---------------------------------------------------------
    SUBROUTINE FFT (X,Y,N,M,WR,WI)
    REAL X(1), Y(1), WR(1), WI(1)
C--------------MAIN FFT LOOPS-----------------------------
C
    N2 = N
    DO 10 K = 1, M
        N1 = N2
        N2 = N2/2
        JT = N2/2 + 1
        DO 1 I = 1, N, N1
        L = I + N2
        T    = X(I) - X(L)
        X(I) = X(I) + X(L)
        X(L) = T
        T    = Y(I) - Y(L)
        Y(I) = Y(I) + Y(L)
        Y(L) = T
   1        CONTINUE
        IF (K.EQ.M) GOTO 10
        IE  = N/N1
        IA  = 1
        DO 20 J = 2, N2
        IA = IA + IE
        IF (J.EQ.JT) GOTO 50
        C = WR(IA)
        S = WI(IA)
        DO 30 I = J, N, N1
            L = I + N2
            T    = X(I) - X(L)
            X(I) = X(I) + X(L)
            TY   = Y(I) - Y(L)
            Y(I) = Y(I) + Y(L)
            X(L) = C*T + S*TY
            Y(L) = C*TY - S*T
   30       CONTINUE
        GOTO 25
   50       DO 40 I = J, N, N1
            L = I + N2
            T    = X(I) - X(L)
            X(I) = X(I) + X(L)
            TY   = Y(I) - Y(L)
            Y(I) = Y(I) + Y(L)
            X(L) = TY
            Y(L) =-T
   40       CONTINUE
   25       A = J*E
   20       CONTINUE
   10   CONTINUE
C------------DIGIT REVERSE COUNTER Goes here----------
    RETURN
    END