      PROGRAM TSAFC                                                                                                                 
C     PROGRAM TSAFC (TSER,       FCOEFS,       INPUT,       OUTPUT,     )       H2                                                  
C    1         TAPE1=TSER, TAPE2=FCOEFS, TAPE5=INPUT, TAPE6=OUTPUT)                                                                 
C     -------------------------------------------------------------             H2                                                  
C                                                                               H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       H2                                                  
C     OCT 28/96 - F. MAJAESS (FORCE PROGRAM TO ABORT IF TIME SERIES LENGTH <4)                                                      
C     APR 26/94 - D. RAMSDEN  (INPUT TYPE TIME RATHER THAN GRID/ZONL/SUBA/COEF)                                                     
C     JUL 22/92 - E. CHAN  (REPLACE UNFORMATTED I/O WITH I/O ROUTINES)                                                              
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     APR 28/89 - B.DUGAS. (CORRECT THE OUTPUT NUMBER OF RECORDS PROCESSED)                                                         
C     AUG 15/85 - F. ZWIERS, B.DUGAS.                                                                                               
C                                                                               H2                                                  
CTSAFC   - FOURIER TRANSFORM A SET OF TIME SERIES RECORDS.              1  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHORS - F. ZWIERS AND B. DUGAS                                               H3                                                  
C                                                                               H3                                                  
CPURPOSE - COMPUTE FINITE FFT'S OF TIME SERIES.                                 H3                                                  
C          NOTE - THE LONGEST TIME SERIES THAT CAN BE PROCESSED HAS             H3                                                  
C                 10240 POINTS.                                                 H3                                                  
C                 PROGRAM WILL ABORT IF TIME SERIES LENGTH < 4.                 H3                                                  
C                                                                               H3                                                  
CINPUT FILE...                                                                  H3                                                  
C                                                                               H3                                                  
C      TSER   = FILE OF TIME SERIES WITH PROPER TIME SERIES LABELS.             H3                                                  
C                                                                               H3                                                  
COUTPUT FILE...                                                                 H3                                                  
C                                                                               H3                                                  
C      FCOEFS = CONTAINS THE TSER TIME SERIES FINITE FOURIER TRANSFORMS.        H3                                                  
C               THE RECORDS IN FCOEFS FILE HAVE LABEL:                          H3                                                  
C                                                                               H3                                                  
C                 IBUF(1) = TIME                                                H3                                                  
C                 IBUF(2) = LOCATION OF REAL POINT IN LATLON FORMAT IF          H3                                                  
C                           IKIND = 1, 2 OR 9 (GRID, ZONL OR SUBA),             H3                                                  
C                         = LOCATION OF COEFFICIENT IN THE REAL LA*2            H3                                                  
C                           LENGTH VECTOR IF IKIND = 3 (SPEC)                   H3                                                  
C                 IBUF(3) = NAME                                                H3                                                  
C                 IBUF(4) = LEVEL                                               H3                                                  
C                 IBUF(5) = LENGTH/2 + 1 WHERE LENGTH IS THE LENGTH OF          H3                                                  
C                           THE TRANSFORMED TIME SERIES .                       H3                                                  
C                 IBUF(6) = 1                                                   H3                                                  
C                 IBUF(7) = DIMENSION OF GRID, KHEM OR DIMENSION OF SPECTRAL    H3                                                  
C                           ARRAY AND IKIND (IN FORMAT CCCRRRKQ OR LRLMTQ)      H3                                                  
C                           WHERE Q IS IKIND (KIND OF DATA).                    H3                                                  
C                 IBUF(8) = PACKING DENSITY                                     H3                                                  
C                                                                                                                                   
CINPUT PARAMETER...                                                                                                                 
C                                                                               H5                                                  
C      LENGTH = THE ACTUAL LENGTH USED IN THE TRANSFORM. THIS NUMBER            H5                                                  
C               IS EVEN. THE REST (IF ANY) OF THE SERIES IS DISCARDED...        H5                                                  
C                                                                               H5                                                  
C               NOTE - IF THE INPUT PARAMETER LENGTH IS EXPRESSIBLE AS          H5                                                  
C                      A PRODUCT OF 2,3,4,5,6 AND/OR 8'S, THIS PROGRAM          H5                                                  
C                      USES  SUBROUTINE VFFT  TO COMPUTE THE  FFT'S OF          H5                                                  
C                      THE SERIES.                                              H5                                                  
C                    - IF ANY OTHER PRIME NUMBER COME INTO THE  DECOM-          H5                                                  
C                      POSITION OF LENGTH, ROUTINE FOURT IS USED TO DO          H5                                                  
C                      THE FFT'S.                                               H5                                                  
C                    - IF THE INPUT VALUE OF LENGTH  IS ZERO, THE PRO-          H5                                                  
C                      GRAM  WILL USE  THE ACTUAL  LENGTH OF THE FIRST          H5                                                  
C                      INPUT SERIES TO DETERMINE WHICH OF THE TWO PRE-          H5                                                  
C                      CEEDING CASES HOLD.                                      H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARD...                                                       H5                                                  
C                                                                               H5                                                  
C*TSAFC    9732                                                                 H5                                                  
C------------------------------------------------------------------------------                                                     
C                                                                                                                                   
                                                                                                                                    
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      DIMENSION X(111168),IX(512),WORK(111186),                                                                                     
     1          TRIGS(18528),IFAX(100),LFAX(7)                                                                                      
      INTEGER   ADD,BASE                                                                                                            
      LOGICAL   OK                                                                                                                  
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      DATA      MAXX/18528/,MAXL/18528/,MAXW/111168/,MAXS/64/,                                                                      
     1          LFAX / 6,8,5,4,3,2,1 /                                                                                              
C     DATA      NFOUR/4HFOUR/,NTIME/4HTIME/,NSPEC/4HSPEC/,                                                                          
                                                                                                                                    
C     * MAXL = LONGEST TIME SERIES THAT CAN BE PROCESSED                                                                            
C     * MAX  = IS LENGTH OF INPUT/OUTPUT BUFFER -- MUST BE AT LEAST                                                                 
C     *        MAXL+2                                                                                                               
C     * MAXW = IS LENGTH OF MEMORY AVAILABLE FOR STORING TIME SERIES                                                                
C     * MAXS = MAXIMUM NUMBER OF TIME SERIES WHICH CAN BE PROCESSED IN                                                              
C     *        ONE BLOCK -- NOTE THAT IX MUST BE DIMENSIONED MAXS*8                                                                 
C-----------------------------------------------------------------------------                                                      
      NFOUR=NC4TO8("FOUR")                                                                                                          
      NTIME=NC4TO8("TIME")                                                                                                          
      NSPEC=NC4TO8("SPEC")                                                                                                          
C                                                                                                                                   
      NF = 4                                                                                                                        
      CALL JCLPNT(NF,1,2,5,6)                                                                                                       
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
                                                                                                                                    
C     * READ UNIT 5 FOR LENGTH.                                                                                                     
                                                                                                                                    
      READ(5,5000,END=9000) LENGTH                                              H4                                                  
                                                                                                                                    
      IF (LENGTH.LT.4)                                         THEN                                                                 
          WRITE(6,6005) LENGTH                                                                                                      
          CALL                                     XIT('TSAFC',-1)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
C     * GET LABEL OF THE FIRST RECORD                                                                                               
                                                                                                                                    
      CALL FBUFFIN(1,IBUF,-8,K,LEN)                                                                                                 
      IF (K.GE.0) GOTO 9010                                                                                                         
      WRITE(6,6020) IBUF                                                                                                            
      REWIND 1                                                                                                                      
                                                                                                                                    
      IF ( IBUF(6).NE.1     .OR.                                                                                                    
     1     IBUF(1).NE.NTIME) THEN                                                                                                   
           WRITE(6,6030)                                                                                                            
           CALL                                    XIT('TSAFC',-2)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      NOUT = NFOUR                                                                                                                  
                                                                                                                                    
      LEN = IBUF(5)                                                                                                                 
                                                                                                                                    
C     * FORCE PROGRAM TO ABORT IF "LEN" < 4.                                                                                        
                                                                                                                                    
      IF (LEN.LT.4)                                            THEN                                                                 
          WRITE(6,6032) LENGTH                                                                                                      
          CALL                                     XIT('TSAFC',-3)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (LENGTH.EQ.0 .OR. LEN.LT.LENGTH)                      THEN                                                                 
          LENGTH = LEN                                                                                                              
          WRITE(6,6035) LENGTH                                                                                                      
      ENDIF                                                                                                                         
                                                                                                                                    
C     * TIME SERIES SHOULD BE OF EVEN LENGTH, SO                                                                                    
C     * DROP ONE OBSERVATION IF THE LENGTH IS ODD.                                                                                  
                                                                                                                                    
      IF (MOD(LENGTH,2).NE.0) LENGTH = LENGTH-1                                                                                     
                                                                                                                                    
      IF (LENGTH.GT.MAXL)                                      THEN                                                                 
          WRITE(6,6040)                                                                                                             
          CALL                                     XIT('TSAFC',-4)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      JUMP    = LENGTH+2                                                                                                            
      LENB2P1 = LENGTH/2 + 1                                                                                                        
      FLEN    = 1.0E0/FLOAT(LENGTH)                                                                                                 
      NSET    = -1                                                                                                                  
                                                                                                                                    
C     * TIME SERIES WILL BE PROCESSED IN GROUPS OF NBLOCK AT ONE                                                                    
C     * TIME TO TAKE ADVANTAGE OF THE VECTORIZATION OF THE FFT.                                                                     
                                                                                                                                    
      NBLOCK = MIN( (MAXW-LEN)/JUMP+1, MAXS)                                                                                        
                                                                                                                                    
C     * FIND THE FACTOR DECOMPOSITION OF LENGTH FOR VFFT.                                                                           
                                                                                                                                    
      NU     = LENGTH                                                                                                               
      IFAC   = 6                                                                                                                    
      K      = 0                                                                                                                    
      L      = 1                                                                                                                    
   20 IF (MOD(NU,IFAC).NE.0) GOTO 30                                                                                                
      NU     = NU/IFAC                                                                                                              
      IF (NU  .EQ.1) GOTO 40                                                                                                        
      IF (IFAC.NE.8) GOTO 20                                                                                                        
   30 L      = L+1                                                                                                                  
      IFAC   = LFAX(L)                                                                                                              
      IF (IFAC.NE.1) GOTO 20                                                                                                        
      ITYPE  = 1                                                                                                                    
      GOTO 1000                                                                                                                     
   40 ITYPE  = 0                                                                                                                    
                                                                                                                                    
      CALL FTSETUP(TRIGS,IFAX,LENGTH)                                                                                               
C-----------------------------------------------------------------------------                                                      
                                                                                                                                    
C     * MAIN LOOP                                                                                                                   
                                                                                                                                    
 1000 NSET = NSET+1                                                                                                                 
                                                                                                                                    
C        * LOAD A BLOCK OF TIME SERIES                                                                                              
                                                                                                                                    
         DO 1050 N=1,NBLOCK                                                                                                         
             ADD = (N-1)*JUMP + 1                                                                                                   
             CALL GETFLD2(1,X(ADD),-1,-1,-1,-1,IBUF,MAXX,OK)                                                                        
             IF (.NOT.OK) GOTO 1100                                                                                                 
                                                                                                                                    
             IF ( IBUF(6).NE.1     .OR.                                                                                             
     1            IBUF(1).NE.NTIME ) THEN                                                                                           
                  WRITE(6,6050) NSET*NBLOCK + N                                                                                     
                  WRITE(6,6020) IBUF                                                                                                
                  WRITE(6,6060)                                                                                                     
                  CALL                             XIT('TSAFC',-5)                                                                  
             ENDIF                                                                                                                  
                                                                                                                                    
             IF (IBUF(5).NE.LEN)                               THEN                                                                 
                 WRITE(6,6050) NSET*NBLOCK + N                                                                                      
                 WRITE(6,6020) IBUF                                                                                                 
                 WRITE(6,6070)                                                                                                      
                 CALL                              XIT('TSAFC',-6)                                                                  
             ENDIF                                                                                                                  
                                                                                                                                    
             BASE=(N-1)*8                                                                                                           
             DO 1050 J=1,8                                                                                                          
                 IX(BASE+J) = IBUF(J)                                                                                               
 1050    CONTINUE                                                                                                                   
                                                                                                                                    
 1100    NS = N-1                                                                                                                   
                                                                                                                                    
C        * PROCESS THIS BLOCK OF TIME SERIES -- THE                                                                                 
C        * CURRENT BLOCK CONTAINS NS TIME SERIES.                                                                                   
                                                                                                                                    
         IF (NS.EQ.0) GOTO 2000                                                                                                     
                                                                                                                                    
C        * PAD THE END OF EACH SERIES WITH TWO ZEROS BECAUSE                                                                        
C        * THE FFT IS RETURNED IN THE SAME SPACE AND IS TWO                                                                         
C        * WORDS LONGER THAN THE INPUTTED TIME SERIES.                                                                              
                                                                                                                                    
         DO 1150 I=LENGTH+1,NS*JUMP-1,JUMP                                                                                          
             X(I)   = 0.0E0                                                                                                         
             X(I+1) = 0.0E0                                                                                                         
 1150    CONTINUE                                                                                                                   
                                                                                                                                    
C        * FFT THE BLOCK                                                                                                            
                                                                                                                                    
         IF (ITYPE.EQ.0)                                       THEN                                                                 
             CALL VFFT(X,WORK,TRIGS,IFAX,1,JUMP,LENGTH,NS,-1)                                                                       
         ELSE                                                                                                                       
             DO 1170 N=1,NS                                                                                                         
                  ADD = (N-1)*JUMP+1                                                                                                
                  CALL FOURT(X(ADD),LENGTH,1,-1,0,WORK,LENGTH)                                                                      
                  DO 1170 I=ADD,ADD+JUMP-1                                                                                          
                      X(I) = X(I)*FLEN                                                                                              
 1170        CONTINUE                                                                                                               
         ENDIF                                                                                                                      
                                                                                                                                    
C        * OUTPUT THE BLOCK                                                                                                         
                                                                                                                                    
         IBUF(1) = NOUT                                                                                                             
         DO 1200 N=1,NS                                                                                                             
             BASE = (N-1)*8                                                                                                         
             DO 1190 J=2,8                                                                                                          
                 IBUF(J) = IX(BASE+J)                                                                                               
 1190        CONTINUE                                                                                                               
             IF (IBUF(1).EQ.NSPEC) IBUF(8) = 1                                                                                      
             IBUF(5) = LENB2P1                                                                                                      
             ADD     = (N-1)*JUMP+1                                                                                                 
             CALL PUTFLD2(2,X(ADD),IBUF,MAXX)                                                                                       
 1200    CONTINUE                                                                                                                   
         IF (NS.EQ.NBLOCK) GOTO 1000                                                                                                
C-----------------------------------------------------------------------------                                                      
                                                                                                                                    
 2000 CONTINUE                                                                                                                      
                                                                                                                                    
C     * ALL DONE.                                                                                                                   
                                                                                                                                    
                   NSO = NS                                                                                                         
      IF (NS.EQ.0) NSO = NBLOCK                                                                                                     
      BASE             = (NSO-1)*8                                                                                                  
      IX(BASE+1)       = NOUT                                                                                                       
      IX(BASE+5)       = LENB2P1                                                                                                    
      IF (NOUT.EQ.NSPEC)                                       THEN                                                                 
          IX(BASE+8)   = 1                                                                                                          
      ENDIF                                                                                                                         
      WRITE(6,6020) (IX(BASE+I),I=1,8)                                                                                              
      WRITE(6,6080) NSET*NBLOCK+NS                                                                                                  
      CALL                                         XIT('TSAFC',0)                                                                   
                                                                                                                                    
C     * PREMATURE EOF ON UNIT 1.                                                                                                    
                                                                                                                                    
 9000 WRITE(6,6000)                                                                                                                 
      CALL                                         XIT('TSAFC',-7)                                                                  
                                                                                                                                    
C     * PREMATURE EOF ON UNIT 1.                                                                                                    
                                                                                                                                    
 9010 WRITE(6,6010)                                                                                                                 
      CALL                                         XIT('TSAFC',-8)                                                                  
C-----------------------------------------------------------------------------                                                      
                                                                                                                                    
 5000 FORMAT(10X,I5)                                                            H4                                                  
 6000 FORMAT('0INPUT FILE IS EMPTY.')                                                                                               
 6005 FORMAT('0LENGTH VALUE IS ',I6)                                                                                                
 6010 FORMAT('0INPUT FILE OF TIME SERIES IS EMPTY.')                                                                                
 6020 FORMAT(1X,A4,I10,1X,A4,5I10)                                                                                                  
 6030 FORMAT('0FIRST RECORD IN INPUT FILE OF TIME SERIES',                                                                          
     1       ' DOES NOT HAVE A TIME SERIES LABEL.')                                                                                 
 6032 FORMAT('0SORRY, LENGTH OF TIME SERIES = ',I6,' < 4!')                                                                         
 6035 FORMAT('0LENGTH SET TO IBUF(5) =',I5)                                                                                         
 6040 FORMAT('0A TIME SERIES IN THE INPUT FILE IS TOO LONG.')                                                                       
 6050 FORMAT('0LABEL FOR RECORD ',I5,':')                                                                                           
 6060 FORMAT(11X,'IS NOT A TIME SERIES LABEL.')                                                                                     
 6070 FORMAT('0CHANGE OF LENGTH DETECTED IN THE INPUT TIME SERIES.')                                                                
 6080 FORMAT('0PROCESSED ',I10,' RECORDS.')                                                                                         
      END                                                                                                                           
