SUBROUTINE F1FT(A,B,NTOT,N,NSPAN,ISN,NFAC,NP,AT,CK,BT,SK)
C
C THIS SUBROUTINE COMPUTES THE MIXED-RADIX FOURIER TRANSFORM OR ITS
C INVERSE FOR COMPLEX DATA.  IT CAN COMPUTE A SINGLE VARIATE TRANS-
C FORM WITH ONE CALL STATEMENT OR A MULTIVARIATE TRANSFORM WITH A
C SERIES OF CALL STATEMENTS.
C
C THE METHOD USED IS SINGLETON S ALGORITHM TO EVALUATE
C
C                      N-1
C                 Y(K)=   X(J)EXP(I2(PI)JK/N) K=0,1,. . .,N-1
C                      J=0
C
C WHERE Y(K) AND X(J) ARE BOTH COMPLEX VALUED VECTORS.
C
C
C     INPUT
C
C       A     -VECTOR OF SIZE NTOT CONTAINING REAL ELEMENTS OF DATA
C       B     -VECTOR OF SIZE NTOT CONTAINING IMAGINARY ELEMENTS OF DATA
C      NTOT   -TOTAL NUMBER OF COMPLEX DATA POINTS
C       N     -NUMBER OF COMPLEX DATA POINTS OF THE CURRENT VARIABLE
C     NSPAN  -NUMBER OF ELEMENTS OF BOTH A AND B NECESSARY TO SPAN
C              ALL VALUES OF THE CURRENT VARIABLE
C              (FOR A SINGLE VARIATE TRANSFORM NTOT=N=NSPAN.  FOR A
C              MULTIVARIATE TRANSFORM N AND NSPAN ARE DIFFERENT IN
C              EACH CALL STATEMENT.  SEE EXAMPLE BELOW).
C      ISN   -DETERMINES BY SIGN THE TYPE OF TRANSFORM BEING COMPUTED,
C              FORWARD OR INVERSE
C             -INDICATES BY ABSOLUTE VALUE THE VECTOR ARRANGEMENT OF
C              INPUT DATA (AND TRANSFORM OUTPUT)
C             =+1 COMPLEX INPUT DATA ARE IN TWO VECTORS, A AND B.
C              REAL COMPONENTS IN A, IMAGINARY COMPONENTS IN B.
C             =-1 FOURIER COEFFICIENT INPUT VALUES ARE IN TWO VECTORS,
C              A AND B.  COSINE VALUES IN A, SINE VALUES IN B.
C             =+2 COMPLEX INPUT DATA ARE STORED ALTERNATELY IN A
C              SINGLE COMPLEX VECTOR, A.  REAL VALUES ARE
C              IN A(1),A(3),... AND IMAGINARY VALUES ARE
C              IN A(2),A(4),....  SECOND ARGUMENT OF
C              F1FT SHOULD BE A(2).  (SEE EXAMPLE.)
C             =-2 FOURIER COEFFICIENT INPUT VALUES ARE STORED ALTER-
C              NATELY IN A SINGLE VECTOR, A.  COSINE VALUES ARE IN
C              A(1),A(3),... AND SINE VALUES ARE IN A(2),A(4),....
C              SECOND ARGUMENT OF F1FT SHOULD BE A(2).
C
C
C     OUTPUT
C
C       A     -VECTOR OF SIZE NTOT CONTAINING COSINE COMPONENTS OF
C              FOURIER COEFFICIENTS
C       B     -VECTOR OF SIZE NTOT CONTAINING SINE COMPONENTS OF FOURIER
C              COEFFICIENTS
C
C IF ISN=1 MULTIPLY OUTPUT BY 2/N FOR UNIT MAGNITUDE.
C
C IF ISN=2 OUTPUT WILL ALTERNATE IN A SINGLE VECTOR, A.  COSINE
C COEFFICIENTS IN A(1),A(3),... AND SINE COEFFICIENTS IN A(2),A(4),....
C OUTPUT SHOULD BE MULTIPLIED BY 4/N FOR UNIT MAGNITUDE.
C
C
C NUMBER OF FACTORS OF N MUST NOT EXCEED 11.
C MAXIMUM PRIME FACTOR OF N MUST NOT EXCEED 23.
C PRODUCT OF THE SQUARE-FREE FACTORS, IF THERE EXISTS MORE THAN ONE,
C MUST NOT EXCEED 210.
C
C
C     ERROR STATES
C          1.   PRIME FACTOR .GT. 23
C          2.   SQUARE-FREE FACTOR PRODUCT .GT. 2.  10
C
C     STACK STORAGE - NONE
C
C
C
C EXAMPLES -
C
C TRANSFORM OF N COMPLEX DATA VALUES STORED IN TWO VECTORS -
C                CALL FFT(A,B,N,N,N,1)
C
C INVERSE TRANSFORM OF FOURIER COEFFICIENTS STORED IN TWO VECTORS -
C                CALL FFT(A,B,N,N,N,-1)
C
C TRANSFORM OF N COMPLEX DATA VALUES STORED ALTERNATELY IN A SINGLE
C COMPLEX VECTOR -
C                CALL FFT(A,A(2),2*N,N,2*N,2)
C
C INVERSE TRANSFORM OF N FOURIER COEFFICIENTS, COSINE AND SINE COMPO-
C NENTS STORED ALTERNATELY IN A SINGLE VECTOR -
C                CALL FFT(A,A(2),2*N,N,2*N,-2)
C
C FOR A MULTIVARIATE TRANSFORM THERE IS NO LIMIT ON THE NUMBER OF
C IMPLIED SUBSCRIPTS.  FFT SHOULD BE CALLED ONCE FOR EACH VARIABLE
C AND THE CALLS MAY BE MADE IN ANY ORDER.
C
C TRANSFORM OF A TRI-VARIATE, A(N1,N2,N3), B(N1,N2,N3) -
C                NTOT=N1*N2*N3
C                CALL FFT(A,B,NTOT,N1,N1,1)
C                CALL FFT(A,B,NTOT,N2,N1*N2,1)
C                CALL FFT(A,B,NTOT,N3,NTOT,1)
C
C
C FOR TRANSFORM OF COMPLETELY REAL DATA USE RLTR IN CONJUNCTION WITH
C FFT.
C
C
C OTHER ROUTINES CALLED -  SETERR
C
C
C THE PORT LIBRARY VERSION OF FFT IS AS SHOWN IN THE SINGLETON
C REFERENCE EXCEPT FOR CHANGES IN THE ERROR HANDLING.
C
C
C REFERENCE -  SINGLETON, R. C.,  AN ALGORITHM FOR COMPUTING THE MIXED
C             RADIX FAST FOURIER TRANSFORM , IEEE TRANSACTIONS ON AUDIO
C             AND ELECTROACOUSTICS, VOL. AU-17, NO. 2, JUNE, 1969
C             PP. 93-103.
C
C
      DIMENSION A(1),B(1)
      DIMENSION NFAC(2),NP(2)
      DIMENSION AT(1),CK(1),BT(1),SK(1)
      EQUIVALENCE (I,II)
C  THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH THE ARRAY DIMENSIONS.
C
      MAXF=23
      MAXP=209
      IF(N .LT. 2) RETURN
      RAD=8.0*ATAN(1.0)
      S72=RAD/5.0
      C72=COS(S72)
      S72=SIN(S72)
      S120=SQRT(0.75)
      INC=ISN
      IF(ISN .GE. 0) GO TO 10
      S72=-S72
      S120=-S120
      RAD=-RAD
      INC=-INC
   10 NT=INC*NTOT
      KS=INC*NSPAN
      KSPAN=KS
      NN=NT-INC
      JC=KS/N
      RADF=RAD*FLOAT(JC)*0.5
      I=0
      JF=0
C  DETERMINE THE FACTORS OF N
      M=0
      K=N
      GO TO 20
   15 M=M+1
      NFAC(M)=4
      K=K/16
   20 IF(K-(K/16)*16 .EQ. 0) GO TO 15
      J=3
      JJ=9
      GO TO 30
   25 M=M+1
      NFAC(M)=J
      K=K/JJ
   30 IF(MOD(K,JJ) .EQ. 0) GO TO 25
      J=J+2
      JJ=J**2
      IF(JJ .LE. K) GO TO 30
      IF(K .GT. 4) GO TO 40
      KT=M
      NFAC(M+1)=K
      IF(K .NE. 1) M=M+1
      GO TO 80
   40 IF(K-(K/4)*4 .NE. 0) GO TO 50
      M=M+1
      NFAC(M)=2
      K=K/4
   50 KT=M
      J=2
   60 IF(MOD(K,J) .NE. 0) GO TO 70
      M=M+1
      NFAC(M)=J
      K=K/J
   70 J=((J+1)/2)*2+1
      IF(J .LE. K) GO TO 60
   80 IF(KT .EQ. 0) GO TO 100
      J=KT
   90 M=M+1
      NFAC(M)=NFAC(J)
      J=J-1
      IF(J .NE. 0) GO TO 90
C  COMPUTE FOURIER TRANSFORM
  100 SD=RADF/FLOAT(KSPAN)
      CD=2.0*SIN(SD)**2
      SD=SIN(SD+SD)
      KK=1
      I=I+1
      IF(NFAC(I) .NE. 2) GO TO 400
C  TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)
      KSPAN=KSPAN/2
      K1=KSPAN+2
  210 K2=KK+KSPAN
      AK=A(K2)
      BK=B(K2)
      A(K2)=A(KK)-AK
      B(K2)=B(KK)-BK
      A(KK)=A(KK)+AK
      B(KK)=B(KK)+BK
      KK=K2+KSPAN
      IF(KK .LE. NN) GO TO 210
      KK=KK-NN
      IF(KK .LE. JC) GO TO 210
      IF(KK .GT. KSPAN) GO TO 800
  220 C1=1.0-CD
      S1=SD
  230 K2=KK+KSPAN
      AK=A(KK)-A(K2)
      BK=B(KK)-B(K2)
      A(KK)=A(KK)+A(K2)
      B(KK)=B(KK)+B(K2)
      A(K2)=C1*AK-S1*BK
      B(K2)=S1*AK+C1*BK
      KK=K2+KSPAN
      IF(KK .LT. NT) GO TO 230
      K2=KK-NT
      C1=-C1
      KK=K1-K2
      IF(KK .GT. K2) GO TO 230
      AK=C1-(CD*C1+SD*S1)
      S1=(SD*C1-CD*S1)+S1
C  THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C  ERROR. IF ROUNDED ARITHMETIC IS USED, ONE CAN SUBSTITUTE
C     C1=AK
      C1=0.5/(AK**2+S1**2)+0.5
      S1=C1*S1
      C1=C1*AK
      KK=KK+JC
      IF(KK .LT. K2) GO TO 230
      K1=K1+INC+INC
      KK=(K1-KSPAN)/2+JC
      IF(KK .LE. JC+JC) GO TO 220
      GO TO 100
C  TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)
  320 K1=KK+KSPAN
      K2=K1+KSPAN
      AK=A(KK)
      BK=B(KK)
      AJ=A(K1)+A(K2)
      BJ=B(K1)+B(K2)
      A(KK)=AK+AJ
      B(KK)=BK+BJ
      AK=-0.5*AJ+AK
      BK=-0.5*BJ+BK
      AJ=(A(K1)-A(K2))*S120
      BJ=(B(K1)-B(K2))*S120
      A(K1)=AK-BJ
      B(K1)=BK+AJ
      A(K2)=AK+BJ
      B(K2)=BK-AJ
      KK=K2+KSPAN
      IF(KK .LT. NN) GO TO 320
      KK=KK-NN
      IF(KK .LE. KSPAN) GO TO 320
      GO TO 700
C  TRANSFORM FOR FACTOR OF 4
  400 IF(NFAC(I) .NE. 4) GO TO 600
      KSPNN=KSPAN
      KSPAN=KSPAN/4
  410 C1=1.0
      S1=0
  420 K1=KK+KSPAN
      K2=K1+KSPAN
      K3=K2+KSPAN
      AKP=A(KK)+A(K2)
      AKM=A(KK)-A(K2)
      AJP=A(K1)+A(K3)
      AJM=A(K1)-A(K3)
      A(KK)=AKP+AJP
      AJP=AKP-AJP
      BKP=B(KK)+B(K2)
      BKM=B(KK)-B(K2)
      BJP=B(K1)+B(K3)
      BJM=B(K1)-B(K3)
      B(KK)=BKP+BJP
      BJP=BKP-BJP
      IF(ISN .LT. 0) GO TO 450
      AKP=AKM-BJM
      AKM=AKM+BJM
      BKP=BKM+AJM
      BKM=BKM-AJM
      IF(S1 .EQ. 0.0) GO TO 460
  430 A(K1)=AKP*C1-BKP*S1
      B(K1)=AKP*S1+BKP*C1
      A(K2)=AJP*C2-BJP*S2
      B(K2)=AJP*S2+BJP*C2
      A(K3)=AKM*C3-BKM*S3
      B(K3)=AKM*S3+BKM*C3
      KK=K3+KSPAN
      IF(KK .LE. NT) GO TO 420
  440 C2=C1-(CD*C1+SD*S1)
      S1=(SD*C1-CD*S1)+S1
C  THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C  ERROR. IF ROUNDED ARITHMETIC IS USED, ONE CAN SUBSTITUTE
C     C1=C2
      C1=0.5/(C2**2+S1**2)+0.5
      S1=C1*S1
      C1=C1*C2
      C2=C1**2-S1**2
      S2=2.0*C1*S1
      C3=C2*C1-S2*S1
      S3=C2*S1+S2*C1
      KK=KK-NT+JC
      IF(KK .LE. KSPAN) GO TO 420
      KK=KK-KSPAN+INC
      IF(KK .LE. JC) GO TO 410
      IF(KSPAN .EQ. JC) GO TO 800
      GO TO 100
  450 AKP=AKM+BJM
      AKM=AKM-BJM
      BKP=BKM-AJM
      BKM=BKM+AJM
      IF(S1 .NE. 0.0) GO TO 430
  460 A(K1)=AKP
      B(K1)=BKP
      A(K2)=AJP
      B(K2)=BJP
      A(K3)=AKM
      B(K3)=BKM
      KK=K3+KSPAN
      IF(KK .LE. NT) GO TO 420
      GO TO 440
C  TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)
  510 C2=C72**2-S72**2
      S2=2.0*C72*S72
  520 K1=KK+KSPAN
      K2=K1+KSPAN
      K3=K2+KSPAN
      K4=K3+KSPAN
      AKP=A(K1)+A(K4)
      AKM=A(K1)-A(K4)
      BKP=B(K1)+B(K4)
      BKM=B(K1)-B(K4)
      AJP=A(K2)+A(K3)
      AJM=A(K2)-A(K3)
      BJP=B(K2)+B(K3)
      BJM=B(K2)-B(K3)
      AA=A(KK)
      BB=B(KK)
      A(KK)=AA+AKP+AJP
      B(KK)=BB+BKP+BJP
      AK=AKP*C72+AJP*C2+AA
      BK=BKP*C72+BJP*C2+BB
      AJ=AKM*S72+AJM*S2
      BJ=BKM*S72+BJM*S2
      A(K1)=AK-BJ
      A(K4)=AK+BJ
      B(K1)=BK+AJ
      B(K4)=BK-AJ
      AK=AKP*C2+AJP*C72+AA
      BK=BKP*C2+BJP*C72+BB
      AJ=AKM*S2-AJM*S72
      BJ=BKM*S2-BJM*S72
      A(K2)=AK-BJ
      A(K3)=AK+BJ
      B(K2)=BK+AJ
      B(K3)=BK-AJ
      KK=K4+KSPAN
      IF(KK .LT. NN) GO TO 520
      KK=KK-NN
      IF(KK .LE. KSPAN) GO TO 520
      GO TO 700
C  TRANSFORM FOR ODD FACTORS
  600 K=NFAC(I)
      KSPNN=KSPAN
      KSPAN=KSPAN/K
      IF(K .EQ. 3) GO TO 320
      IF(K .EQ. 5) GO TO 510
      IF(K .EQ. JF) GO TO 640
      JF=K
      S1=RAD/FLOAT(K)
      C1=COS(S1)
      S1=SIN(S1)
      IF(JF .LE. MAXF)  GO TO 625
      CALL ERROFF
C/6S
C     CALL SETERR(27HF1FT - PRIME FACTOR .GT. 23,27,1,1)
C/7S
      CALL SETERR('F1FT - PRIME FACTOR .GT. 23',27,1,1)
C/
      RETURN
  625 CK(JF)=1.0
      SK(JF)=0.0
      J=1
  630 CK(J)=CK(K)*C1+SK(K)*S1
      SK(J)=CK(K)*S1-SK(K)*C1
      K=K-1
      CK(K)=CK(J)
      SK(K)=-SK(J)
      J=J+1
      IF(J .LT. K) GO TO 630
  640 K1=KK
      K2=KK+KSPNN
      AA=A(KK)
      BB=B(KK)
      AK=AA
      BK=BB
      J=1
      K1=K1+KSPAN
  650 K2=K2-KSPAN
      J=J+1
      AT(J)=A(K1)+A(K2)
      AK=AT(J)+AK
      BT(J)=B(K1)+B(K2)
      BK=BT(J)+BK
      J=J+1
      AT(J)=A(K1)-A(K2)
      BT(J)=B(K1)-B(K2)
      K1=K1+KSPAN
      IF(K1 .LT. K2) GO TO 650
      A(KK)=AK
      B(KK)=BK
      K1=KK
      K2=KK+KSPNN
      J=1
  660 K1=K1+KSPAN
      K2=K2-KSPAN
      JJ=J
      AK=AA
      BK=BB
      AJ=0.0
      BJ=0.0
      K=1
  670 K=K+1
      AK=AT(K)*CK(JJ)+AK
      BK=BT(K)*CK(JJ)+BK
      K=K+1
      AJ=AT(K)*SK(JJ)+AJ
      BJ=BT(K)*SK(JJ)+BJ
      JJ=JJ+J
      IF(JJ .GT. JF) JJ=JJ-JF
      IF(K .LT. JF) GO TO 670
      K=JF-J
      A(K1)=AK-BJ
      B(K1)=BK+AJ
      A(K2)=AK+BJ
      B(K2)=BK-AJ
      J=J+1
      IF(J .LT. K) GO TO 660
      KK=KK+KSPNN
      IF(KK .LE. NN) GO TO 640
      KK=KK-NN
      IF(KK .LE. KSPAN) GO TO 640
C  MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4)
  700 IF(I .EQ. M) GO TO 800
      KK=JC+1
  710 C2=1.0-CD
      S1=SD
  720 C1=C2
      S2=S1
      KK=KK+KSPAN
  730 AK=A(KK)
      A(KK)=C2*AK-S2*B(KK)
      B(KK)=S2*AK+C2*B(KK)
      KK=KK+KSPNN
      IF(KK .LE. NT) GO TO 730
      AK=S1*S2
      S2=S1*C2+C1*S2
      C2=C1*C2-AK
      KK=KK-NT+KSPAN
      IF(KK .LE. KSPNN) GO TO 730
      C2=C1-(CD*C1+SD*S1)
      S1=S1+(SD*C1-CD*S1)
C  THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C  ERROR. IF ROUNDED ARITHMETIC IS USED, THEY MAY BE DELETED.
      C1=0.5/(C2**2+S1**2)+0.5
      S1=C1*S1
      C2=C1*C2
      KK=KK-KSPNN+JC
      IF(KK .LE. KSPAN) GO TO 720
      KK=KK-KSPAN+JC+INC
      IF(KK .LE. JC+JC) GO TO 710
      GO TO 100
C  PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES
C  PERMUTATION FOR SQUARE FACTORS OF N
  800 NP(1)=KS
      IF(KT .EQ. 0) GO TO 890
      K=KT+KT+1
      IF(M .LT. K) K=K-1
      J=1
      NP(K+1)=JC
  810 NP(J+1)=NP(J)/NFAC(J)
      NP(K)=NP(K+1)*NFAC(J)
      J=J+1
      K=K-1
      IF(J .LT. K) GO TO 810
      K3=NP(K+1)
      KSPAN=NP(2)
      KK=JC+1
      K2=KSPAN+1
      J=1
      IF(N .NE. NTOT) GO TO 850
C  PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE)
  820 AK=A(KK)
      A(KK)=A(K2)
      A(K2)=AK
      BK=B(KK)
      B(KK)=B(K2)
      B(K2)=BK
      KK=KK+INC
      K2=KSPAN+K2
      IF(K2 .LT. KS) GO TO 820
  830 K2=K2-NP(J)
      J=J+1
      K2=NP(J+1)+K2
      IF(K2 .GT. NP(J)) GO TO 830
      J=1
  840 IF(KK .LT. K2) GO TO 820
      KK=KK+INC
      K2=KSPAN+K2
      IF(K2 .LT. KS) GO TO 840
      IF(KK .LT. KS) GO TO 830
      JC=K3
      GO TO 890
C  PERMUTATION FOR MULTIVARIATE TRANSFORM
  850 K=KK+JC
  860 AK=A(KK)
      A(KK)=A(K2)
      A(K2)=AK
      BK=B(KK)
      B(KK)=B(K2)
      B(K2)=BK
      KK=KK+INC
      K2=K2+INC
      IF(KK .LT. K) GO TO 860
      KK=KK+KS-JC
      K2=K2+KS-JC
      IF(KK .LT. NT) GO TO 850
      K2=K2-NT+KSPAN
      KK=KK-NT+JC
      IF(K2 .LT. KS) GO TO 850
  870 K2=K2-NP(J)
      J=J+1
      K2=NP(J+1)+K2
      IF(K2 .GT. NP(J)) GO TO 870
      J=1
  880 IF(KK .LT. K2) GO TO 850
      KK=KK+JC
      K2=KSPAN+K2
      IF(K2 .LT. KS) GO TO 880
      IF(KK .LT. KS) GO TO 870
      JC=K3
  890 IF(2*KT+1 .GE. M) RETURN
      KSPNN=NP(KT+1)
C  PERMUTATION FOR SQUARE-FREE FACTORS OF N
      J=M-KT
      NFAC(J+1)=1
  900 NFAC(J)=NFAC(J)*NFAC(J+1)
      J=J-1
      IF(J .NE. KT) GO TO 900
      KT=KT+1
      NN=NFAC(KT)-1
      IF(NN .LE. MAXP)  GO TO 901
      CALL ERROFF
C/6S
C     CALL SETERR(42HF1FT - SQUARE-FREE FACTOR PRODUCT .GT. 210,42,2,1)
C/7S
      CALL SETERR('F1FT - SQUARE-FREE FACTOR PRODUCT .GT. 210',42,2,1)
C/
      RETURN
  901 JJ=0
      J=0
      GO TO 906
  902 JJ=JJ-K2
      K2=KK
      K=K+1
      KK=NFAC(K)
  904 JJ=KK+JJ
      IF(JJ .GE. K2) GO TO 902
      NP(J)=JJ
  906 K2=NFAC(KT)
      K=KT+1
      KK=NFAC(K)
      J=J+1
      IF(J .LE. NN) GO TO 904
C  DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1
      J=0
      GO TO 914
 910  K=KK
      KK=NP(K)
      NP(K)=-KK
      IF(KK .NE. J) GO TO 910
      K3=KK
  914 J=J+1
      KK=NP(J)
      IF(KK .LT. 0) GO TO 914
      IF(KK .NE. J) GO TO 910
      NP(J)=-J
      IF(J .NE. NN) GO TO 914
      MAXF=INC*MAXF
C  REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES
      GO TO 950
  924 J=J-1
      IF(NP(J) .LT. 0) GO TO 924
      JJ=JC
  926 KSPAN=JJ
      IF(JJ .GT. MAXF) KSPAN=MAXF
      JJ=JJ-KSPAN
      K=NP(J)
      KK=JC*K+II+JJ
      K1=KK+KSPAN
      K2=0
  928 K2=K2+1
      AT(K2)=A(K1)
      BT(K2)=B(K1)
      K1=K1-INC
      IF(K1 .NE. KK) GO TO 928
  932 K1=KK+KSPAN
      K2=K1-JC*(K+NP(K))
      K=-NP(K)
  936 A(K1)=A(K2)
      B(K1)=B(K2)
      K1=K1-INC
      K2=K2-INC
      IF(K1 .NE. KK) GO TO 936
      KK=K2
      IF(K .NE. J) GO TO 932
      K1=KK+KSPAN
      K2=0
  940 K2=K2+1
      A(K1)=AT(K2)
      B(K1)=BT(K2)
      K1=K1-INC
      IF(K1 .NE. KK) GO TO 940
      IF(JJ .NE. 0) GO TO 926
      IF(J .NE. 1) GO TO 924
  950 J=K3+1
      NT=NT-KSPNN
      II=NT-INC+1
      IF(NT .GE. 0) GO TO 924
      RETURN
       END

.