C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * *            
C     SUBROUTINE FOR DIFFERENTIATION BY FOURIER TRANSFORM               
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * *            
C                                                                       
C                                       CODED BY Y.OHSAKI               
C                                                                       
C     PURPOSE                                                           
C       TO DIFFERENTIATE AN EQUI-SPACED TIME HISTORY BY APPLICATION OF  
C       FAST FOURIER TRANSFORM                                          
C                                                                       
C     USAGE                                                             
C       CALL DIFR(N,X,ND,DT)                                            
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       N     - TOTAL NUMBER OF DATA  N.LE.8192                         
C       X(ND) - ORIGINAL/DIFFERENTIATED DATA AT CALL/RETURN             
C       ND    - DIMENSION OF X IN CALLING PROGRAM                       
C       DT    - TIME INCREMENT IN DATA                                  
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       FAST                                                            
C                                                                       
      SUBROUTINE DIFR(N,X,ND,DT)                                        
C                                                                       
      COMPLEX   C(8192)                                                 
      DIMENSION X(ND)                                                   
      PARAMETER (P2=6.283185)                                           
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      DO 110 M=1,N                                                      
      C(M)=CMPLX(X(M),0.)                                               
  110 CONTINUE                                                          
      NT=2                                                              
  120 IF(NT.GE.N) GO TO 130                                             
      NT=NT*2                                                           
      GO TO 120                                                         
  130 IF(NT.EQ.N) GO TO 150                                             
      DO 140 M=N+1,NT                                                   
      C(M)=(0.,0.)                                                      
  140 CONTINUE                                                          
  150 NFOLD=NT/2+1                                                      
C                                                                       
C     FOURIER TRANSFORM                                                 
C                                                                       
      CALL FAST(NT,C,8192,-1)                                           
C                                                                       
C     DIFFERENTIATION                                                   
C                                                                       
      C(1)=(0.,0.)                                                      
      DO 160 K=2,NFOLD-1                                                
      C(K)=CMPLX(0.,REAL(K-1))*C(K)                                     
      C(NT-K+2)=CONJG(C(K))                                             
  160 CONTINUE                                                          
      C(NFOLD)=(0.,0.)                                                  
C                                                                       
C     FOURIER INVERSE TRANSFORM                                         
C                                                                       
      CALL FAST(NT,C,8192,+1)                                           
      P=P2/REAL(NT)**2/DT                                               
      DO 170 M=1,N                                                      
      X(M)=REAL(C(M))*P                                                 
  170 CONTINUE                                                          
      RETURN                                                            
      END                                                               
                                                                        
