**-------------------------------------------**                         
CCC             MIGRATION PROGRAM                                       
**-------------------------------------------**                         
      PROGRAM MAIN                                                      
      COMMON/GENRAL/NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0, 
     .              NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT   
      COMMON/FFT/TRIGX(1000),TRIGY(1000),IFAX(15),IFAY(15),RKX(400),    
     .           RKY(400),RKXC(800),RKYC(800)                           
      COMMON/IO/NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES                
      COMMON/PP1/PP1( 61000),PM1( 60000)                                
      COMMON/CDPP1/CDPP1( 60000),CDPM1( 60000)                          
      COMMON BUF( 120320),AW( 61000)                                    
      COMMON/VEL/VEL( 15000),DSECT( 15000)                              
      COMMON/FF1/RTRIG(1500),CTRIG(1500),IFXA(15),IFYA(15)              
C MAXIMUM RECORD LENGTH FOR DIRECT ACCESS FILES IN BYTES. HERE, NTT*8   
C IT IS ASSUMEM MAXRL.LE.962560                                         
      PARAMETER (MAXRL=962560)                                          
C                                                                       
      REAL T1,T2,T3,T4                                                  
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
C                                                                       
      OPEN(UNIT=41,FILE='SMI41',STATUS='OLD')                           
      OPEN(UNIT=10,FILE='SMV',STATUS='NEW')                             
                                                                        
C-------                                                                
C  READ COMMON BLOCKS                                                   
C-------                                                                
C                                                                       
C.....CALL TIMING ROUTINES                                              
C                                                                       
      CALL CPUTIM(T1)                                                   
      CALL ELAPSE(T2)                                                   
C                                                                       
      READ(41,201)NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0    
      READ(41,202)NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT     
      READ(41,203)(TRIGX(I),TRIGY(I),I=1,1000)                          
      READ(41,204)(IFAX(I),IFAY(I),I=1,15)                              
      READ(41,203)(RKX(I),RKY(I),I=1,400)                               
      READ(41,203)(RKXC(I),RKYC(I),I=1,800)                             
      READ(41,205)NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES              
201   FORMAT(3(I5,1X),3(E16.8,1X)/I5,1X,2(E16.8,1X),I10,1X,E16.8)       
202   FORMAT(3(I5,1X),E16.8,1X,2(I5,1X)/2(E16.8,1X),I5)                 
203   FORMAT(4(E16.8,1X)/)                                              
204   FORMAT(5(I10,1X)/)                                                
205   FORMAT(4(I10,1X)/)                                                
***                                                                     
* -------  OPEN FILES                                                   
***                                                                     
C... NOTE:  ON CRAYS, FOR DIRECT ACCESS FILES THE RECORD LENGTH         
C... 'RECL' IS IN BYTES                                                 
                                                                        
C     OPEN(ISSD,FILE='SMT1',FORM='UNFORMATTED',                         
C    .     ACCESS='DIRECT',RECL=8*NTT)                                  
C REPLACE L. POINTER 10-89 TO HANDLE DIFFERENT RECORD LENGTH if NESS.   
C MAXRL=962560=8*NTT                                                    
      OPEN(UNIT=ISSD,FILE='SMT1',FORM='UNFORMATTED',                    
     .     ACCESS='DIRECT',RECL=MAXRL,STATUS='NEW')                     
      OPEN(UNIT=IDISK1,FILE='SMT2',FORM='UNFORMATTED',STATUS='NEW')     
      OPEN(UNIT=IDISK2,FILE='SMO6',STATUS='NEW')                        
C                                                                       
      CALL FFTFAX(NY,IFYA,RTRIG)                                        
      CALL CFTFAX(NX,IFXA,CTRIG)                                        
      NXY=NX*NY                                                         
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
      DO 77 I=1,NSTEPS                                                  
       DO 78 K=1,NXY                                                    
78      VEL(K)=2000. + FLOAT(K)                                         
77     WRITE(IDISK1)(VEL(K),K=1,NXY)                                    
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
      REWIND(IDISK1)                                                    
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
      DO 80 I=1,NTT                                                     
80     BUF(I)=1.                                                        
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
      DO 90 I=1,IWCUT,NWJUMP                                            
       CALL WSSD(BUF,I)                                                 
90     CONTINUE                                                         
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
C     REWIND(ISSD)                                                      
                                                                        
C     CALL REMARK('START MIGRAT')                                       
      CALL MIGRAT(BUF,DSECT,VEL,PP1,PM1,CDPP1,CDPM1,AW)                 
C                                                                       
      CALL CPUTIM(T3)                                                   
      CALL ELAPSE(T4)                                                   
C                                                                       
      WRITE(10,*)                                                       
      WRITE(10,*) 'ELAPSED CPU TIME IN SECONDS: ',T3-T1                 
C MFLOPS = MFLOP 1 CPU CRAY X-MP/CPU TIME                               
      WRITE(10,*) 'MFLOPS RATE: ', 10558.609/(T3-T1)                    
      WRITE(10,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',                
     *            MOD(T4-T2+1000000.0,1000000.0)                        
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
C                                                                       
C     CALL REMARK('END MIGRAT')                                         
      STOP                                                              
      END                                                               
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$                             
      SUBROUTINE MIGRAT(BUF,DSECT,VEL,PP1,PM1,CDPP1,CDPM1,AW)           
      COMMON/GENRAL/NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0, 
     .              NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT   
      COMMON/FFT/TRIGX(1000),TRIGY(1000),IFAX(15),IFAY(15),RKX(400),    
     .           RKY(400),RKXC(800),RKYC(800)                           
      COMMON/IO/NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES                
      DIMENSION BUF(NX,NY,2,1),DSECT(NX,1),VEL(1)                       
      DIMENSION PP1(NX,NY+2,1),PM1(NX,NY,1),CDPP1(NX,NY,1),             
     .          CDPM1(NX,NY,1),AW(NX,NY+2,1)                            
CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (3)                                                  
CCCC                                                                    
      KK=1                                                              
      NXY=NX*NY                                                         
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
                                                                        
      DO 100 NZ=1,NSTEPS,NZJUMP                                         
C                                                                       
C READ VELOCITY FOR DEPTH LEVEL                                         
      READ(IDISK1)(VEL(I),I=1,NXY)                                      
      VELMIN=10000000.                                                  
      VELMAX=0.                                                         
      DO 30 I=1,NXY                                                     
      DSECT(I,1)=0.                                                     
      VELMAX=MAX(VEL(I),VELMAX)                                         
30    VELMIN=MIN(VEL(I),VELMIN)                                         
      DO 200 NW=1,IWCUT,NWJUMP                                          
      NN=MIN(NWJUMP,IWCUT-NW+1)                                         
C READ P & CDP FROM SSD                                                 
      CALL RSSD(BUF,NW)                                                 
*--                                                                     
      DO  20 IW=1,NN                                                    
      W=(NW+IW-1)*DW                                                    
*--                                                                     
      IF(NZ.EQ.9)THEN                                                   
                                                                        
      WC=W*W/(VEL0*VEL0)                                                
         DO 112 J=1,NY                                                  
         DO 112 I=1,NX                                                  
         PP1(I,J,KK)=BUF(I,J,1,IW)                                      
112      PM1(I,J,KK)=BUF(I,J,2,IW)                                      
      CALL FFF(PP1(1,1,KK),AW(1,1,KK),NX,NY)                            
      CALL FFF(PM1(1,1,KK),PP1(1,1,KK),NX,NY)                           
CC                                                                      
      ZERO=1.E-40                                                       
      IJ=0                                                              
      DO 8 J=1,NY+2,2                                                   
      IJ=IJ+1                                                           
      DO 8 I=1,NX                                                       
      FILT=1.                                                           
      TEMP=ABS(WC-RKX(I)-RKY(IJ))                                       
      PP1(I,J,KK)=-SQRT(TEMP)*PP1(I,J,KK)*VEL0                          
      PP1(I,J+1,KK)=-SQRT(TEMP)*PP1(I,J+1,KK)*VEL0                      
      FF=WC-RKX(I)-RKY(IJ)                                              
C      FILT=CVMGT(0.,FILT,FF.LT.ZERO)                                   
      IF(FF.LT.ZERO) FILT = 0.                                          
      AW(I,J,KK)=AW(I,J,KK)*FILT                                        
      PP1(I,J,KK)=PP1(I,J,KK)*FILT                                      
      AW(I,J+1,KK)=AW(I,J+1,KK)*FILT                                    
      PP1(I,J+1,KK)=PP1(I,J+1,KK)*FILT                                  
8     CONTINUE                                                          
C                                                                       
CC                                                                      
      CALL FFB(PM1(1,1,KK),PP1(1,1,KK),NX,NY)                           
      CALL FFB(PP1(1,1,KK),AW(1,1,KK),NX,NY)                            
CC                                                                      
         DO 113 J=1,NY                                                  
         DO 113 I=1,NX                                                  
         BUF(I,J,1,IW)=PP1(I,J,KK)                                      
113      BUF(I,J,2,IW)=PM1(I,J,KK)                                      
                                                                        
      ENDIF                                                             
                                                                        
*********************                                                   
      CALL MARCH(BUF,VEL,PP1(1,1,KK),PM1(1,1,KK),CDPP1(1,1,KK),         
     .           CDPM1(1,1,KK),AW(1,1,KK),NW,IW)                        
*********************                                                   
        DO 55 J=1,NY                                                    
        DO 55 I=1,NX                                                    
        PP1(I,J,KK)=BUF(I,J,1,IW)                                       
55      PM1(I,J,KK)=BUF(I,J,2,IW)                                       
                                                                        
      WC=W*W/(VELMAX*VELMAX)                                            
      CALL FFF(PP1(1,1,KK),AW(1,1,KK),NX,NY)                            
      CALL FFF(PM1(1,1,KK),PP1(1,1,KK),NX,NY)                           
CC                                                                      
      ZERO=1.E-40                                                       
      IJ=0                                                              
      DO 9 J=1,NY+2,2                                                   
      IJ=IJ+1                                                           
      DO 9 I=1,NX                                                       
      FILT=1.                                                           
      FF=WC-RKX(I)-RKY(IJ)                                              
C      FILT=CVMGT(0.,FILT,FF.LT.ZERO)                                   
      IF(FF.LT.ZERO) FILT = 0.                                          
      AW(I,J,KK)=AW(I,J,KK)*FILT                                        
      PP1(I,J,KK)=PP1(I,J,KK)*FILT                                      
      AW(I,J+1,KK)=AW(I,J+1,KK)*FILT                                    
      PP1(I,J+1,KK)=PP1(I,J+1,KK)*FILT                                  
9     CONTINUE                                                          
CC                                                                      
      CALL FFB(PM1(1,1,KK),PP1(1,1,KK),NX,NY)                           
      CALL FFB(PP1(1,1,KK),AW(1,1,KK),NX,NY)                            
                                                                        
         DO 66 J=1,NY                                                   
         DO 66 I=1,NX                                                   
         BUF(I,J,1,IW)=PP1(I,J,KK)                                      
66       BUF(I,J,2,IW)=PM1(I,J,KK)                                      
      DO 21 J=1,NY                                                      
      DO 21 I=1,NX                                                      
21    DSECT(I,J)=DSECT(I,J)+BUF(I,J,1,IW)                               
*--                                                                     
20    CONTINUE                                                          
C WRITE P & CDP TO SSD                                                  
      CALL WSSD(BUF,NW)                                                 
200   CONTINUE                                                          
                                                                        
100   CONTINUE                                                          
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
                                                                        
C... CALL THE VERIFICATION ROUTINE                                      
                                                                        
      CALL VERIFY(IDISK2,NX,DSECT)                                      
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
                                                                        
      RETURN                                                            
      END                                                               
**********************************                                      
      SUBROUTINE MARCH(BUF,VEL,PP1,PM1,CDPP1,CDPM1,AW,NW,IW)            
      COMMON/GENRAL/NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0, 
     .              NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT   
      COMMON/FFT/TRIGX(1000),TRIGY(1000),IFAX(15),IFAY(15),RKX(400),    
     .           RKY(400),RKXC(800),RKYC(800)                           
      COMMON/IO/NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES                
      DIMENSION BUF(NX,NY,2,1),VEL(NX,1),BES(150)                       
      DIMENSION PP1(NX,1),PM1(NX,1),CDPP1(NX,1),CDPM1(NX,1),AW(NX,1)    
      DATA BES/150*1.0/                                                 
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
      W=(NW+IW-1)*DW                                                    
      II=NW+IW-1                                                        
      DZR = 1.                                                          
      CALL BESSEL(NBES,II)                                              
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
C                                                                       
      DO 30 J=1,NY                                                      
      DO 30 I=1,NX                                                      
      PP1(I,J)=0.                                                       
      CDPP1(I,J)=0.                                                     
      PM1(I,J)=BUF(I,J,1,IW)                                            
30    CDPM1(I,J)=BUF(I,J,2,IW)                                          
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
CC                                                                      
      CALL MUL(PM1,CDPM1,PP1,CDPP1,VEL,DZR,W,AW)                        
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
CC                                                                      
      DO 40 J=1,NY                                                      
      DO 40 I=1,NX                                                      
      BUF(I,J,1,IW)=BES(1)*PM1(I,J)+BES(2)*PP1(I,J)                     
40    BUF(I,J,2,IW)=BES(1)*CDPM1(I,J)+BES(2)*CDPP1(I,J)                 
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
**                                                                      
      DZR=DZR*2.                                                        
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
**                                                                      
      DO 500 NB=3,NBES                                                  
CC                                                                      
      CALL MUL(PP1,CDPP1,PM1,CDPM1,VEL,DZR,W,AW)                        
CC                                                                      
      DO 60 J=1,NY                                                      
      DO 60 I=1,NX                                                      
      BUF(I,J,1,IW)=BUF(I,J,1,IW)+BES(NB)*PM1(I,J)                      
      BUF(I,J,2,IW)=BUF(I,J,2,IW)+BES(NB)*CDPM1(I,J)                    
      CC=PM1(I,J)                                                       
      CCC=CDPM1(I,J)                                                    
      PM1(I,J)=PP1(I,J)                                                 
      CDPM1(I,J)=CDPP1(I,J)                                             
      PP1(I,J)=CC                                                       
      CDPP1(I,J)=CCC                                                    
60    CONTINUE                                                          
500   CONTINUE                                                          
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
CC                                                                      
      RETURN                                                            
      END                                                               
*****************************                                           
      SUBROUTINE MUL(PIN,CDPIN,POUT,CDPOUT,VEL,DZR,W,AW)                
      COMMON/GENRAL/NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0, 
     .              NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT   
      COMMON/FFT/TRIGX(1000),TRIGY(1000),IFAX(15),IFAY(15),RKX(400),    
     .           RKY(400),RKXC(800),RKYC(800)                           
      COMMON/IO/NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES                
      DIMENSION PIN(NX,1),CDPIN(NX,1),POUT(NX,1),CDPOUT(NX,1),VEL(NX,1) 
      DIMENSION AW(1)                                                   
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
      W2=W*W                                                            
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
C                                                                       
      DO 1 J=1,NY                                                       
      DO 1 I=1,NX                                                       
      POUT(I,J)=POUT(I,J)+DZR*CDPIN(I,J)/VEL(I,J)                       
1     CDPOUT(I,J)=CDPOUT(I,J)-DZR*PIN(I,J)*W2/VEL(I,J)                  
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
C                                                                       
      CALL DIFX(PIN,CDPOUT,AW,AW,NX,NY,VEL,DZR)                         
C                                                                       
      CALL DIFY(PIN,CDPOUT,AW,AW,NX,NY,VEL,DZR)                         
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
C--------------------------------------------------------------         
C                                                                       
      SUBROUTINE DIFX(VIN,VOUT,AEQ,A,NNXT,NNYT,VEL,DZR)                 
      COMMON/FFT/TRIGX(1000),TRIGY(1000),IFAX(15),IFAY(15),RKX(400),    
     .           RKY(400),RKXC(800),RKYC(800)                           
      DIMENSION VIN(NNXT,1),VOUT(NNXT,1),VEL(NNXT,1)                    
      DIMENSION AEQ(NNXT,1),A(1)                                        
      COMMON/WORK/WORK(20000)                                           
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
      JC=(NNYT-1)/2+1                                                   
      JC2=JC*2                                                          
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
      DO 1 I=1,NNXT*(NNYT+2)                                            
1     A(I)=0.                                                           
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
      DO 2 NY=1,NNYT                                                    
      DO 2 NX=1,NNXT                                                    
2     AEQ(NX,NY)=VIN(NX,NY)                                             
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
****                                                                    
      DO 600 IS=1,JC,64                                                 
      LN=MIN(64,JC-IS+1)                                                
      IP=2*NNXT*(IS-1)+1                                                
      CALL CFFT9X(A(IP),A(IP+NNXT),WORK,TRIGX,IFAX,1,2*NNXT,NNXT,LN,-1) 
600   CONTINUE                                                          
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
****                                                                    
      DO 106 J=1,JC2,2                                                  
      IND=1                                                             
      DO 105 I=1,NNXT                                                   
      TR=AEQ(I,J)*RKXC(IND)-AEQ(I,J+1)*RKXC(IND+1)                      
      AEQ(I,J+1)=AEQ(I,J)*RKXC(IND+1)+AEQ(I,J+1)*RKXC(IND)              
      AEQ(I,J)=TR                                                       
105   IND=IND+2                                                         
106   CONTINUE                                                          
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
****                                                                    
      DO 700 IS=1,JC,64                                                 
      LN=MIN(64,JC-IS+1)                                                
      IP=2*NNXT*(IS-1)+1                                                
      CALL CFFT9X(A(IP),A(IP+NNXT),WORK,TRIGX,IFAX,1,2*NNXT,NNXT,LN,1)  
700   CONTINUE                                                          
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
****                                                                    
         DO 8 NY=1,NNYT                                                 
         DO 8 NX=1,NNXT                                                 
8        VOUT(NX,NY)=VOUT(NX,NY)-AEQ(NX,NY)*VEL(NX,NY)*DZR              
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      END                                                               
C--------------------------------------------------------------         
C                                                                       
      SUBROUTINE DIFY(VIN,VOUT,AEQ,A,NNXT,NNYT,VEL,DZR)                 
      COMMON/FFT/TRIGX(1000),TRIGY(1000),IFAX(15),IFAY(15),RKX(400),    
     .           RKY(400),RKXC(800),RKYC(800)                           
      DIMENSION VIN(NNXT,1),VOUT(NNXT,1),VEL(NNXT,1)                    
      DIMENSION AEQ(2*((NNXT-1)/2+1),1),A(1)                            
      COMMON/WORK/WORK(20000)                                           
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
      JC=(NNXT-1)/2+1                                                   
      JC2=JC*2                                                          
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
      DO 1 I=1,NNXT*(NNYT+2)                                            
1     A(I)=0.                                                           
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
      DO 2 NY=1,NNYT                                                    
      DO 2 NX=1,NNXT                                                    
2     AEQ(NX,NY)=VIN(NX,NY)                                             
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
****                                                                    
      DO 600 IS=1,JC,64                                                 
      LN=MIN(64,JC-IS+1)                                                
      IP=(IS-1)*2+1                                                     
      CALL CFFT9X(A(IP),A(IP+1),WORK,TRIGY,IFAY,JC2,2,NNYT,LN,-1)       
600   CONTINUE                                                          
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
****                                                                    
      DO 106 I=1,JC2,2                                                  
      IND=1                                                             
      DO 105 J=1,NNYT                                                   
      TR=AEQ(I,J)*RKYC(IND)-AEQ(I+1,J)*RKYC(IND+1)                      
      AEQ(I+1,J)=AEQ(I,J)*RKYC(IND+1)+AEQ(I+1,J)*RKYC(IND)              
      AEQ(I,J)=TR                                                       
105   IND=IND+2                                                         
106   CONTINUE                                                          
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
****                                                                    
      DO 700 IS=1,JC,64                                                 
      LN=MIN(64,JC-IS+1)                                                
      IP=(IS-1)*2+1                                                     
      CALL CFFT9X(A(IP),A(IP+1),WORK,TRIGY,IFAY,JC2,2,NNYT,LN,1)        
700   CONTINUE                                                          
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
****                                                                    
         DO 12 NY=1,NNYT                                                
         DO 12 NX=1,NNXT                                                
12       VOUT(NX,NY)=VOUT(NX,NY)-AEQ(NX,NY)*VEL(NX,NY)*DZR              
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
      RETURN                                                            
      END                                                               
*****************************                                           
       SUBROUTINE RSSD(BUF,NW)                                          
      COMMON/GENRAL/NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0, 
     .              NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT   
      COMMON/IO/NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES                
      DIMENSION BUF(1)                                                  
      PARAMETER (MAXRL = 962560)                                        
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
C--                                                                     
      IREC=(NW-1)/NWJUMP+1                                              
C MODIFIED L. POINTER CSRD 10-89 TO HANDLE SMALLER RECORD LENGTH.       
C CALCULATE NUMBER OF FULL MAXRL BYTE RECORDS                           
      NRECS = (NTT*8)/MAXRL                                             
C CALCULATE HOW MANY BYTES LEFT                                         
      NLEFT = MOD(NTT*8,MAXRL)                                          
      IF (NLEFT.NE.0) THEN                                              
         IREC=(IREC-1)*(NRECS+1) + 1                                    
      ELSE                                                              
         IREC=(IREC-1)*NRECS + 1                                        
      ENDIF                                                             
      CALL GETWAV(ISSD,BUF,IREC,NTT,IER,NRECS,NLEFT)                    
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
                                                                        
       RETURN                                                           
       END                                                              
C----------------------------------------------C                        
      SUBROUTINE WSSD(BUF,NW)                                           
      COMMON/GENRAL/NTRACE,NSAMP,NLINE,DX,DY,DZ,NFFT,DT,FCUT,NTOT,VEL0, 
     .              NWJUMP,NZJUMP,NSTEPS,DW,NX,NY,VELMIN,VELMAX,IWCUT   
      COMMON/IO/NTT,ISSD,IDISK1,IDISK2,IIN,NIO,NLIN,NRES                
      DIMENSION BUF(1)                                                  
      PARAMETER (MAXRL =962560)                                         
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
                                                                        
      IREC=(NW-1)/NWJUMP+1                                              
C MODIFIED L. POINTER CSRD 10-89 TO HANDLE SMALLER RECORD LENGTH.       
C CALCULATE NUMBER OF FULL MAXRL BYTE RECORDS                           
      NRECS = (NTT*8)/MAXRL                                             
C CALCULATE HOW MANY BYTES LEFT                                         
      NLEFT = MOD(NTT*8,MAXRL)                                          
      IF (NLEFT.NE.0) THEN                                              
         IREC=(IREC-1)*(NRECS+1) + 1                                    
      ELSE                                                              
         IREC=(IREC-1)*NRECS + 1                                        
      ENDIF                                                             
      CALL PUTWAV(ISSD,BUF,IREC,NTT,IER,NRECS,NLEFT)                    
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
                                                                        
       RETURN                                                           
       END                                                              
****************************************************                    
      SUBROUTINE BESSEL(NB,I)                                           
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
      IF(I.LE.6)NB=4                                                    
      IF(I.GT.6.AND.I.LE.26)NB=5                                        
      IF(I.GT.26.AND.I.LE.64)NB=6                                       
      IF(I.GT.64.AND.I.LE.120)NB=7                                      
      IF(I.GT.120.AND.I.LE.189)NB=8                                     
      IF(I.GT.189.AND.I.LE.281)NB=9                                     
      IF(I.GT.281)NB=10                                                 
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      RETURN                                                            
      END                                                               
C     #######################################                           
      SUBROUTINE FFF(VIN,A,NX,NY)                                       
      COMMON/FF1/RTRIG(1500),CTRIG(1500),IFXA(15),IFYA(15)              
      COMMON/WORK/WORK(20000)                                           
      DIMENSION VIN(NX,1),A(1)                                          
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
      JC=(NY-1)/2+1                                                     
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
       DO 2 J=1,NX*(NY+2)                                               
2      A(J)=0.                                                          
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
       DO 1 J=1,NX*NY                                                   
1      A(J)=VIN(J,1)                                                    
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
      CALL FFT991(A,WORK,RTRIG,IFYA,NX,1,NY,NX,-1)                      
C                                                                       
      CALL CFFT9X(A(1),A(1+NX),WORK,CTRIG,IFXA,1,2*NX,NX,JC+1,-1)       
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C     #######################################                           
      SUBROUTINE FFB(VIN,A,NX,NY)                                       
      COMMON/FF1/RTRIG(1500),CTRIG(1500),IFXA(15),IFYA(15)              
      DIMENSION VIN(NX,1),A(1)                                          
      COMMON/WORK/WORK(20000)                                           
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
      JC=(NY-1)/2+1                                                     
C                                                                       
      CALL CFFT9X(A(1),A(1+NX),WORK,CTRIG,IFXA,1,2*NX,NX,JC+1,1)        
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
C                                                                       
       DO 9 J=1,NX*(NY+2)                                               
9      A(J)=A(J)/NX                                                     
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
C                                                                       
      CALL FFT991(A,WORK,RTRIG,IFYA,NX,1,NY,NX,1)                       
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
       DO 5 J=1,NX*NY                                                   
5      VIN(J,1)=A(J)                                                    
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      RETURN                                                            
      END                                                               
**********************************************************************  
**    UTIL                                                              
*****************************************************                   
      SUBROUTINE CFFT9X(AR,AI,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)     
C                                                                       
C PURPOSE      PERFORMS MULTIPLE FAST FOURIER TRANSFORMS.  THIS PACKAGE 
C              WILL PERFORM A NUMBER OF SIMULTANEOUS COMPLEX PERIODIC   
C              FOURIER TRANSFORMS OR CORRESPONDING INVERSE TRANSFORMS.  
C              THAT IS, GIVEN A SET OF COMPLEX GRIDPOINT VECTORS, THE   
C              PACKAGE RETURNS A SET OF COMPLEX FOURIER                 
C              COEFFICIENT VECTORS, OR VICE VERSA.  THE LENGTH OF THE   
C              TRANSFORMS MUST BE A NUMBER GREATER THAN 1 THAT HAS      
C              NO PRIME FACTORS OTHER THAN 2, 3, AND 5.                 
C                                                                       
C              THE PACKAGE CFFT9X CONTAINS SEVERAL USER-LEVEL ROUTINES: 
C                                                                       
C            SUBROUTINE CFTFAX                                          
C                AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE     
C                BEFORE A SEQUENCE OF CALLS TO CFFT9X                   
C                (PROVIDED THAT N IS NOT CHANGED).                      
C                                                                       
C            SUBROUTINE CFFT9X                                          
C                THE ACTUAL TRANSFORM ROUTINE ROUTINE, CABABLE OF       
C                PERFORMING BOTH THE TRANSFORM AND ITS INVERSE.         
C                HOWEVER, AS THE TRANSFORMS ARE NOT NORMALIZED,         
C                THE APPLICATION OF A TRANSFORM FOLLOWED BY ITS         
C                INVERSE WILL YIELD THE ORIGINAL VALUES MULTIPLIED      
C                                                                       
C USAGE        LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 0, 
C              Q .GE. 0, AND R .GE. 0.  THEN A TYPICAL SEQUENCE OF      
C              CALLS TO TRANSFORM A GIVEN SET OF COMPLEX VECTORS OF     
C              LENGTH N TO A SET OF (UNSCALED) COMPLEX FOURIER          
C              COEFFICIENT VECTORS OF LENGTH N IS                       
C                                                                       
C                   DIMENSION IFAX(13),TRIGS(2*N)                       
C                   COMPLEX A(...), WORK(...)                           
C                                                                       
C                   CALL CFTFAX (N, IFAX, TRIGS)                        
C                   CALL CFFT9X (A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C                                                                       
C              THE OUTPUT VECTORS OVERWRITE THE INPUT VECTORS, AND      
C              THESE ARE STORED IN A.  WITH APPROPRIATE CHOICES FOR     
C              THE OTHER ARGUMENTS, THESE VECTORS MAY BE CONSIDERED     
C              EITHER THE ROWS OR THE COLUMNS OF THE ARRAY A.           
C              SEE THE INDIVIDUAL WRITE-UPS FOR CFTFAX AND              
C              CFFT9X BELOW, FOR A DETAILED DESCRIPTION OF THE          
C              ARGUMENTS.                                               
C                                                                       
C HISTORY      THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN   
C              NOVEMBER, 1978.  IT WAS MODIFIED, DOCUMENTED, AND TESTED 
C              FOR NCAR BY RUSS REW IN SEPTEMBER, 1980.  IT WAS         
C              FURTHER MODIFIED FOR THE FULLY COMPLEX CASE BY DAVE      
C              FULKER IN NOVEMBER, 1980.                                
C                                                                       
C-----------------------------------------------------------------------
      DIMENSION AR(1),AI(1),WORK(1),TRIGS(1),IFAX(1)                    
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
C                                                                       
C     SUBROUTINE "CFFT9X" - MULTIPLE FAST COMPLEX FOURIER TRANSFORM     
C                                                                       
C     A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA                   
C     WORK IS AN AREA OF SIZE N*LOT                                     
C     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES       
C     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N                
C     INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'                    
C         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)                    
C     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR       
C     N IS THE LENGTH OF THE DATA VECTORS                               
C     LOT IS THE NUMBER OF DATA VECTORS                                 
C     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT               
C           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL               
C                                                                       
C                                                                       
C     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN      
C     PARALLEL.                                                         
      LL=LOT+LOT                                                        
      INK=INC                                                           
      JUM=JUMP                                                          
      NFAX=IFAX(1)                                                      
C                                                                       
C                                                                       
C     PERFORM THE TRANSFORM PASSES, ONE PASS FOR EACH FACTOR.  DURING   
C     EACH PASS THE DATA ARE MOVED FROM A TO WORK OR FROM WORK TO A.    
C     IF NFAX IS ODD, THE LAST FACTOR IS PROCESSED IN PLACE.            
C                                                                       
      LA=1                                                              
      NMOD = MOD(NFAX,2)                                                
      IF(ISIGN.GT.0) THEN                                               
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
                                                                        
          DO 10 K=1+NMOD,NFAX,2                                         
          CALL CPASS (AR(1),AI(1),WORK(1),WORK(1+LOT),TRIGS,            
     *    INK,LL,JUM,1,LOT,N,IFAX(K+1),LA)                              
          LA=LA*IFAX(K+1)                                               
          CALL CPASS (WORK(1),WORK(1+LOT),AR(1),AI(1),TRIGS,            
     *    LL,INK,1,JUM,LOT,N,IFAX(K+2),LA)                              
          LA=LA*IFAX(K+2)                                               
   10     CONTINUE                                                      
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
          IF( NMOD .NE. 0 ) THEN                                        
CXXX NEW MOD                                                            
           IF (IFAX(2) .NE. 5) THEN                                     
              CALL CPASS (AR(1),AI(1),AR(1),AI(1),TRIGS,                
     *        INK,INK,JUM,JUM,LOT,N,IFAX(2),LA)                         
           ELSE                                                         
            I = 1                                                       
CX10551 CALL DOINX( 551 )                                               
CY10551 CALL DOINY( 551 )                                               
            DO 12 L=1,2*N*LOT,2                                         
             WORK(L) = AR(I)                                            
             WORK(L+1) = AI(I)                                          
   12        I = I + INK                                                
CX20552 CALL DOOUTX( 552 )                                              
CY20552 CALL DOOUTY( 552 )                                              
              CALL CPASS (WORK(1),WORK(2),AR(1),AI(1),TRIGS,            
     *        2,INK,JUM,JUM,LOT,N,IFAX(2),LA)                           
           ENDIF                                                        
          ENDIF                                                         
CXXX END NEW MOD                                                        
      ELSE                                                              
CX10553 CALL DOINX( 553 )                                               
CY10553 CALL DOINY( 553 )                                               
                                                                        
          DO 30 K=1+NMOD,NFAX,2                                         
          CALL CPASSM(AR(1),AI(1),WORK(1),WORK(1+LOT),TRIGS,            
     *    INK,LL,JUM,1,LOT,N,IFAX(K+1),LA)                              
          LA=LA*IFAX(K+1)                                               
          CALL CPASSM(WORK(1),WORK(1+LOT),AR(1),AI(1),TRIGS,            
     *    LL,INK,1,JUM,LOT,N,IFAX(K+2),LA)                              
          LA=LA*IFAX(K+2)                                               
   30     CONTINUE                                                      
CX20554 CALL DOOUTX( 554 )                                              
CY20554 CALL DOOUTY( 554 )                                              
          IF( NMOD .NE. 0 ) THEN                                        
CXXX NEW MOD                                                            
           IF (IFAX(2) .NE. 5) THEN                                     
              CALL CPASSM(AR(1),AI(1),AR(1),AI(1),TRIGS,                
     *        INK,INK,JUM,JUM,LOT,N,IFAX(2),LA)                         
           ELSE                                                         
            I = 1                                                       
CX10555 CALL DOINX( 555 )                                               
CY10555 CALL DOINY( 555 )                                               
            DO 32 L=1,2*N*LOT,2                                         
             WORK(L) = AR(I)                                            
             WORK(L+1) = AI(I)                                          
   32        I = I + INK                                                
CX20556 CALL DOOUTX( 556 )                                              
CY20556 CALL DOOUTY( 556 )                                              
              CALL CPASSM(WORK(1),WORK(2),AR(1),AI(1),TRIGS,            
     *        2,INK,JUM,JUM,LOT,N,IFAX(2),LA)                           
           ENDIF                                                        
CXXX END NEW MOD                                                        
          ENDIF                                                         
      ENDIF                                                             
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE CPASS (A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA)
      DIMENSION A(N),B(N),C(N),D(N),TRIGS(N)                            
C                                                                       
C     SUBROUTINE "CPASS" - MULTIPLE VERSION OF "CPASSA"                 
C     PERFORMS ONE PASS THROUGH DATA                                    
C     AS PART OF MULTIPLE COMPLEX (INVERSE) FFT ROUTINE                 
C     A IS FIRST REAL INPUT VECTOR                                      
C     B IS FIRST IMAGINARY INPUT VECTOR                                 
C     C IS FIRST REAL OUTPUT VECTOR                                     
C     D IS FIRST IMAGINARY OUTPUT VECTOR                                
C     TRIGS IS PRECALCULATED TABLE OF SINES & COSINES                   
C     INC1 IS ADDRESSING INCREMENT FOR A AND B                          
C     INC2 IS ADDRESSING INCREMENT FOR C AND D                          
C     INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S                    
C     INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S                    
C     LOT IS THE NUMBER OF VECTORS                                      
C     N IS LENGTH OF VECTORS                                            
C     IFAC IS CURRENT FACTOR OF N                                       
C     LA IS PRODUCT OF PREVIOUS FACTORS                                 
C                                                                       
      DATA SIN36/0.587785252292473/,COS36/0.809016994374947/,           
     *     SIN72/0.951056516295154/,COS72/0.309016994374947/,           
     *     SIN60/0.866025403784437/                                     
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
C                                                                       
      M=N/IFAC                                                          
      IINK=M*INC1                                                       
      JINK=LA*INC2                                                      
      JUMP=(IFAC-1)*JINK                                                
      IBASE=0                                                           
      JBASE=0                                                           
      IGO=IFAC-1                                                        
      IF (IGO.GT.4) THEN                                                
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      ENDIF                                                             
      GO TO (10,50,90,130),IGO                                          
C                                                                       
C     CODING FOR FACTOR 2                                               
C                                                                       
   10 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
CX10557 CALL DOINX( 557 )                                               
CY10557 CALL DOINY( 557 )                                               
      DO 20 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 15 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+A(IB+I)                                           
      D(JA+J)=B(IA+I)+B(IB+I)                                           
      C(JB+J)=A(IA+I)-A(IB+I)                                           
      D(JB+J)=B(IA+I)-B(IB+I)                                           
      I=I+INC3                                                          
      J=J+INC4                                                          
   15 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   20 CONTINUE                                                          
CX20558 CALL DOOUTX( 558 )                                              
CY20558 CALL DOOUTY( 558 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10559 CALL DOINX( 559 )                                               
CY10559 CALL DOINY( 559 )                                               
      DO 40 K=LA1,M,LA                                                  
      KB=K+K-2                                                          
      C1=TRIGS(KB+1)                                                    
      S1=TRIGS(KB+2)                                                    
      DO 30 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 25 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+A(IB+I)                                           
      D(JA+J)=B(IA+I)+B(IB+I)                                           
      C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))                 
      D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))                 
      I=I+INC3                                                          
      J=J+INC4                                                          
   25 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   30 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
   40 CONTINUE                                                          
CX20560 CALL DOOUTX( 560 )                                              
CY20560 CALL DOOUTY( 560 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
C                                                                       
C     CODING FOR FACTOR 3                                               
C                                                                       
   50 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
      IC=IB+IINK                                                        
      JC=JB+JINK                                                        
CX10561 CALL DOINX( 561 )                                               
CY10561 CALL DOINY( 561 )                                               
      DO 60 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 55 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))                                 
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))                                 
      C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))) 
      C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))) 
      D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))) 
      D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))) 
      I=I+INC3                                                          
      J=J+INC4                                                          
   55 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   60 CONTINUE                                                          
CX20562 CALL DOOUTX( 562 )                                              
CY20562 CALL DOOUTY( 562 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10563 CALL DOINX( 563 )                                               
CY10563 CALL DOINY( 563 )                                               
      DO 80 K=LA1,M,LA                                                  
      KB=K+K-2                                                          
      KC=KB+KB                                                          
      C1=TRIGS(KB+1)                                                    
      S1=TRIGS(KB+2)                                                    
      C2=TRIGS(KC+1)                                                    
      S2=TRIGS(KC+2)                                                    
      DO 70 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 65 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))                                 
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))                                 
      C(JB+J)=                                                          
     *    C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      D(JB+J)=                                                          
     *    S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      C(JC+J)=                                                          
     *    C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      D(JC+J)=                                                          
     *    S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      I=I+INC3                                                          
      J=J+INC4                                                          
   65 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   70 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
   80 CONTINUE                                                          
CX20564 CALL DOOUTX( 564 )                                              
CY20564 CALL DOOUTY( 564 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
C                                                                       
C     CODING FOR FACTOR 4                                               
C                                                                       
   90 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
      IC=IB+IINK                                                        
      JC=JB+JINK                                                        
      ID=IC+IINK                                                        
      JD=JC+JINK                                                        
CX10565 CALL DOINX( 565 )                                               
CY10565 CALL DOINY( 565 )                                               
      DO 100 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 95 IJK=1,LOT                                                   
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))                       
      C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))                       
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))                       
      D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))                       
      C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))                       
      C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))                       
      D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))                       
      D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))                       
      I=I+INC3                                                          
      J=J+INC4                                                          
   95 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  100 CONTINUE                                                          
CX20566 CALL DOOUTX( 566 )                                              
CY20566 CALL DOOUTY( 566 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10567 CALL DOINX( 567 )                                               
CY10567 CALL DOINY( 567 )                                               
      DO 120 K=LA1,M,LA                                                 
      KB=K+K-2                                                          
      KC=KB+KB                                                          
      KD=KC+KB                                                          
      C1=TRIGS(KB+1)                                                    
      S1=TRIGS(KB+2)                                                    
      C2=TRIGS(KC+1)                                                    
      S2=TRIGS(KC+2)                                                    
      C3=TRIGS(KD+1)                                                    
      S3=TRIGS(KD+2)                                                    
      DO 110 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 105 IJK=1,LOT                                                  
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))                       
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))                       
      C(JC+J)=                                                          
     *    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))                      
     *   -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))                      
      D(JC+J)=                                                          
     *    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))                      
     *   +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))                      
      C(JB+J)=                                                          
     *    C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))                      
     *   -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))                      
      D(JB+J)=                                                          
     *    S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))                      
     *   +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))                      
      C(JD+J)=                                                          
     *    C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))                      
     *   -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))                      
      D(JD+J)=                                                          
     *    S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))                      
     *   +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))                      
      I=I+INC3                                                          
      J=J+INC4                                                          
  105 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  110 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
  120 CONTINUE                                                          
CX20568 CALL DOOUTX( 568 )                                              
CY20568 CALL DOOUTY( 568 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
C                                                                       
C     CODING FOR FACTOR 5                                               
C                                                                       
  130 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
      IC=IB+IINK                                                        
      JC=JB+JINK                                                        
      ID=IC+IINK                                                        
      JD=JC+JINK                                                        
      IE=ID+IINK                                                        
      JE=JD+JINK                                                        
CX10569 CALL DOINX( 569 )                                               
CY10569 CALL DOINY( 569 )                                               
      DO 140 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 135 IJK=1,LOT                                                  
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))               
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))               
      C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *  -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))              
      C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *  +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))              
      D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *  +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))              
      D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *  -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))              
      C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *  -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))              
      C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *  +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))              
      D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *  +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))              
      D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *  -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))              
      I=I+INC3                                                          
      J=J+INC4                                                          
  135 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  140 CONTINUE                                                          
CX20570 CALL DOOUTX( 570 )                                              
CY20570 CALL DOOUTY( 570 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10571 CALL DOINX( 571 )                                               
CY10571 CALL DOINY( 571 )                                               
      DO 160 K=LA1,M,LA                                                 
      KB=K+K-2                                                          
      KC=KB+KB                                                          
      KD=KC+KB                                                          
      KE=KD+KB                                                          
      C1=TRIGS(KB+1)                                                    
      S1=TRIGS(KB+2)                                                    
      C2=TRIGS(KC+1)                                                    
      S2=TRIGS(KC+2)                                                    
      C3=TRIGS(KD+1)                                                    
      S3=TRIGS(KD+2)                                                    
      C4=TRIGS(KE+1)                                                    
      S4=TRIGS(KE+2)                                                    
      DO 150 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 145 IJK=1,LOT                                                  
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))               
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))               
      C(JB+J)=                                                          
     *    C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      D(JB+J)=                                                          
     *    S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      C(JE+J)=                                                          
     *    C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      D(JE+J)=                                                          
     *    S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      C(JC+J)=                                                          
     *    C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      D(JC+J)=                                                          
     *    S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      C(JD+J)=                                                          
     *    C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      D(JD+J)=                                                          
     *    S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      I=I+INC3                                                          
      J=J+INC4                                                          
  145 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  150 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
  160 CONTINUE                                                          
CX20572 CALL DOOUTX( 572 )                                              
CY20572 CALL DOOUTY( 572 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE CPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA)
      DIMENSION A(N),B(N),C(N),D(N),TRIGS(N)                            
C                                                                       
C     SUBROUTINE "CPASSM" - MULTIPLE VERSION                            
C     PERFORMS ONE PASS THROUGH DATA                                    
C     AS PART OF MULTIPLE COMPLEX (INVERSE) FFT ROUTINE                 
C     A IS FIRST REAL INPUT VECTOR                                      
C     B IS FIRST IMAGINARY INPUT VECTOR                                 
C     C IS FIRST REAL OUTPUT VECTOR                                     
C     D IS FIRST IMAGINARY OUTPUT VECTOR                                
C     TRIGS IS PRECALCULATED TABLE OF SINES & COSINES                   
C     INC1 IS ADDRESSING INCREMENT FOR A AND B                          
C     INC2 IS ADDRESSING INCREMENT FOR C AND D                          
C     INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S                    
C     INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S                    
C     LOT IS THE NUMBER OF VECTORS                                      
C     N IS LENGTH OF VECTORS                                            
C     IFAC IS CURRENT FACTOR OF N                                       
C     LA IS PRODUCT OF PREVIOUS FACTORS                                 
C                                                                       
      DATA SIN36/-0.587785252292473/,COS36/0.809016994374947/,          
     *     SIN72/-0.951056516295154/,COS72/0.309016994374947/,          
     *     SIN60/-0.866025403784437/                                    
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
C                                                                       
      M=N/IFAC                                                          
      IINK=M*INC1                                                       
      JINK=LA*INC2                                                      
      JUMP=(IFAC-1)*JINK                                                
      IBASE=0                                                           
      JBASE=0                                                           
      IGO=IFAC-1                                                        
      IF (IGO.GT.4) THEN                                                
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      ENDIF                                                             
      GO TO (10,50,90,130),IGO                                          
C                                                                       
C     CODING FOR FACTOR 2                                               
C                                                                       
   10 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
CX10573 CALL DOINX( 573 )                                               
CY10573 CALL DOINY( 573 )                                               
      DO 20 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 15 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+A(IB+I)                                           
      D(JA+J)=B(IA+I)+B(IB+I)                                           
      C(JB+J)=A(IA+I)-A(IB+I)                                           
      D(JB+J)=B(IA+I)-B(IB+I)                                           
      I=I+INC3                                                          
      J=J+INC4                                                          
   15 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   20 CONTINUE                                                          
CX20574 CALL DOOUTX( 574 )                                              
CY20574 CALL DOOUTY( 574 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10575 CALL DOINX( 575 )                                               
CY10575 CALL DOINY( 575 )                                               
      DO 40 K=LA1,M,LA                                                  
      KB=K+K-2                                                          
      C1=TRIGS(KB+1)                                                    
      S1=-TRIGS(KB+2)                                                   
      DO 30 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 25 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+A(IB+I)                                           
      D(JA+J)=B(IA+I)+B(IB+I)                                           
      C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))                 
      D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))                 
      I=I+INC3                                                          
      J=J+INC4                                                          
   25 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   30 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
   40 CONTINUE                                                          
CX20576 CALL DOOUTX( 576 )                                              
CY20576 CALL DOOUTY( 576 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
C                                                                       
C     CODING FOR FACTOR 3                                               
C                                                                       
   50 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
      IC=IB+IINK                                                        
      JC=JB+JINK                                                        
CX10577 CALL DOINX( 577 )                                               
CY10577 CALL DOINY( 577 )                                               
      DO 60 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 55 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))                                 
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))                                 
      C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))) 
      C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))) 
      D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))) 
      D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))) 
      I=I+INC3                                                          
      J=J+INC4                                                          
   55 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   60 CONTINUE                                                          
CX20578 CALL DOOUTX( 578 )                                              
CY20578 CALL DOOUTY( 578 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10579 CALL DOINX( 579 )                                               
CY10579 CALL DOINY( 579 )                                               
      DO 80 K=LA1,M,LA                                                  
      KB=K+K-2                                                          
      KC=KB+KB                                                          
      C1=TRIGS(KB+1)                                                    
      S1=-TRIGS(KB+2)                                                   
      C2=TRIGS(KC+1)                                                    
      S2=-TRIGS(KC+2)                                                   
      DO 70 L=1,LA                                                      
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 65 IJK=1,LOT                                                   
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))                                 
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))                                 
      C(JB+J)=                                                          
     *    C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      D(JB+J)=                                                          
     *    S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      C(JC+J)=                                                          
     *    C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      D(JC+J)=                                                          
     *    S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      I=I+INC3                                                          
      J=J+INC4                                                          
   65 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
   70 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
   80 CONTINUE                                                          
CX20580 CALL DOOUTX( 580 )                                              
CY20580 CALL DOOUTY( 580 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
C                                                                       
C     CODING FOR FACTOR 4                                               
C                                                                       
   90 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
      IC=IB+IINK                                                        
      JC=JB+JINK                                                        
      ID=IC+IINK                                                        
      JD=JC+JINK                                                        
CX10581 CALL DOINX( 581 )                                               
CY10581 CALL DOINY( 581 )                                               
      DO 100 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 95 IJK=1,LOT                                                   
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))                       
      C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))                       
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))                       
      D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))                       
      C(JB+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))                       
      C(JD+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))                       
      D(JB+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))                       
      D(JD+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))                       
      I=I+INC3                                                          
      J=J+INC4                                                          
   95 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  100 CONTINUE                                                          
CX20582 CALL DOOUTX( 582 )                                              
CY20582 CALL DOOUTY( 582 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10583 CALL DOINX( 583 )                                               
CY10583 CALL DOINY( 583 )                                               
      DO 120 K=LA1,M,LA                                                 
      KB=K+K-2                                                          
      KC=KB+KB                                                          
      KD=KC+KB                                                          
      C1=TRIGS(KB+1)                                                    
      S1=-TRIGS(KB+2)                                                   
      C2=TRIGS(KC+1)                                                    
      S2=-TRIGS(KC+2)                                                   
      C3=TRIGS(KD+1)                                                    
      S3=-TRIGS(KD+2)                                                   
      DO 110 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 105 IJK=1,LOT                                                  
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))                       
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))                       
      C(JC+J)=                                                          
     *    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))                      
     *   -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))                      
      D(JC+J)=                                                          
     *    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))                      
     *   +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))                      
      C(JB+J)=                                                          
     *    C1*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))                      
     *   -S1*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))                      
      D(JB+J)=                                                          
     *    S1*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))                      
     *   +C1*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))                      
      C(JD+J)=                                                          
     *    C3*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))                      
     *   -S3*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))                      
      D(JD+J)=                                                          
     *    S3*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))                      
     *   +C3*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))                      
      I=I+INC3                                                          
      J=J+INC4                                                          
  105 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  110 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
  120 CONTINUE                                                          
CX20584 CALL DOOUTX( 584 )                                              
CY20584 CALL DOOUTY( 584 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
C                                                                       
C     CODING FOR FACTOR 5                                               
C                                                                       
  130 IA=1                                                              
      JA=1                                                              
      IB=IA+IINK                                                        
      JB=JA+JINK                                                        
      IC=IB+IINK                                                        
      JC=JB+JINK                                                        
      ID=IC+IINK                                                        
      JD=JC+JINK                                                        
      IE=ID+IINK                                                        
      JE=JD+JINK                                                        
CX10585 CALL DOINX( 585 )                                               
CY10585 CALL DOINY( 585 )                                               
      DO 140 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 135 IJK=1,LOT                                                  
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))               
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))               
      C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *  -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))              
      C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *  +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))              
      D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *  +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))              
      D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *  -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))              
      C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *  -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))              
      C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *  +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))              
      D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *  +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))              
      D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *  -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))              
      I=I+INC3                                                          
      J=J+INC4                                                          
  135 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  140 CONTINUE                                                          
CX20586 CALL DOOUTX( 586 )                                              
CY20586 CALL DOOUTY( 586 )                                              
      IF (LA.EQ.M) THEN                                                 
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      ENDIF                                                             
      LA1=LA+1                                                          
      JBASE=JBASE+JUMP                                                  
CX10587 CALL DOINX( 587 )                                               
CY10587 CALL DOINY( 587 )                                               
      DO 160 K=LA1,M,LA                                                 
      KB=K+K-2                                                          
      KC=KB+KB                                                          
      KD=KC+KB                                                          
      KE=KD+KB                                                          
      C1=TRIGS(KB+1)                                                    
      S1=-TRIGS(KB+2)                                                   
      C2=TRIGS(KC+1)                                                    
      S2=-TRIGS(KC+2)                                                   
      C3=TRIGS(KD+1)                                                    
      S3=-TRIGS(KD+2)                                                   
      C4=TRIGS(KE+1)                                                    
      S4=-TRIGS(KE+2)                                                   
      DO 150 L=1,LA                                                     
      I=IBASE                                                           
      J=JBASE                                                           
COMMENTED OUT DIR$ IVDEP                                                
      DO 145 IJK=1,LOT                                                  
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))               
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))               
      C(JB+J)=                                                          
     *    C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      D(JB+J)=                                                          
     *    S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      C(JE+J)=                                                          
     *    C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      D(JE+J)=                                                          
     *    S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) 
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))         
     *   +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) 
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))         
      C(JC+J)=                                                          
     *    C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      D(JC+J)=                                                          
     *    S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      C(JD+J)=                                                          
     *    C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      D(JD+J)=                                                          
     *    S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) 
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))         
     *   +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) 
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))         
      I=I+INC3                                                          
      J=J+INC4                                                          
  145 CONTINUE                                                          
      IBASE=IBASE+INC1                                                  
      JBASE=JBASE+INC2                                                  
  150 CONTINUE                                                          
      JBASE=JBASE+JUMP                                                  
  160 CONTINUE                                                          
CX20588 CALL DOOUTX( 588 )                                              
CY20588 CALL DOOUTY( 588 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE CFTFAX(N,IFAX,TRIGS)                                   
      DIMENSION IFAX(1),TRIGS(1)                                        
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
C                                                                       
C     THIS ROUTINE WAS MODIFIED FROM TEMPERTON"S ORIGINAL               
C     BY DAVE FULKER.  IT NO LONGER PRODUCES FACTORS IN ASCENDING       
C     ORDER, AND THERE ARE NONE OF THE ORIGINAL 'MODE' OPTIONS.         
C                                                                       
C ON INPUT     N                                                        
C               THE LENGTH OF EACH COMPLEX TRANSFORM TO BE PERFORMED    
C                                                                       
C               N MUST BE GREATER THAN 1 AND CONTAIN NO PRIME           
C               FACTORS GREATER THAN 5.                                 
C                                                                       
C ON OUTPUT    IFAX                                                     
C               IFAX(1)                                                 
C                 THE NUMBER OF FACTORS CHOSEN OR -99 IN CASE OF ERROR  
C               IFAX(2) THRU IFAX( IFAX(1)+1 )                          
C                 THE FACTORS OF N IN THE FOLLOWIN ORDER:  APPEARING    
C                 FIRST ARE AS MANY FACTORS OF 4 AS CAN BE OBTAINED.    
C                 SUBSEQUENT FACTORS ARE PRIMES, AND APPEAR IN          
C                 ASCENDING ORDER, EXCEPT FOR MULTIPLE FACTORS.         
C                                                                       
C              TRIGS                                                    
C               2N SIN AND COS VALUES FOR USE BY THE TRANSFORM ROUTINE  
C                                                                       
      CALL FACT(N,IFAX)                                                 
      K = IFAX(1)                                                       
      IF (K .LT. 1 .OR. IFAX(K+1) .GT. 5) IFAX(1) = -99                 
C     IF (IFAX(1) .LE. 0 )CALL ABORT  ***ABORT NOT ANSI                 
      IF (IFAX(1) .LE. 0) THEN                                          
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
      STOP 'ABORT'                                                      
      ENDIF                                                             
      CALL CFTRIG (N, TRIGS)                                            
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FACT(N,IFAX)                                           
C     FACTORIZATION ROUTINE THAT FIRST EXTRACTS ALL FACTORS OF 4        
      DIMENSION IFAX(1)                                                 
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
      IF (N.GT.1) GO TO 10                                              
      IFAX(1) = 0                                                       
      IF (N.LT.1) IFAX(1) = -99                                         
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
      RETURN                                                            
   10 NN=N                                                              
      K=1                                                               
C     TEST FOR FACTORS OF 4                                             
   20 IF (MOD(NN,4).NE.0) GO TO 30                                      
      K=K+1                                                             
      IFAX(K)=4                                                         
      NN=NN/4                                                           
      IF (NN.EQ.1) GO TO 80                                             
      GO TO 20                                                          
C     TEST FOR EXTRA FACTOR OF 2                                        
   30 IF (MOD(NN,2).NE.0) GO TO 40                                      
      K=K+1                                                             
      IFAX(K)=2                                                         
      NN=NN/2                                                           
      IF (NN.EQ.1) GO TO 80                                             
C     TEST FOR FACTORS OF 3                                             
   40 IF (MOD(NN,3).NE.0) GO TO 50                                      
      K=K+1                                                             
      IFAX(K)=3                                                         
      NN=NN/3                                                           
      IF (NN.EQ.1) GO TO 80                                             
      GO TO 40                                                          
C     NOW FIND REMAINING FACTORS                                        
   50 L=5                                                               
      MAX = SQRT(FLOAT(NN))                                             
      INC=2                                                             
C     INC ALTERNATELY TAKES ON VALUES 2 AND 4                           
   60 IF (MOD(NN,L).NE.0) GO TO 70                                      
      K=K+1                                                             
      IFAX(K)=L                                                         
      NN=NN/L                                                           
      IF (NN.EQ.1) GO TO 80                                             
      GO TO 60                                                          
   70 IF (L.GT.MAX) GO TO 75                                            
      L=L+INC                                                           
      INC=6-INC                                                         
      GO TO 60                                                          
   75 K = K+1                                                           
      IFAX(K) = NN                                                      
   80 IFAX(1)=K-1                                                       
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
C     IFAX(1) NOW CONTAINS NUMBER OF FACTORS                            
      RETURN                                                            
      END                                                               
      SUBROUTINE CFTRIG(N,TRIGS)                                        
      DIMENSION TRIGS(1)                                                
CX30033 CALL SBINX (33)                                                 
CY30033 CALL SBINY (33)                                                 
      PI=2.0*ASIN(1.0)                                                  
      DEL=(PI+PI)/FLOAT(N)                                              
      L=N+N                                                             
CX10589 CALL DOINX( 589 )                                               
CY10589 CALL DOINY( 589 )                                               
      DO 10 I=1,L,2                                                     
      ANGLE=0.5*FLOAT(I-1)*DEL                                          
      TRIGS(I)=COS(ANGLE)                                               
      TRIGS(I+1)=SIN(ANGLE)                                             
   10 CONTINUE                                                          
CX20590 CALL DOOUTX( 590 )                                              
CY20590 CALL DOOUTY( 590 )                                              
CX40034 CALL SBOUTX (34)                                                
CY40034 CALL SBOUTY (34)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GETWAV(IDEV,BUF,LSTART,LENGTH,ISTAT,NRECS,NLEFT)       
C     SUBROUTINE GETWAV(IDEV,BUF,LSTART,LENGTH,ISTAT)                   
      DIMENSION BUF(1)                                                  
      PARAMETER (MAXRL=962560)                                          
      INTEGER NWORD, I, J                                               
C CALCULATE NUMBER OF WORDS IN MAX RECORD LENGTH                        
CX30035 CALL SBINX (35)                                                 
CY30035 CALL SBINY (35)                                                 
      NWORD = MAXRL/8                                                   
C READ FULL RECORDS                                                     
CX10591 CALL DOINX( 591 )                                               
CY10591 CALL DOINY( 591 )                                               
      DO 100 I = 1,NRECS                                                
         READ(IDEV,REC=LSTART+I-1)(BUF((I-1)*NWORD+J),J=1,NWORD)        
C     READ(IDEV,REC=LSTART)(BUF(I),I=1,LENGTH)                          
 100  CONTINUE                                                          
C READ REMAINING WORDS                                                  
CX20592 CALL DOOUTX( 592 )                                              
CY20592 CALL DOOUTY( 592 )                                              
      IF (NLEFT.NE.0) THEN                                              
         READ(IDEV,REC=LSTART+NRECS) (BUF(NRECS*NWORD+J),J=1,NLEFT/8)   
      END IF                                                            
CX40036 CALL SBOUTX (36)                                                
CY40036 CALL SBOUTY (36)                                                
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PUTWAV(IDEV,BUF,LSTART,LENGTH,ISTAT,NRECS,NLEFT)       
C     SUBROUTINE PUTWAV(IDEV,BUF,LSTART,LENGTH,ISTAT)                   
      DIMENSION BUF(1)                                                  
      PARAMETER (MAXRL=962560)                                          
      INTEGER NWORD, I, J                                               
C CALCULATE NUMBER OF WORDS IN MAX RECORD LENGTH                        
CX30037 CALL SBINX (37)                                                 
CY30037 CALL SBINY (37)                                                 
      NWORD = MAXRL/8                                                   
C WRITE FULL RECORDS                                                    
CX10593 CALL DOINX( 593 )                                               
CY10593 CALL DOINY( 593 )                                               
      DO 100 I = 1,NRECS                                                
         WRITE(IDEV,REC=LSTART+I-1)(BUF((I-1)*NWORD+J),J=1,NWORD)       
C     WRITE(IDEV,REC=LSTART)(BUF(I),I=1,LENGTH)                         
 100  CONTINUE                                                          
CX20594 CALL DOOUTX( 594 )                                              
CY20594 CALL DOOUTY( 594 )                                              
C                                                                       
C WRITE REMAINING WORDS                                                 
      IF (NLEFT.NE.0) THEN                                              
         WRITE(IDEV,REC=LSTART+NRECS)(BUF(NRECS*NWORD+J),J=1,NLEFT/8)   
      END IF                                                            
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      RETURN                                                            
      END                                                               
C                                                                       
       SUBROUTINE FFT991(VIN,WORK,TR,IFF,INC1,INC2,N,NUM,ISIGN)         
       DIMENSION VIN(1),WORK(1),TR(1),IFF(1),IFAX(15)                   
CX30039 CALL SBINX (39)                                                 
CY30039 CALL SBINY (39)                                                 
CX10595 CALL DOINX( 595 )                                               
CY10595 CALL DOINY( 595 )                                               
C       NFAC=IFF(1)                                                     
C       DO 1 K=1,NFAC                                                   
C1      IFAX(K)=IFF(K+1)                                                
C---- PERFORM NUM FFTS                                                  
       DO 100 NN=1,NUM                                                  
       IP=(NN-1)*INC2+1                                                 
       IE2=N*INC1+IP                                                    
       IE1=IP+2*INC1                                                    
       IF (ISIGN.EQ.1) THEN                                             
       DO 220 I=IE1,IE2,INC1                                            
220     VIN(I-INC1)=VIN(I)                                              
       ENDIF                                                            
C      CALL RFFT(VIN(IP),WORK,N,TR,IFAX,NFAC,INC1,ISIGN)                
       CALL RFFT(VIN(IP),WORK,N,TR,IFF(2),IFF(1),INC1,ISIGN)            
       IF (ISIGN.EQ.-1) THEN                                            
       DO 110 I=IE2,IE1,-INC1                                           
C REARRANGE DATA TO FFT991 ORDER                                        
110    VIN(I)=VIN(I-INC1)                                               
C IMAGINARY NYQUIST                                                     
       VIN(IE2+INC1)=0.                                                 
C IMAGINARY DC                                                          
       VIN(IE1-INC1)=0.                                                 
       ENDIF                                                            
100    CONTINUE                                                         
CX20596 CALL DOOUTX( 596 )                                              
CY20596 CALL DOOUTY( 596 )                                              
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
       RETURN                                                           
       END                                                              
C---------------------------------------------------------------------  
        SUBROUTINE RFFT(A,C,N,TRIGS,IFAX,NFAC,ISKIP,ISIGN)              
        DIMENSION TRIGS(1)                                              
        INTEGER IFAX(1)                                                 
        DIMENSION A(1),C(1)                                             
CX30041 CALL SBINX (41)                                                 
CY30041 CALL SBINY (41)                                                 
        IF(ISIGN.GT.0) GO TO 100                                        
C        REAL TO COMPLEX                                                
        MA=N                                                            
        LA=1                                                            
         II=1                                                           
        NF=NFAC+1                                                       
CX10597 CALL DOINX( 597 )                                               
CY10597 CALL DOINY( 597 )                                               
        DO 102 I=1,NFAC                                                 
         IFAC=IFAX(NF-I)                                                
         MA=MA/IFAC                                                     
         MASKIP=MA*ISKIP                                                
         IF(II.GT.0) CALL RPASS(A,A(1+MASKIP*IFAC),C,C(1+MA),           
     .       TRIGS,TRIGS(2),IFAC,MA,LA,ISKIP,1,N)                       
         IF(II.LT.0) CALL RPASS(C,C(1+MA*IFAC),A,A(1+MASKIP),           
     .       TRIGS,TRIGS(2),IFAC,MA,LA,1,ISKIP,N)                       
         LA=LA*IFAC                                                     
         II=-II                                                         
102     CONTINUE                                                        
CX20598 CALL DOOUTX( 598 )                                              
CY20598 CALL DOOUTY( 598 )                                              
C       CALL REMARK(' AFTER RPASS')                                     
        IF(II.LT.0) THEN                                                
         JJ=0                                                           
CX10599 CALL DOINX( 599 )                                               
CY10599 CALL DOINY( 599 )                                               
         DO 101 J=1,N*ISKIP,ISKIP                                       
          JJ=JJ+1                                                       
          A(J)=C(JJ)                                                    
101      CONTINUE                                                       
CX20600 CALL DOOUTX( 600 )                                              
CY20600 CALL DOOUTY( 600 )                                              
        END IF                                                          
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
        RETURN                                                          
100        CONTINUE                                                     
C       ------ COMPLEX TO REAL ---------                                
        LA=1                                                            
         II=1                                                           
C        .. MULTIPLY DC (AND NYQUIST) BY 0.5 ACCORDING TO TEMPERTON. .. 
        A(1)=0.5*A(1)                                                   
        N1=(N-1)*ISKIP+1                                                
        IF(MOD(N,2).EQ.0)A(N1)=0.5*A(N1)                                
        CON=2./FLOAT(N)                                                 
        MA=N                                                            
CX10601 CALL DOINX( 601 )                                               
CY10601 CALL DOINY( 601 )                                               
        DO 103 I=1,NFAC                                                 
         IFAC=IFAX(I)                                                   
         MA=MA/IFAC                                                     
         LASKIP=LA*ISKIP                                                
         IF(II.GT.0) CALL RPASSM(A,A(1+LASKIP),C,C(1+LA*IFAC),          
     .        TRIGS,TRIGS(2),IFAC,LA,MA,ISKIP,1,N)                      
         IF(II.LT.0) CALL RPASSM(C,C(1+LA),A,A(1+LASKIP*IFAC),          
     .        TRIGS,TRIGS(2),IFAC,LA,MA,1,ISKIP,N)                      
         LA=LA*IFAC                                                     
         II=-II                                                         
103     CONTINUE                                                        
CX20602 CALL DOOUTX( 602 )                                              
CY20602 CALL DOOUTY( 602 )                                              
C       CALL REMARK(' AFTER RPASSM')                                    
        IF(II.LT.0) THEN                                                
         JJ=0                                                           
CX10603 CALL DOINX( 603 )                                               
CY10603 CALL DOINY( 603 )                                               
         DO 104 J=1,N*ISKIP,ISKIP                                       
          JJ=JJ+1                                                       
          A(J)=C(JJ)*CON                                                
104      CONTINUE                                                       
CX20604 CALL DOOUTX( 604 )                                              
CY20604 CALL DOOUTY( 604 )                                              
        ELSE                                                            
CX10605 CALL DOINX( 605 )                                               
CY10605 CALL DOINY( 605 )                                               
        DO 105 J=1,N*ISKIP,ISKIP                                        
         A(J)=A(J)*CON                                                  
105     CONTINUE                                                        
CX20606 CALL DOOUTX( 606 )                                              
CY20606 CALL DOOUTY( 606 )                                              
        END IF                                                          
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
        RETURN                                                          
        END                                                             
C-----------------------------------------------------------------      
C                                                                       
C       PERFORM COMPLEX TO REAL RFFT PASS BASED ON DECIMATION IN        
C        FREQUENCY                                                      
C                                                                       
C-----------------------------------------------------------------      
        SUBROUTINE RPASSM(A,B,C,D,TRIGS,TRIGSI,IFAC,LA,MA,              
     .      ISKIP,ISKIP1,N)                                             
        DIMENSION A(N),C(N),TRIGS(1),TRIGSI(1),B(1),D(1)                
        SAVE SIN45,SIN60,SIN72,SIN36,SQ54,SN3672,IFLAG                  
        DATA IFLAG /0/                                                  
CX30043 CALL SBINX (43)                                                 
CY30043 CALL SBINY (43)                                                 
        IF(IFLAG.EQ.0) THEN                                             
         IFLAG=1                                                        
         SIN45=0.5*SQRT(2.)                                             
         SIN60=SIN(3.141592653589793/3.)                                
         SIN72=SIN(3.141592653589793*72./180.)                          
         SIN36=SIN(3.141592653589793*36./180.)                          
         SQ54=0.25*SQRT(5.)                                             
         SN3672=SIN36/SIN72                                             
        END IF                                                          
        M=N/IFAC                                                        
        M2=M+M                                                          
        MSKIP=M*ISKIP                                                   
        M2SKIP=MSKIP+MSKIP                                              
        MHLF=M/2                                                        
        LASKIP=LA*ISKIP                                                 
        LASKP1=LA*ISKIP1                                                
        IB=LASKP1                                                       
        JUMP=(IFAC-1)*LASKP1                                            
        JUMP3=JUMP*3                                                    
        IJ=(N-LA)*ISKIP                                                 
        GO TO (100,200,300,400,500,600),IFAC                            
100        PRINT *,'**** RADIX ONE ENCOUNTERED ****'                    
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
        STOP 999                                                        
200        CONTINUE                                                     
C       ------- FACTOR TWO ------------                                 
C         CALCULATE DC                                                  
         L1=1                                                           
CX10607 CALL DOINX( 607 )                                               
CY10607 CALL DOINY( 607 )                                               
         DO 10 L=1,LASKIP,ISKIP                                         
          Y0=A(L)                                                       
          Y1=A(L+IJ)                                                    
          C(L1)=Y0+Y1                                                   
          C(L1+IB)=Y0-Y1                                                
          L1=L1+ISKIP1                                                  
10       CONTINUE                                                       
CX20608 CALL DOOUTX( 608 )                                              
CY20608 CALL DOOUTY( 608 )                                              
      IF(LA.EQ.M) THEN                                                  
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
         I=LASKIP+1                                                     
         J=LASKP1+1+JUMP                                                
         MM=MHLF                                                        
         IF(MOD(MA,2).EQ.0) MM=MM-LA                                    
CX10609 CALL DOINX( 609 )                                               
CY10609 CALL DOINY( 609 )                                               
        DO 11 K=LA,MM,LA                                                
         KK=K+K                                                         
         JN=IJ-KK*ISKIP                                                 
         KK=KK+1                                                        
         TREAL=TRIGS(KK)                                                
         TIMAG=TRIGSI(KK)                                               
         DO 12 L=1,LASKIP,ISKIP                                         
          A0REAL=A(I)                                                   
          LJN=L+JN                                                      
          ALJNRE=A(LJN)                                                 
          A0IMAG=B(I)                                                   
          ALJNIM=B(LJN)                                                 
          C(J)=A0REAL+ALJNRE                                            
          D(J)=A0IMAG-ALJNIM                                            
          JIB=J+IB                                                      
          C(JIB)=TREAL*(A0REAL-ALJNRE)-TIMAG*(A0IMAG+                   
     .        ALJNIM)                                                   
          D(JIB)=TIMAG*(A0REAL-ALJNRE)+                                 
     .       TREAL*(A0IMAG+ ALJNIM)                                     
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
12       CONTINUE                                                       
         J=J+JUMP3                                                      
         I=I+LASKIP                                                     
11      CONTINUE                                                        
CX20610 CALL DOOUTX( 610 )                                              
CY20610 CALL DOOUTY( 610 )                                              
CX10611 CALL DOINX( 611 )                                               
CY10611 CALL DOINY( 611 )                                               
        DO 13 L=1,LASKIP,ISKIP                                          
         C(J)=A(I)                                                      
         C(J+IB)=-B(I)                                                  
         I=I+ISKIP                                                      
         J=J+ISKIP1                                                     
13      CONTINUE                                                        
CX20612 CALL DOOUTX( 612 )                                              
CY20612 CALL DOOUTY( 612 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
        RETURN                                                          
C        --------- FACTOR THREE ------------                            
300        CONTINUE                                                     
        LA3SKP=LASKP1*3                                                 
        KEND=(N-LA)/2+LA                                                
        IF(MOD(MA,2).EQ.0) KEND=(N-2*LA)/2+LA+LA                        
        KDIF=KEND+KEND-M2-LA                                            
        IF(MOD(MA,2).EQ.0) KDIF=KDIF-LA                                 
        JMPLA3=JUMP+LA3SKP                                              
        IJ2=M2SKIP                                                      
        IJ1=IJ2-LASKIP                                                  
        IC=IB+LASKP1                                                    
C         CALCULATE DC                                                  
        L1=1                                                            
CX10613 CALL DOINX( 613 )                                               
CY10613 CALL DOINY( 613 )                                               
        DO 14 L=1,LASKIP,ISKIP                                          
         A0=A(L)                                                        
         A1=A(L+IJ1)                                                    
         B1=A(L+IJ2)                                                    
         T2=A0-0.5*A1                                                   
         T3=SIN60*B1                                                    
         C(L1)=A0+A1                                                    
         C(L1+IB)=T2-T3                                                 
         C(L1+IC)=T2+T3                                                 
         L1=L1+ISKIP1                                                   
14      CONTINUE                                                        
CX20614 CALL DOOUTX( 614 )                                              
CY20614 CALL DOOUTY( 614 )                                              
      IF(LA.EQ.M) THEN                                                  
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
        I=LASKIP+1                                                      
        J=LASKP1+1+JUMP                                                 
         MM=MHLF                                                        
         IF(MOD(MA,2).EQ.0) MM=MM-LA                                    
CX10615 CALL DOINX( 615 )                                               
CY10615 CALL DOINY( 615 )                                               
        DO 15 K=LA,MM,LA                                                
          K2=K+K                                                        
          KJ=KDIF-K2                                                    
          KJ=(KJ+KJ)*ISKIP                                              
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=TRIGSI(KK)                                             
          DO 16 L=1,LASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IM2SKP=I+M2SKIP                                              
           A1REAL=A(IM2SKP)                                             
           A1IMAG=B(IM2SKP)                                             
           JNL=KJ+I                                                     
           A2REAL=A(JNL)                                                
           A2IMAG=-B(JNL)                                               
           T1REAL=A1REAL+A2REAL                                         
           T1IMAG=A1IMAG+A2IMAG                                         
           T2REAL=A0REAL-0.5*T1REAL                                     
           T2IMAG=A0IMAG-0.5*T1IMAG                                     
           T3REAL=SIN60*(A1REAL-A2REAL)                                 
           T3IMAG=SIN60*(A1IMAG-A2IMAG)                                 
           C(J)=A0REAL+T1REAL                                           
           D(J)=A0IMAG+T1IMAG                                           
           X1REAL=T2REAL-T3IMAG                                         
           X1IMAG=T2IMAG+T3REAL                                         
           X2REAL=T2REAL+T3IMAG                                         
           X2IMAG=T2IMAG-T3REAL                                         
           JIB=J+IB                                                     
           C(JIB)=TR1REA*X1REAL-TR1IMA*X1IMAG                           
           D(JIB)=TR1IMA*X1REAL+TR1REA*X1IMAG                           
           JIC=J+IC                                                     
           C(JIC)=TR2REA*X2REAL-TR2IMA*X2IMAG                           
           D(JIC)=TR2IMA*X2REAL+TR2REA*X2IMAG                           
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
16        CONTINUE                                                      
          J=J+JMPLA3                                                    
          I=I+LASKIP                                                    
15       CONTINUE                                                       
CX20616 CALL DOOUTX( 616 )                                              
CY20616 CALL DOOUTY( 616 )                                              
      IF(MOD(MA,2).NE.0) THEN                                           
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
CX10617 CALL DOINX( 617 )                                               
CY10617 CALL DOINY( 617 )                                               
         DO 17 L=1,LASKIP,ISKIP                                         
          A0=A(I)                                                       
          B0=A(I+LASKIP)                                                
          IM2SKP=I+M2SKIP                                               
          A1=A(IM2SKP)                                                  
          T1=0.5*A0-A1                                                  
          C(J)=A0+A1                                                    
          C(J+IB)=T1-SIN60*B0                                           
          C(J+IC)=-T1-SIN60*B0                                          
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
17       CONTINUE                                                       
CX20618 CALL DOOUTX( 618 )                                              
CY20618 CALL DOOUTY( 618 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
        RETURN                                                          
C       ------ FACTOR 4 -----------                                     
400        CONTINUE                                                     
        LA4SKP=LASKP1*4                                                 
        JMPLA4=JUMP+LA4SKP                                              
        IJ2=M2SKIP                                                      
        IJ1=IJ2-LASKIP                                                  
        IC=IB+LASKP1                                                    
        ID=IC+LASKP1                                                    
        L1=1                                                            
CX10619 CALL DOINX( 619 )                                               
CY10619 CALL DOINY( 619 )                                               
        DO 20 L=1,LASKIP,ISKIP                                          
         Y0=A(L)                                                        
         Y1=A(L+IJ1)                                                    
         Y3=A(L+IJ2)                                                    
         Y2=A(L+IJ)                                                     
         T1=Y0+Y2                                                       
         T3=Y0-Y2                                                       
         C(L1)=T1+Y1                                                    
         C(L1+IB)=T3-Y3                                                 
         C(L1+IC)=T1-Y1                                                 
         C(L1+ID)=T3+Y3                                                 
         L1=L1+ISKIP1                                                   
20      CONTINUE                                                        
CX20620 CALL DOOUTX( 620 )                                              
CY20620 CALL DOOUTY( 620 )                                              
      IF(LA.EQ.M) THEN                                                  
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
        I=LASKIP+1                                                      
        J=LASKP1+1+JUMP                                                 
         MM=MHLF                                                        
         IF(MOD(MA,2).EQ.0) MM=MM-LA                                    
CX10621 CALL DOINX( 621 )                                               
CY10621 CALL DOINY( 621 )                                               
        DO 21 K=LA,MM,LA                                                
          K2=K+K                                                        
         JN=IJ-K2*ISKIP                                                 
         JM1=JN-M2SKIP                                                  
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR3REA=TRIGS(KK)                                              
          TR3IMA=TRIGSI(KK)                                             
          DO 22 L=1,LASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IM2SKP=I+M2SKIP                                              
           A1REAL=A(IM2SKP)                                             
           A1IMAG=B(IM2SKP)                                             
           JNL=JN+L                                                     
           A2REAL=A(JNL)                                                
           A2IMAG=-B(JNL)                                               
           JM1L=JM1+L                                                   
           A3REAL=A(JM1L)                                               
           A3IMAG=-B(JM1L)                                              
           T1REAL=A0REAL+A2REAL                                         
           T1IMAG=A0IMAG+A2IMAG                                         
           T2REAL=A1REAL+A3REAL                                         
           T2IMAG=A1IMAG+A3IMAG                                         
           T3REAL=A0REAL-A2REAL                                         
           T3IMAG=A0IMAG-A2IMAG                                         
           T4REAL=A1REAL-A3REAL                                         
           T4IMAG=A1IMAG-A3IMAG                                         
           C(J)=T1REAL+T2REAL                                           
           D(J  )=T1IMAG+T2IMAG                                         
           X1REAL=T3REAL-T4IMAG                                         
           X1IMAG=T3IMAG+T4REAL                                         
           X2REAL=T1REAL-T2REAL                                         
           X2IMAG=T1IMAG-T2IMAG                                         
           X3REAL=T3REAL+T4IMAG                                         
           X3IMAG=T3IMAG-T4REAL                                         
           JIB=J+IB                                                     
           C(JIB)=TR1REA*X1REAL-TR1IMA*X1IMAG                           
           D(JIB)=TR1IMA*X1REAL+TR1REA*X1IMAG                           
           JIC=J+IC                                                     
           C(JIC)=TR2REA*X2REAL-TR2IMA*X2IMAG                           
           D(JIC)=TR2IMA*X2REAL+TR2REA*X2IMAG                           
           JID=J+ID                                                     
           C(JID)=TR3REA*X3REAL-TR3IMA*X3IMAG                           
           D(JID)=TR3IMA*X3REAL+TR3REA*X3IMAG                           
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
22        CONTINUE                                                      
          J=J+JMPLA4                                                    
          I=I+LASKIP                                                    
21       CONTINUE                                                       
CX20622 CALL DOOUTX( 622 )                                              
CY20622 CALL DOOUTY( 622 )                                              
CX10623 CALL DOINX( 623 )                                               
CY10623 CALL DOINY( 623 )                                               
         DO 23 L=1,LASKIP,ISKIP                                         
          A0=A(I)                                                       
          B0=B(I)                                                       
          IM2SKP=I+M2SKIP                                               
          A1=A(IM2SKP)                                                  
          B1=B(IM2SKP)                                                  
          T1=SIN45*(B0+B1)                                              
          T2=SIN45*(A0-A1)                                              
          C(J)=A0+A1                                                    
          C(J+IB)=T2-T1                                                 
          C(J+IC)=B1-B0                                                 
          C(J+ID)=-(T2+T1)                                              
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
23       CONTINUE                                                       
CX20624 CALL DOOUTX( 624 )                                              
CY20624 CALL DOOUTY( 624 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
        RETURN                                                          
C        ------- FACTOR 5 ---------                                     
500        CONTINUE                                                     
        LA5SKP=LASKP1*5                                                 
        KEND=(N-LA)/2+LA                                                
        IF(MOD(MA,2).EQ.0) KEND=(N-2*LA)/2+LA+LA                        
        KDIF=KEND+KEND-M2-LA-M                                          
        IF(MOD(MA,2).EQ.0) KDIF=KDIF-LA                                 
        JMPLA5=JUMP+LA5SKP                                              
        IJ2=M2SKIP                                                      
        IJ1=IJ2-LASKIP                                                  
        IJ4=M2SKIP+M2SKIP                                               
        IJ3=IJ4-LASKIP                                                  
        IC=IB+LASKP1                                                    
        ID=IC+LASKP1                                                    
        IE=ID+LASKP1                                                    
C         CALCULATE DC                                                  
        L1=1                                                            
CX10625 CALL DOINX( 625 )                                               
CY10625 CALL DOINY( 625 )                                               
        DO 30 L=1,LASKIP,ISKIP                                          
         A0=A(L)                                                        
         A1=A(L+IJ1)                                                    
         B1=A(L+IJ2)                                                    
         A2=A(L+IJ3)                                                    
         B2=A(L+IJ4)                                                    
         T3=SIN72*B1                                                    
         T4=SIN72*B2                                                    
         T5=A1+A2                                                       
         T6=SQ54*(A1-A2)                                                
         T7=A0-0.25*T5                                                  
         T8=T7+T6                                                       
         T9=T7-T6                                                       
         T10=T3+SN3672*T4                                               
         T11=SN3672*T3-T4                                               
         C(L1)=A0+T5                                                    
         C(L1+IB)=T8-T10                                                
         C(L1+IC)=T9-T11                                                
         C(L1+ID)=T9+T11                                                
         C(L1+IE)=T8+T10                                                
         L1=L1+ISKIP1                                                   
30      CONTINUE                                                        
CX20626 CALL DOOUTX( 626 )                                              
CY20626 CALL DOOUTY( 626 )                                              
      IF(LA.EQ.M) THEN                                                  
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
        I=LASKIP+1                                                      
        J=LASKP1+1+JUMP                                                 
         MM=MHLF                                                        
         IF(MOD(MA,2).EQ.0) MM=MM-LA                                    
CX10627 CALL DOINX( 627 )                                               
CY10627 CALL DOINY( 627 )                                               
        DO 31 K=LA,MM,LA                                                
          K2=K+K                                                        
          KJ=KDIF-K2                                                    
          KJ=(KJ+KJ)*ISKIP                                              
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR3REA=TRIGS(KK)                                              
          TR3IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR4REA=TRIGS(KK)                                              
          TR4IMA=TRIGSI(KK)                                             
          DO 32 L=1,LASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IM2SKP=I+M2SKIP                                              
           A1REAL=A(IM2SKP)                                             
           A1IMAG=B(IM2SKP)                                             
           IM4=IM2SKP+M2SKIP                                            
           A2REAL=A(IM4)                                                
           A2IMAG=B(IM4)                                                
           JNL=KJ+I                                                     
           A3REAL=A(JNL)                                                
           A3IMAG=-B(JNL)                                               
           JNLM2S=JNL-M2SKIP                                            
           A4REAL=A(JNLM2S)                                             
           A4IMAG=-B(JNLM2S)                                            
           T1REAL=A1REAL+A4REAL                                         
           T1IMAG=A1IMAG+A4IMAG                                         
           T2REAL=A2REAL+A3REAL                                         
           T2IMAG=A2IMAG+A3IMAG                                         
           T3REAL=SIN72*(A1REAL-A4REAL)                                 
           T3IMAG=SIN72*(A1IMAG-A4IMAG)                                 
           T4REAL=SIN72*(A2REAL-A3REAL)                                 
           T4IMAG=SIN72*(A2IMAG-A3IMAG)                                 
           T5REAL=T1REAL+T2REAL                                         
           T5IMAG=T1IMAG+T2IMAG                                         
           T6REAL=SQ54*(T1REAL-T2REAL)                                  
           T6IMAG=SQ54*(T1IMAG-T2IMAG)                                  
           T7REAL=A0REAL-T5REAL*0.25                                    
           T7IMAG=A0IMAG-T5IMAG*0.25                                    
           T8REAL=T7REAL+T6REAL                                         
           T8IMAG=T7IMAG+T6IMAG                                         
           T9REAL=T7REAL-T6REAL                                         
           T9IMAG=T7IMAG-T6IMAG                                         
           T10REA=T3REAL+SN3672*T4REAL                                  
           T10IMA=T3IMAG+SN3672*T4IMAG                                  
           T11REA=SN3672*T3REAL-T4REAL                                  
           T11IMA=SN3672*T3IMAG-T4IMAG                                  
           C(J)=A0REAL+T5REAL                                           
           D(J)=A0IMAG+T5IMAG                                           
           X1REAL=T8REAL-T10IMA                                         
           X1IMAG=T8IMAG+T10REA                                         
           X2REAL=T9REAL-T11IMA                                         
           X2IMAG=T9IMAG+T11REA                                         
           X3REAL=T9REAL+T11IMA                                         
           X3IMAG=T9IMAG-T11REA                                         
           X4REAL=T8REAL+T10IMA                                         
           X4IMAG=T8IMAG-T10REA                                         
           JIB=J+IB                                                     
           C(JIB)=TR1REA*X1REAL-TR1IMA*X1IMAG                           
           D(JIB)=TR1IMA*X1REAL+TR1REA*X1IMAG                           
           JIC=J+IC                                                     
           C(JIC)=TR2REA*X2REAL-TR2IMA*X2IMAG                           
           D(JIC)=TR2IMA*X2REAL+TR2REA*X2IMAG                           
           JID=J+ID                                                     
           C(JID)=TR3REA*X3REAL-TR3IMA*X3IMAG                           
           D(JID)=TR3IMA*X3REAL+TR3REA*X3IMAG                           
           JIE=J+IE                                                     
           C(JIE)=TR4REA*X4REAL-TR4IMA*X4IMAG                           
           D(JIE)=TR4IMA*X4REAL+TR4REA*X4IMAG                           
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
32        CONTINUE                                                      
          J=J+JMPLA5                                                    
          I=I+LASKIP                                                    
31       CONTINUE                                                       
CX20628 CALL DOOUTX( 628 )                                              
CY20628 CALL DOOUTY( 628 )                                              
      IF(MOD(MA,2).NE.0) THEN                                           
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
CX10629 CALL DOINX( 629 )                                               
CY10629 CALL DOINY( 629 )                                               
         DO 33 L=1,LASKIP,ISKIP                                         
          A0=A(I)                                                       
          B0=B(I)                                                       
          IM2SKP=I+M2SKIP                                               
          A1=A(IM2SKP)                                                  
          B1=B(IM2SKP)                                                  
          IM4=IM2SKP+M2SKIP                                             
          A2=A(IM4)                                                     
          T1=A0+A1                                                      
          T2=0.25*T1-A2                                                 
          T3=SQ54*(A0-A1)                                               
          T4=SIN36*B0+SIN72*B1                                          
          T5=SIN72*B0-SIN36*B1                                          
          T6=T3+T2                                                      
          T7=T3-T2                                                      
          C(J)=T1+A2                                                    
          C(J+IB)=T6-T4                                                 
          C(J+IC)=T7-T5                                                 
          C(J+ID)=-T7-T5                                                
          C(J+IE)=-T6-T4                                                
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
33       CONTINUE                                                       
CX20630 CALL DOOUTX( 630 )                                              
CY20630 CALL DOOUTY( 630 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
        RETURN                                                          
C        ----- FACTOR 6 -----------------                               
600        CONTINUE                                                     
        LA6SKP=LASKP1*6                                                 
        JMPLA6=JUMP+LA6SKP                                              
        IJ2=M2SKIP                                                      
        IJ1=IJ2-LASKIP                                                  
        IC=IB+LASKP1                                                    
        ID=IC+LASKP1                                                    
        IE=ID+LASKP1                                                    
        IG=IE+LASKP1                                                    
        L1=1                                                            
CX10631 CALL DOINX( 631 )                                               
CY10631 CALL DOINY( 631 )                                               
         DO 41 L=1,LASKIP,ISKIP                                         
          A0=A(L)                                                       
          LIJ1=L+IJ1                                                    
          A1=A(LIJ1)                                                    
          B1=B(LIJ1)                                                    
          LIJ3=LIJ1+M2SKIP                                              
          A2=A(LIJ3)                                                    
          B2=B(LIJ3)                                                    
          A3=A(L+IJ)                                                    
           T2REAL=A0-A2*0.5                                             
           T3IMAG=SIN60*B2                                              
           Y0REAL=A0+A2                                                 
           Y4REAL=T2REAL-T3IMAG                                         
           Y2REAL=T2REAL+T3IMAG                                         
           T2REAL=A3-A1*0.5                                             
           T3IMAG=-SIN60*B1                                             
           Y3REAL=A3+A1                                                 
           Y1REAL=T2REAL-T3IMAG                                         
           Y5REAL=T2REAL+T3IMAG                                         
           X0REAL=Y0REAL+Y3REAL                                         
           X1REAL=Y4REAL-Y1REAL                                         
           X2REAL=Y2REAL+Y5REAL                                         
           X3REAL=Y0REAL-Y3REAL                                         
           X4REAL=Y4REAL+Y1REAL                                         
           X5REAL=Y2REAL-Y5REAL                                         
           C(L1   )=X0REAL                                              
           C(L1+IB)=X1REAL                                              
           C(L1+IC)=X2REAL                                              
           C(L1+ID)=X3REAL                                              
           C(L1+IE)=X4REAL                                              
           C(L1+IG)=X5REAL                                              
           L1=L1+ISKIP1                                                 
41       CONTINUE                                                       
CX20632 CALL DOOUTX( 632 )                                              
CY20632 CALL DOOUTY( 632 )                                              
      IF(LA.EQ.M) THEN                                                  
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
        I=LASKIP+1                                                      
        J=LASKP1+1+JUMP                                                 
         MM=MHLF                                                        
         IF(MOD(MA,2).EQ.0) MM=MM-LA                                    
CX10633 CALL DOINX( 633 )                                               
CY10633 CALL DOINY( 633 )                                               
        DO 42 K=LA,MM,LA                                                
          K2=K+K                                                        
         JN=IJ-K2*ISKIP                                                 
         JM1=JN-M2SKIP                                                  
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR3REA=TRIGS(KK)                                              
          TR3IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR4REA=TRIGS(KK)                                              
          TR4IMA=TRIGSI(KK)                                             
          KK=KK+K2                                                      
          TR5REA=TRIGS(KK)                                              
          TR5IMA=TRIGSI(KK)                                             
          DO 43 L=1,LASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IM2SKP=I+M2SKIP                                              
           A1REAL=A(IM2SKP)                                             
           A1IMAG=B(IM2SKP)                                             
           IM4=IM2SKP+M2SKIP                                            
           A2REAL=A(IM4)                                                
           A2IMAG=B(IM4)                                                
           JNL=JN+L                                                     
           A3REAL=A(JNL)                                                
           A3IMAG=-B(JNL)                                               
           JNLM2S=JNL-M2SKIP                                            
           A4REAL=A(JNLM2S)                                             
           A4IMAG=-B(JNLM2S)                                            
           JNLM4=JNLM2S-M2SKIP                                          
           A5REAL=A(JNLM4)                                              
           A5IMAG=-B(JNLM4)                                             
           T1REAL=A2REAL+A4REAL                                         
           T1IMAG=A2IMAG+A4IMAG                                         
           T2REAL=A0REAL-0.5*T1REAL                                     
           T2IMAG=A0IMAG-0.5*T1IMAG                                     
           T3REAL=SIN60*(A2REAL-A4REAL)                                 
           T3IMAG=SIN60*(A2IMAG-A4IMAG)                                 
           Y0REAL=A0REAL+T1REAL                                         
           Y0IMAG=A0IMAG+T1IMAG                                         
           Y4REAL=T2REAL-T3IMAG                                         
           Y4IMAG=T2IMAG+T3REAL                                         
           Y2REAL=T2REAL+T3IMAG                                         
           Y2IMAG=T2IMAG-T3REAL                                         
           T1REAL=A5REAL+A1REAL                                         
           T1IMAG=A5IMAG+A1IMAG                                         
           T2REAL=A3REAL-0.5*T1REAL                                     
           T2IMAG=A3IMAG-0.5*T1IMAG                                     
           T3REAL=SIN60*(A5REAL-A1REAL)                                 
           T3IMAG=SIN60*(A5IMAG-A1IMAG)                                 
           Y3REAL=A3REAL+T1REAL                                         
           Y3IMAG=A3IMAG+T1IMAG                                         
           Y1REAL=T2REAL-T3IMAG                                         
           Y1IMAG=T2IMAG+T3REAL                                         
           Y5REAL=T2REAL+T3IMAG                                         
           Y5IMAG=T2IMAG-T3REAL                                         
           X0REAL=Y0REAL+Y3REAL                                         
           X0IMAG=Y0IMAG+Y3IMAG                                         
           X4REAL=Y4REAL+Y1REAL                                         
           X4IMAG=Y4IMAG+Y1IMAG                                         
           X2REAL=Y2REAL+Y5REAL                                         
           X2IMAG=Y2IMAG+Y5IMAG                                         
           X3REAL=Y0REAL-Y3REAL                                         
           X3IMAG=Y0IMAG-Y3IMAG                                         
           X1REAL=Y4REAL-Y1REAL                                         
           X1IMAG=Y4IMAG-Y1IMAG                                         
           X5REAL=Y2REAL-Y5REAL                                         
           X5IMAG=Y2IMAG-Y5IMAG                                         
           C(J    )=X0REAL                                              
           D(J      )=X0IMAG                                            
           JIB=J+IB                                                     
           C(JIB)=TR1REA*X1REAL-TR1IMA*X1IMAG                           
           D(JIB)=TR1IMA*X1REAL+TR1REA*X1IMAG                           
           JIC=J+IC                                                     
           C(JIC)=TR2REA*X2REAL-TR2IMA*X2IMAG                           
           D(JIC)=TR2IMA*X2REAL+TR2REA*X2IMAG                           
           JID=J+ID                                                     
           C(JID)=TR3REA*X3REAL-TR3IMA*X3IMAG                           
           D(JID)=TR3IMA*X3REAL+TR3REA*X3IMAG                           
           JIE=J+IE                                                     
           C(JIE)=TR4REA*X4REAL-TR4IMA*X4IMAG                           
           D(JIE)=TR4IMA*X4REAL+TR4REA*X4IMAG                           
           JIG=J+IG                                                     
           C(JIG)=TR5REA*X5REAL-TR5IMA*X5IMAG                           
           D(JIG)=TR5IMA*X5REAL+TR5REA*X5IMAG                           
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
43        CONTINUE                                                      
          J=J+JMPLA6                                                    
          I=I+LASKIP                                                    
42       CONTINUE                                                       
CX20634 CALL DOOUTX( 634 )                                              
CY20634 CALL DOOUTY( 634 )                                              
      IF(MOD(MA,2).NE.0) THEN                                           
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
CX10635 CALL DOINX( 635 )                                               
CY10635 CALL DOINY( 635 )                                               
         DO 44 L=1,LASKIP,ISKIP                                         
          A0=A(I)                                                       
          B0=B(I)                                                       
          IM2SKP=I+M2SKIP                                               
          A1=A(IM2SKP)                                                  
          B1=B(IM2SKP)                                                  
          IM4=IM2SKP+M2SKIP                                             
          A2=A(IM4)                                                     
          B2=B(IM4)                                                     
          T1=A0+A2                                                      
          T2=B0+B2                                                      
          T3=SIN60*(A0-A2)                                              
          T4=SIN60*(B0-B2)                                              
          T5=0.5*T1-A1                                                  
          T6=0.5*T2+B1                                                  
          C(J)=A1+T1                                                    
          C(J+IB)=T3-T6                                                 
          C(J+IC)=T5-T4                                                 
          C(J+ID)=B1-T2                                                 
          C(J+IE)=-(T4+T5)                                              
          C(J+IG)=-(T3+T6)                                              
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
44       CONTINUE                                                       
CX20636 CALL DOOUTX( 636 )                                              
CY20636 CALL DOOUTY( 636 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
         RETURN                                                         
        END                                                             
C----------------------------------------------------------             
C                                                                       
C       PERFORM REAL TO COMPLEX RFFT PASS                               
C                                                                       
C----------------------------------------------------------             
        SUBROUTINE RPASS(A,B,C,D,TRIGS,TRIGSI,IFAC,MA,LA,               
     .    ISKIP,ISKIP1,N)                                               
        DIMENSION A(N),C(N),TRIGS(1),TRIGSI(1),B(1),D(1)                
        SAVE SIN45,SIN60,SIN72,SIN36,SQ54,SN3672,IFLAG                  
        DATA IFLAG /0/                                                  
CX30045 CALL SBINX (45)                                                 
CY30045 CALL SBINY (45)                                                 
        IF(IFLAG.EQ.0) THEN                                             
         IFLAG=1                                                        
         SIN45=0.5*SQRT(2.)                                             
         SIN60=SIN(3.141592653589793/3.)                                
         SIN72=SIN(3.141592653589793*72./180.)                          
         SIN36=SIN(3.141592653589793*36./180.)                          
         SQ54=0.25*SQRT(5.)                                             
         SN3672=SIN36/SIN72                                             
        END IF                                                          
        M=N/IFAC                                                        
        M2=M+M                                                          
        MSKIP=M*ISKIP                                                   
        M2SKIP=MSKIP+MSKIP                                              
        MSKIP1=M*ISKIP1                                                 
        M2SKP1=MSKIP1+MSKIP1                                            
        MHLF=M/2                                                        
        MASKIP=MA*ISKIP                                                 
        MASKP1=MA*ISKIP1                                                
        IB=MASKIP                                                       
        JUMP=(IFAC-1)*MASKIP                                            
        JUMP3=JUMP*3                                                    
        IJ=(N-MA)*ISKIP1                                                
        GO TO (100,200,300,400,500 ,600),IFAC                           
100        PRINT *,'**** RADIX ONE ENCOUNTERED ****'                    
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
        STOP 999                                                        
200        CONTINUE                                                     
C       -------- FACTOR TWO ----------                                  
         L1=1                                                           
CX10637 CALL DOINX( 637 )                                               
CY10637 CALL DOINY( 637 )                                               
         DO 10 L=1,MASKIP,ISKIP                                         
          A0REAL=A(L)                                                   
          A1REAL=A(L+IB)                                                
          C(L1)=A0REAL+A1REAL                                           
          C(L1+IJ)=A0REAL-A1REAL                                        
          L1=L1+ISKIP1                                                  
10       CONTINUE                                                       
CX20638 CALL DOOUTX( 638 )                                              
CY20638 CALL DOOUTY( 638 )                                              
         J=MASKP1+1                                                     
         I=MASKIP+1+JUMP                                                
      IF(MA.EQ.M) THEN                                                  
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
         MM=MHLF                                                        
         IF(MOD(LA,2).EQ.0) MM=MM-MA                                    
CX10639 CALL DOINX( 639 )                                               
CY10639 CALL DOINY( 639 )                                               
        DO 11 K=MA,MM,MA                                                
         KK=K+K                                                         
         JN=IJ-KK*ISKIP1                                                
         KK=KK+1                                                        
         TREAL=TRIGS(KK)                                                
         TIMAG=TRIGSI(KK)                                               
         L1=1                                                           
         DO 12 L=1,MASKIP,ISKIP                                         
          A0REAL=A(I)                                                   
          A0IMAG=B(I)                                                   
          IIB=I+IB                                                      
          A1REAL=A(IIB)                                                 
          A1IMAG=B(IIB)                                                 
          C(J)=A0REAL+A1REAL*TREAL+A1IMAG*TIMAG                         
          D(J)=A0IMAG-A1REAL*TIMAG+A1IMAG*TREAL                         
          JNL=JN+L1                                                     
          C(JNL)=A0REAL-A1REAL*TREAL-A1IMAG*TIMAG                       
          D(JNL)=-A0IMAG-A1REAL*TIMAG+                                  
     .      A1IMAG*TREAL                                                
          L1=L1+ISKIP1                                                  
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
12       CONTINUE                                                       
         I=I+JUMP3                                                      
         J=J+MASKP1                                                     
11      CONTINUE                                                        
CX20640 CALL DOOUTX( 640 )                                              
CY20640 CALL DOOUTY( 640 )                                              
      IF(MOD(LA,2).NE.0) THEN                                           
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
CX10641 CALL DOINX( 641 )                                               
CY10641 CALL DOINY( 641 )                                               
        DO 13 L=1,MASKIP,ISKIP                                          
         C(J)=A(I)                                                      
         D(J)=-A(I+IB)                                                  
         I=I+ISKIP                                                      
         J=J+ISKIP1                                                     
13      CONTINUE                                                        
CX20642 CALL DOOUTX( 642 )                                              
CY20642 CALL DOOUTY( 642 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
        RETURN                                                          
C       ------ FACTOR 3 -----------                                     
300        CONTINUE                                                     
        MA3SKP=MASKIP*3                                                 
        KEND=(N-MA)/2+MA                                                
        IF(MOD(LA,2).EQ.0) KEND=(N-2*MA)/2+MA*2                         
        KDIF=KEND+KEND-M2-MA                                            
        IF(MOD(LA,2).EQ.0) KDIF=KDIF-MA                                 
        JMPMA3=JUMP+MA3SKP                                              
        IJ2=M2SKP1                                                      
        IJ1=IJ2-MASKP1                                                  
        IC=IB+MASKIP                                                    
        L1=1                                                            
CX10643 CALL DOINX( 643 )                                               
CY10643 CALL DOINY( 643 )                                               
         DO 20 L=1,MASKIP,ISKIP                                         
          X0=A(L)                                                       
          X1=A(L+IB)                                                    
          X2=A(L+IC)                                                    
          T1=X1+X2                                                      
          T2=X0-0.5*T1                                                  
          T3=SIN60*(X1-X2)                                              
          C(L1)=X0+T1                                                   
          C(L1+IJ1)=T2                                                  
          C(L1+IJ2)=-T3                                                 
          L1=L1+ISKIP1                                                  
20       CONTINUE                                                       
CX20644 CALL DOOUTX( 644 )                                              
CY20644 CALL DOOUTY( 644 )                                              
      IF(MA.EQ.M) THEN                                                  
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
         J=MASKP1+1                                                     
         I=MASKIP+1+JUMP                                                
         MM=MHLF                                                        
         IF(MOD(LA,2).EQ.0) MM=MM-MA                                    
CX10645 CALL DOINX( 645 )                                               
CY10645 CALL DOINY( 645 )                                               
        DO 21 K=MA,MM,MA                                                
          K2=K+K                                                        
          KJ=KDIF-K2                                                    
          KJ=(KJ+KJ)*ISKIP1                                             
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=-TRIGSI(KK)                                            
          DO 22 L=1,MASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IIB=I+IB                                                     
           A1REAL=A(IIB)*TR1REA-B(IIB)*TR1IMA                           
           A1IMAG=B(IIB)*TR1REA+A(IIB)*TR1IMA                           
           IIC=I+IC                                                     
           A2REAL=A(IIC)*TR2REA-B(IIC)*TR2IMA                           
           A2IMAG=B(IIC)*TR2REA+A(IIC)*TR2IMA                           
           T1REAL=A1REAL+A2REAL                                         
           T1IMAG=A1IMAG+A2IMAG                                         
           T2REAL=A0REAL-T1REAL*0.5                                     
           T2IMAG=A0IMAG-T1IMAG*0.5                                     
           T3REAL=SIN60*(A1REAL-A2REAL)                                 
           T3IMAG=SIN60*(A1IMAG-A2IMAG)                                 
           C(J)=A0REAL+T1REAL                                           
           D(J  )=A0IMAG+T1IMAG                                         
           JM2SKP=J+M2SKP1                                              
           C(JM2SKP)=T2REAL+T3IMAG                                      
           D(JM2SKP)=T2IMAG-T3REAL                                      
           JNL=KJ+J                                                     
           C(JNL   )=T2REAL-T3IMAG                                      
           D(JNL)=-(T2IMAG+T3REAL)                                      
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
22        CONTINUE                                                      
          I=I+JMPMA3                                                    
          J=J+MASKP1                                                    
21       CONTINUE                                                       
CX20646 CALL DOOUTX( 646 )                                              
CY20646 CALL DOOUTY( 646 )                                              
      IF(MOD(LA,2).NE.0) THEN                                           
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
CX10647 CALL DOINX( 647 )                                               
CY10647 CALL DOINY( 647 )                                               
         DO 23 L=1,MASKIP,ISKIP                                         
          X0=A(I)                                                       
          X1=A(I+IB)                                                    
          X2=A(I+IC)                                                    
          T1=X1-X2                                                      
          C(J)=X0+T1*0.5                                                
          C(J+MASKP1)=-SIN60*(X1+X2)                                    
          JM2SKP=J+M2SKIP                                               
          C(JM2SKP)=X0-T1                                               
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
23       CONTINUE                                                       
CX20648 CALL DOOUTX( 648 )                                              
CY20648 CALL DOOUTY( 648 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
         RETURN                                                         
C       ------ FACTOR 4 -----------                                     
400        CONTINUE                                                     
        MA4SKP=MASKIP*4                                                 
        JMPMA4=JUMP+MA4SKP                                              
        IJ2=M2SKP1                                                      
        IJ1=IJ2-MASKP1                                                  
        IC=IB+MASKIP                                                    
        ID=IC+MASKIP                                                    
        L1=1                                                            
CX10649 CALL DOINX( 649 )                                               
CY10649 CALL DOINY( 649 )                                               
         DO 30 L=1,MASKIP,ISKIP                                         
          A0REAL=A(L)                                                   
          A1REAL=A(L+IB)                                                
          A2REAL=A(L+IC)                                                
          A3REAL=A(L+ID)                                                
          T1=A0REAL+A2REAL                                              
          T2=A1REAL+A3REAL                                              
          T3=A0REAL-A2REAL                                              
          T4=A3REAL-A1REAL                                              
          C(L1)=T1+T2                                                   
          C(L1+IJ1)=T3                                                  
          C(L1+IJ2)=T4                                                  
          C(L1+IJ)=T1-T2                                                
          L1=L1+ISKIP1                                                  
30       CONTINUE                                                       
CX20650 CALL DOOUTX( 650 )                                              
CY20650 CALL DOOUTY( 650 )                                              
      IF(MA.EQ.M) THEN                                                  
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
         J=MASKP1+1                                                     
         I=MASKIP+1+JUMP                                                
         MM=MHLF                                                        
         IF(MOD(LA,2).EQ.0) MM=MM-MA                                    
CX10651 CALL DOINX( 651 )                                               
CY10651 CALL DOINY( 651 )                                               
        DO 31 K=MA,MM,MA                                                
          K2=K+K                                                        
         JN=IJ-K2 *ISKIP1                                               
         JM1=JN-M2SKP1                                                  
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR3REA=TRIGS(KK)                                              
          TR3IMA=-TRIGSI(KK)                                            
          L1=1                                                          
          DO 32 L=1,MASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IIB=I+IB                                                     
           A1REAL=A(IIB)*TR1REA-B(IIB)*TR1IMA                           
           A1IMAG=B(IIB)*TR1REA+A(IIB)*TR1IMA                           
           IIC=I+IC                                                     
           A2REAL=A(IIC)*TR2REA-B(IIC)*TR2IMA                           
           A2IMAG=B(IIC)*TR2REA+A(IIC)*TR2IMA                           
           IID=I+ID                                                     
           A3REAL=A(IID)*TR3REA-B(IID)*TR3IMA                           
           A3IMAG=B(IID)*TR3REA+A(IID)*TR3IMA                           
           T1REAL=A0REAL+A2REAL                                         
           T1IMAG=A0IMAG+A2IMAG                                         
           T2REAL=A1REAL+A3REAL                                         
           T2IMAG=A1IMAG+A3IMAG                                         
           T3REAL=A0REAL-A2REAL                                         
           T3IMAG=A0IMAG-A2IMAG                                         
           T4REAL=A1REAL-A3REAL                                         
           T4IMAG=A1IMAG-A3IMAG                                         
           C(J)=T1REAL+T2REAL                                           
           D(J  )=T1IMAG+T2IMAG                                         
           JM2SKP=J+M2SKP1                                              
           C(JM2SKP)=T3REAL+T4IMAG                                      
           D(JM2SKP)=T3IMAG-T4REAL                                      
           JNL=JN+L1                                                    
           C(JNL   )=T1REAL-T2REAL                                      
           D(JNL)=T2IMAG-T1IMAG                                         
           JM1L=JM1+L1                                                  
           C(JM1L   )=T3REAL-T4IMAG                                     
           D(JM1L)=-T4REAL-T3IMAG                                       
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
           L1=L1+ISKIP1                                                 
32        CONTINUE                                                      
          I=I+JMPMA4                                                    
          J=J+MASKP1                                                    
31       CONTINUE                                                       
CX20652 CALL DOOUTX( 652 )                                              
CY20652 CALL DOOUTY( 652 )                                              
      IF(MOD(LA,2).NE.0) THEN                                           
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
CX10653 CALL DOINX( 653 )                                               
CY10653 CALL DOINY( 653 )                                               
         DO 33 L=1,MASKIP,ISKIP                                         
          X0=A(I)                                                       
          X1=A(I+IB)                                                    
          X2=A(I+IC)                                                    
          X3=A(I+ID)                                                    
          T1=SIN45*(X1-X3)                                              
          T2=SIN45*(X1+X3)                                              
          C(J)=X0+T1                                                    
          D(J)=-T2-X2                                                   
          JM2SKP=J+M2SKP1                                               
          C(JM2SKP)=X0-T1                                               
          D(JM2SKP)=X2-T2                                               
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
33       CONTINUE                                                       
CX20654 CALL DOOUTX( 654 )                                              
CY20654 CALL DOOUTY( 654 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
         RETURN                                                         
C        ------- FACTOR 5 ------------                                  
500        CONTINUE                                                     
        MA5SKP=MASKIP*5                                                 
        KEND=(N-MA)/2+MA                                                
        IF(MOD(LA,2).EQ.0) KEND=(N-2*MA)/2+MA*2                         
        KDIF=KEND+KEND-M2-MA-M                                          
        IF(MOD(LA,2).EQ.0) KDIF=KDIF-MA                                 
        JMPMA5=JUMP+MA5SKP                                              
        IJ2=M2SKP1                                                      
        IJ1=IJ2-MASKP1                                                  
        IJ4=IJ2+M2SKP1                                                  
        IJ3=IJ4-MASKP1                                                  
        IC=IB+MASKIP                                                    
        ID=IC+MASKIP                                                    
        IE=ID+MASKIP                                                    
        L1=1                                                            
CX10655 CALL DOINX( 655 )                                               
CY10655 CALL DOINY( 655 )                                               
         DO 40 L=1,MASKIP,ISKIP                                         
          X0=A(L)                                                       
          X1=A(L+IB)                                                    
          X2=A(L+IC)                                                    
          X3=A(L+ID)                                                    
          X4=A(L+IE)                                                    
          T1=X1+X4                                                      
          T2=X2+X3                                                      
          T3=SIN72*(X1-X4)                                              
          T4=SIN72*(X2-X3)                                              
          T5=T1+T2                                                      
          T6=SQ54*(T1-T2)                                               
          T7=X0-0.25*T5                                                 
          C(L1)=X0+T5                                                   
          C(L1+IJ1)=T7+T6                                               
          C(L1+IJ2)=-(T3+SN3672*T4)                                     
          C(L1+IJ3)=T7-T6                                               
          C(L1+IJ4)=-(SN3672*T3-T4)                                     
          L1=L1+ISKIP1                                                  
40       CONTINUE                                                       
CX20656 CALL DOOUTX( 656 )                                              
CY20656 CALL DOOUTY( 656 )                                              
      IF(MA.EQ.M) THEN                                                  
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
         J=MASKP1+1                                                     
         I=MASKIP+1+JUMP                                                
         MM=MHLF                                                        
         IF(MOD(LA,2).EQ.0) MM=MM-MA                                    
CX10657 CALL DOINX( 657 )                                               
CY10657 CALL DOINY( 657 )                                               
        DO 41 K=MA,MM,MA                                                
          K2=K+K                                                        
          KJ=KDIF-K2                                                    
          KJ=(KJ+KJ)*ISKIP1                                             
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR3REA=TRIGS(KK)                                              
          TR3IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR4REA=TRIGS(KK)                                              
          TR4IMA=-TRIGSI(KK)                                            
          DO 42 L=1,MASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IIB=I+IB                                                     
           A1REAL=A(IIB)*TR1REA-B(IIB)*TR1IMA                           
           A1IMAG=B(IIB)*TR1REA+A(IIB)*TR1IMA                           
           IIC=I+IC                                                     
           A2REAL=A(IIC)*TR2REA-B(IIC)*TR2IMA                           
           A2IMAG=B(IIC)*TR2REA+A(IIC)*TR2IMA                           
           IID=I+ID                                                     
           A3REAL=A(IID)*TR3REA-B(IID)*TR3IMA                           
           A3IMAG=B(IID)*TR3REA+A(IID)*TR3IMA                           
           IIE=I+IE                                                     
           A4REAL=A(IIE)*TR4REA-B(IIE)*TR4IMA                           
           A4IMAG=B(IIE)*TR4REA+A(IIE)*TR4IMA                           
           T1REAL=A1REAL+A4REAL                                         
           T1IMAG=A1IMAG+A4IMAG                                         
           T2REAL=A2REAL+A3REAL                                         
           T2IMAG=A2IMAG+A3IMAG                                         
           T3REAL=SIN72*(A1REAL-A4REAL)                                 
           T3IMAG=SIN72*(A1IMAG-A4IMAG)                                 
           T4REAL=SIN72*(A2REAL-A3REAL)                                 
           T4IMAG=SIN72*(A2IMAG-A3IMAG)                                 
           T5REAL=T1REAL+T2REAL                                         
           T5IMAG=T1IMAG+T2IMAG                                         
           T6REAL=SQ54*(T1REAL-T2REAL)                                  
           T6IMAG=SQ54*(T1IMAG-T2IMAG)                                  
           T7REAL=A0REAL-T5REAL*0.25                                    
           T7IMAG=A0IMAG-T5IMAG*0.25                                    
           T8REAL=T7REAL+T6REAL                                         
           T8IMAG=T7IMAG+T6IMAG                                         
           T9REAL=T7REAL-T6REAL                                         
           T9IMAG=T7IMAG-T6IMAG                                         
           T10REA=T3REAL+SN3672*T4REAL                                  
           T10IMA=T3IMAG+SN3672*T4IMAG                                  
           T11REA=SN3672*T3REAL-T4REAL                                  
           T11IMA=SN3672*T3IMAG-T4IMAG                                  
           C(J)=A0REAL+T5REAL                                           
           D(J)=A0IMAG+T5IMAG                                           
           JM2SKP=J+M2SKP1                                              
           C(JM2SKP)=T8REAL+T10IMA                                      
           D(JM2SKP)=T8IMAG-T10REA                                      
           JM4=JM2SKP+M2SKP1                                            
           C(JM4)=T9REAL+T11IMA                                         
           D(JM4)=T9IMAG-T11REA                                         
           JNL=KJ+J                                                     
           C(JNL   )=T9REAL-T11IMA                                      
           D(JNL)=-(T9IMAG+T11REA)                                      
           JNLM2S=JNL-M2SKP1                                            
           C(JNLM2S)=T8REAL-T10IMA                                      
           D(JNLM2S)=-(T8IMAG+T10REA)                                   
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
42        CONTINUE                                                      
          I=I+JMPMA5                                                    
          J=J+MASKP1                                                    
41       CONTINUE                                                       
CX20658 CALL DOOUTX( 658 )                                              
CY20658 CALL DOOUTY( 658 )                                              
      IF(MOD(LA,2).NE.0) THEN                                           
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
CX10659 CALL DOINX( 659 )                                               
CY10659 CALL DOINY( 659 )                                               
         DO 43 L=1,MASKIP,ISKIP                                         
          X0=A(I)                                                       
          X1=A(I+IB)                                                    
          X2=A(I+IC)                                                    
          X3=A(I+ID)                                                    
          X4=A(I+IE)                                                    
          T1=X1-X4                                                      
          T2=X1+X4                                                      
          T3=X2-X3                                                      
          T4=X2+X3                                                      
          T5=T1-T3                                                      
          T6=X0+0.25*T5                                                 
          T7=SQ54*(T1+T3)                                               
          C(J)=T6+T7                                                    
          D(J)=-SIN36*T2-SIN72*T4                                       
          JM2SKP=J+M2SKP1                                               
          C(JM2SKP)=T6-T7                                               
          D(JM2SKP)=-SIN72*T2+SIN36*T4                                  
          C(JM2SKP+M2SKP1)=X0-T5                                        
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
43       CONTINUE                                                       
CX20660 CALL DOOUTX( 660 )                                              
CY20660 CALL DOOUTY( 660 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
         RETURN                                                         
C        ---- FACTOR 6 -----------                                      
600        CONTINUE                                                     
        MA6SKP=MASKIP*6                                                 
        JMPMA6=JUMP+MA6SKP                                              
        IJ2=M2SKP1                                                      
        IJ1=IJ2-MASKP1                                                  
        IC=IB+MASKIP                                                    
        ID=IC+MASKIP                                                    
        IE=ID+MASKIP                                                    
        IG=IE+MASKIP                                                    
        L1=1                                                            
CX10661 CALL DOINX( 661 )                                               
CY10661 CALL DOINY( 661 )                                               
         DO 50 L=1,MASKIP,ISKIP                                         
          X0=A(L)                                                       
          X1=A(L+IB)                                                    
          X2=A(L+IC)                                                    
          X3=A(L+ID)                                                    
          X4=A(L+IE)                                                    
          X5=A(L+IG)                                                    
           T1REAL=X2+X4                                                 
           T2REAL=X0-0.5*T1REAL                                         
           T3REAL=SIN60*(X2-X4)                                         
           Y0REAL=X0+T1REAL                                             
           Y4REAL=T2REAL                                                
           Y4IMAG=-T3REAL                                               
           Y2REAL=T2REAL                                                
           Y2IMAG=T3REAL                                                
           T1REAL=X5+X1                                                 
           T2REAL=X3-0.5*T1REAL                                         
           T3REAL=SIN60*(X5-X1)                                         
           Y3REAL=X3+T1REAL                                             
           Y1REAL=T2REAL                                                
           Y1IMAG=-T3REAL                                               
           Y5REAL=T2REAL                                                
           Y5IMAG=T3REAL                                                
           C(L1    )=Y0REAL+Y3REAL                                      
           LIJ1=L1+IJ1                                                  
           C(LIJ1)=Y4REAL-Y1REAL                                        
           D(LIJ1)=Y4IMAG-Y1IMAG                                        
           LIJ3=LIJ1+M2SKP1                                             
           C(LIJ3)=Y2REAL+Y5REAL                                        
           D(LIJ3)=Y2IMAG+Y5IMAG                                        
           C(L1+IJ)=Y0REAL-Y3REAL                                       
           L1=L1+ISKIP1                                                 
50       CONTINUE                                                       
CX20662 CALL DOOUTX( 662 )                                              
CY20662 CALL DOOUTY( 662 )                                              
      IF(MA.EQ.M) THEN                                                  
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
         J=MASKP1+1                                                     
         I=MASKIP+1+JUMP                                                
         MM=MHLF                                                        
         IF(MOD(LA,2).EQ.0) MM=MM-MA                                    
CX10663 CALL DOINX( 663 )                                               
CY10663 CALL DOINY( 663 )                                               
        DO 51 K=MA,MM,MA                                                
          K2=K+K                                                        
         JN=IJ-K2*ISKIP1                                                
         JM1=JN-M2SKP1                                                  
          KK=K2+1                                                       
          TR1REA=TRIGS(KK)                                              
          TR1IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR2REA=TRIGS(KK)                                              
          TR2IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR3REA=TRIGS(KK)                                              
          TR3IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR4REA=TRIGS(KK)                                              
          TR4IMA=-TRIGSI(KK)                                            
          KK=KK+K2                                                      
          TR5REA=TRIGS(KK)                                              
          TR5IMA=-TRIGSI(KK)                                            
          DO 52 L=1,MASKIP,ISKIP                                        
           A0REAL=A(I)                                                  
           A0IMAG=B(I)                                                  
           IIB=I+IB                                                     
           A1REAL=A(IIB)*TR1REA-B(IIB)*TR1IMA                           
           A1IMAG=B(IIB)*TR1REA+A(IIB)*TR1IMA                           
           IIC=I+IC                                                     
           A2REAL=A(IIC)*TR2REA-B(IIC)*TR2IMA                           
           A2IMAG=B(IIC)*TR2REA+A(IIC)*TR2IMA                           
           IID=I+ID                                                     
           A3REAL=A(IID)*TR3REA-B(IID)*TR3IMA                           
           A3IMAG=B(IID)*TR3REA+A(IID)*TR3IMA                           
           IIE=I+IE                                                     
           A4REAL=A(IIE)*TR4REA-B(IIE)*TR4IMA                           
           A4IMAG=B(IIE)*TR4REA+A(IIE)*TR4IMA                           
           IIG=I+IG                                                     
           A5REAL=A(IIG)*TR5REA-B(IIG)*TR5IMA                           
           A5IMAG=B(IIG)*TR5REA+A(IIG)*TR5IMA                           
           T1REAL=A2REAL+A4REAL                                         
           T1IMAG=A2IMAG+A4IMAG                                         
           T2REAL=A0REAL-0.5*T1REAL                                     
           T2IMAG=A0IMAG-0.5*T1IMAG                                     
           T3REAL=SIN60*(A2REAL-A4REAL)                                 
           T3IMAG=SIN60*(A2IMAG-A4IMAG)                                 
           Y0REAL=A0REAL+T1REAL                                         
           Y0IMAG=A0IMAG+T1IMAG                                         
           Y4REAL=T2REAL+T3IMAG                                         
           Y4IMAG=T2IMAG-T3REAL                                         
           Y2REAL=T2REAL-T3IMAG                                         
           Y2IMAG=T2IMAG+T3REAL                                         
           T1REAL=A5REAL+A1REAL                                         
           T1IMAG=A5IMAG+A1IMAG                                         
           T2REAL=A3REAL-0.5*T1REAL                                     
           T2IMAG=A3IMAG-0.5*T1IMAG                                     
           T3REAL=SIN60*(A5REAL-A1REAL)                                 
           T3IMAG=SIN60*(A5IMAG-A1IMAG)                                 
           Y3REAL=A3REAL+T1REAL                                         
           Y3IMAG=A3IMAG+T1IMAG                                         
           Y1REAL=T2REAL+T3IMAG                                         
           Y1IMAG=T2IMAG-T3REAL                                         
           Y5REAL=T2REAL-T3IMAG                                         
           Y5IMAG=T2IMAG+T3REAL                                         
           C(J    )=Y0REAL+Y3REAL                                       
           D(J    )=Y0IMAG+Y3IMAG                                       
           JM2SKP=J+M2SKP1                                              
           C(JM2SKP    )=Y4REAL-Y1REAL                                  
           D(JM2SKP    )=Y4IMAG-Y1IMAG                                  
           JM4=JM2SKP+M2SKP1                                            
           C(JM4    )=Y2REAL+Y5REAL                                     
           D(JM4    )=Y2IMAG+Y5IMAG                                     
           JNL=JN+L                                                     
           C(JNL    )=Y0REAL-Y3REAL                                     
           D(JNL    )=-(Y0IMAG-Y3IMAG)                                  
           JNLM2S=JNL-M2SKP1                                            
           C(JNLM2S)=Y4REAL+Y1REAL                                      
           D(JNLM2S)=-(Y4IMAG+Y1IMAG)                                   
           JNLM4=JNLM2S-M2SKP1                                          
           C(JNLM4)=Y2REAL-Y5REAL                                       
           D(JNLM4)=-(Y2IMAG-Y5IMAG)                                    
           I=I+ISKIP                                                    
           J=J+ISKIP1                                                   
52        CONTINUE                                                      
          I=I+JMPMA6                                                    
          J=J+MASKP1                                                    
51       CONTINUE                                                       
CX20664 CALL DOOUTX( 664 )                                              
CY20664 CALL DOOUTY( 664 )                                              
      IF(MOD(LA,2).NE.0) THEN                                           
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      ENDIF                                                             
CX10665 CALL DOINX( 665 )                                               
CY10665 CALL DOINY( 665 )                                               
         DO 53 L=1,MASKIP,ISKIP                                         
          X0=A(I)                                                       
          X1=A(I+IB)                                                    
          X2=A(I+IC)                                                    
          X3=A(I+ID)                                                    
          X4=A(I+IE)                                                    
          X5=A(I+IG)                                                    
          T1=SIN60*(X5-X1)                                              
          T2=SIN60*(X2+X4)                                              
          T3=X2-X4                                                      
          T4=X1+X5                                                      
          T5=X0+0.5*T3                                                  
          T6=-X3-0.5*T4                                                 
          C(J)=T5-T1                                                    
          D(J)=T6-T2                                                    
          JM2SKP=J+M2SKP1                                               
          C(JM2SKP)=X0-T3                                               
          D(JM2SKP)=X3-T4                                               
          JM4=JM2SKP+M2SKP1                                             
          C(JM4)=T5+T1                                                  
          D(JM4)=T6+T2                                                  
          I=I+ISKIP                                                     
          J=J+ISKIP1                                                    
53       CONTINUE                                                       
CX20666 CALL DOOUTX( 666 )                                              
CY20666 CALL DOOUTY( 666 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
         RETURN                                                         
         END                                                            
      SUBROUTINE FFTFAX(N,IFAX,TRIGS)                                   
      DIMENSION IFAX(13),TRIGS(1)                                       
C                                                                       
      DATA MODE /3/                                                     
CX30047 CALL SBINX (47)                                                 
CY30047 CALL SBINY (47)                                                 
      CALL FAX (IFAX, N, MODE)                                          
      I = IFAX(1)                                                       
      IF (IFAX(I+1) .GT. 5) IFAX(1) = -99                               
C     IF (IFAX(1) .LE. 0 ) CALL ABORT  *** ABORT NOT ANSI               
      IF (IFAX(1) .LE. 0 ) THEN                                         
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
      STOP 'ABORT'                                                      
      ENDIF                                                             
      CALL FFTRIG (TRIGS, N, MODE)                                      
      CALL FFTTRG(N,TRIGS)                                              
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FAX(IFAX,N,MODE)                                       
      DIMENSION IFAX(10)                                                
CX30049 CALL SBINX (49)                                                 
CY30049 CALL SBINY (49)                                                 
      NN=N                                                              
      IF (ABS(MODE).EQ.1) GO TO 10                                      
      IF (ABS(MODE).EQ.8) GO TO 10                                      
      NN1=N/2                                                           
      IF ((NN1+NN1).EQ.N) GO TO 10                                      
      IFAX(1)=-99                                                       
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
      RETURN                                                            
   10 K=1                                                               
C     TEST FOR FACTORS OF 4                                             
   20 IF (MOD(NN,4).NE.0) GO TO 30                                      
      K=K+1                                                             
      IFAX(K)=4                                                         
      NN=NN/4                                                           
      IF (NN.EQ.1) GO TO 80                                             
      GO TO 20                                                          
C     TEST FOR EXTRA FACTOR OF 2                                        
   30 IF (MOD(NN,2).NE.0) GO TO 40                                      
      K=K+1                                                             
      IFAX(K)=2                                                         
      NN=NN/2                                                           
      IF (NN.EQ.1) GO TO 80                                             
C     TEST FOR FACTORS OF 3                                             
   40 IF (MOD(NN,3).NE.0) GO TO 50                                      
      K=K+1                                                             
      IFAX(K)=3                                                         
      NN=NN/3                                                           
      IF (NN.EQ.1) GO TO 80                                             
      GO TO 40                                                          
C     NOW FIND REMAINING FACTORS                                        
   50 L=5                                                               
      INC=2                                                             
C     INC ALTERNATELY TAKES ON VALUES 2 AND 4                           
   60 IF (MOD(NN,L).NE.0) GO TO 70                                      
      K=K+1                                                             
      IFAX(K)=L                                                         
      NN=NN/L                                                           
      IF (NN.EQ.1) GO TO 80                                             
      GO TO 60                                                          
   70 L=L+INC                                                           
      INC=6-INC                                                         
      GO TO 60                                                          
   80 IFAX(1)=K-1                                                       
C     IFAX(1) CONTAINS NUMBER OF FACTORS                                
      NFAX=IFAX(1)                                                      
C     SORT FACTORS INTO ASCENDING ORDER                                 
      IF (NFAX.EQ.1) GO TO 110                                          
CX10667 CALL DOINX( 667 )                                               
CY10667 CALL DOINY( 667 )                                               
      DO 100 II=2,NFAX                                                  
      ISTOP=NFAX+2-II                                                   
      DO 90 I=2,ISTOP                                                   
      IF (IFAX(I+1).GE.IFAX(I)) GO TO 90                                
      ITEM=IFAX(I)                                                      
      IFAX(I)=IFAX(I+1)                                                 
      IFAX(I+1)=ITEM                                                    
   90 CONTINUE                                                          
  100 CONTINUE                                                          
CX20668 CALL DOOUTX( 668 )                                              
CY20668 CALL DOOUTY( 668 )                                              
  110 CONTINUE                                                          
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FFTRIG(TRIGS,N,MODE)                                   
      DIMENSION TRIGS(1)                                                
CX30051 CALL SBINX (51)                                                 
CY30051 CALL SBINY (51)                                                 
      PI=2.0*ASIN(1.0)                                                  
      IMODE=ABS(MODE)                                                   
      NN=N                                                              
      IF (IMODE.GT.1.AND.IMODE.LT.6) NN=N/2                             
      DEL=(PI+PI)/FLOAT(NN)                                             
      L=NN+NN                                                           
CX10669 CALL DOINX( 669 )                                               
CY10669 CALL DOINY( 669 )                                               
      DO 10 I=1,L,2                                                     
      ANGLE=0.5*FLOAT(I-1)*DEL                                          
      TRIGS(I)=COS(ANGLE)                                               
      TRIGS(I+1)=SIN(ANGLE)                                             
   10 CONTINUE                                                          
CX20670 CALL DOOUTX( 670 )                                              
CY20670 CALL DOOUTY( 670 )                                              
      IF (IMODE.EQ.1) THEN                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      ENDIF                                                             
      IF (IMODE.EQ.8) THEN                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      ENDIF                                                             
      DEL=0.5*DEL                                                       
      NH=(NN+1)/2                                                       
      L=NH+NH                                                           
      LA=NN+NN                                                          
CX10671 CALL DOINX( 671 )                                               
CY10671 CALL DOINY( 671 )                                               
      DO 20 I=1,L,2                                                     
      ANGLE=0.5*FLOAT(I-1)*DEL                                          
      TRIGS(LA+I)=COS(ANGLE)                                            
      TRIGS(LA+I+1)=SIN(ANGLE)                                          
   20 CONTINUE                                                          
CX20672 CALL DOOUTX( 672 )                                              
CY20672 CALL DOOUTY( 672 )                                              
      IF (IMODE.LE.3) THEN                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      ENDIF                                                             
      DEL=0.5*DEL                                                       
      LA=LA+NN                                                          
      IF (MODE.EQ.5) GO TO 40                                           
CX10673 CALL DOINX( 673 )                                               
CY10673 CALL DOINY( 673 )                                               
      DO 30 I=2,NN                                                      
      ANGLE=FLOAT(I-1)*DEL                                              
      TRIGS(LA+I)=2.0*SIN(ANGLE)                                        
   30 CONTINUE                                                          
CX20674 CALL DOOUTX( 674 )                                              
CY20674 CALL DOOUTY( 674 )                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
   40 CONTINUE                                                          
      DEL=0.5*DEL                                                       
CX10675 CALL DOINX( 675 )                                               
CY10675 CALL DOINY( 675 )                                               
      DO 50 I=2,N                                                       
      ANGLE=FLOAT(I-1)*DEL                                              
      TRIGS(LA+I)=SIN(ANGLE)                                            
   50 CONTINUE                                                          
CX20676 CALL DOOUTX( 676 )                                              
CY20676 CALL DOOUTY( 676 )                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FFTTRG(N,TR)                                           
      DIMENSION TR(1)                                                   
      REAL CC                                                           
CX30053 CALL SBINX (53)                                                 
CY30053 CALL SBINY (53)                                                 
      DN = 3.141592653589793/N                                          
CX10677 CALL DOINX( 677 )                                               
CY10677 CALL DOINY( 677 )                                               
      DO 1 I=1,2*N,2                                                    
         CC = (I-1)*DN                                                  
         TR(I) = COS(CC)                                                
    1    TR(I+1) = SIN(CC)                                              
CX20678 CALL DOOUTX( 678 )                                              
CY20678 CALL DOOUTY( 678 )                                              
CX40054 CALL SBOUTX (54)                                                
CY40054 CALL SBOUTY (54)                                                
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE VERIFY(IDISK2,NX,DSECT)                                
      DIMENSION DSECT(NX,1)                                             
CX30055 CALL SBINX (55)                                                 
CY30055 CALL SBINY (55)                                                 
C                                                                       
      WRITE(10,1)                                                       
1     FORMAT(1X,
     * 'OUTPUT FOR PERFECT CLUB BENCHMARK: MG3D INSTRUMENTED VERSION'/ 
     *       1X,'$Revision: 1.A.1 $ $Author: psharma $'/)      
      WRITE(10,2)                                                       
2     FORMAT(1X,'VALIDATION PARAMETERS:'/)                              
C                                                                       
C     WRITE DEPTH LEVEL                                                 
C                                                                       
      WRITE(IDISK2,231)(DSECT(I, 61),I=1,NX)                            
231   FORMAT(1X,5E15.8)                                                 
      WRITE(10,231)(DSECT(I, 61),I=1,3)                                 
C                                                                       
      IVALID=0                                                          
C                                                                       
      IF(ABS((DSECT(1,61)-.70482932E+04)/.70482932E+04).GE.1E-2)        
     .   IVALID=1                                                       
      IF(ABS((DSECT(2,61)-.70231017E+04)/.70231017E+04).GE.1E-2)        
     .   IVALID=1                                                       
      IF(ABS((DSECT(3,61)-.69990724E+04)/.69990724E+04).GE.1E-2)        
     .   IVALID=1                                                       
C                                                                       
      IF(IVALID.EQ.0)THEN                                               
        WRITE(10,3)                                                     
3       FORMAT(//1X,'RESULTS FOR THIS RUN ARE:   VALID')                
      ELSE                                                              
        WRITE(10,4)                                                     
4       FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')                
      ENDIF                                                             
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
C                                                                       
      RETURN                                                            
      END                                                               
