      PROGRAM FCATS                                                                                                                 
C     PROGRAM FCATS (FCOEFS,       TSER,       INPUT,       OUTPUT,     )       H2                                                  
C    1         TAPE1=FCOEFS, TAPE2=TSER, TAPE5=INPUT, TAPE6=OUTPUT)                                                                 
C     -------------------------------------------------------------             H2                                                  
C                                                                               H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       H2                                                  
C     APR 26/94 - D. RAMSDEN  (OUTPUT 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     AUG 26/85 - F. ZWIERS, B.DUGAS.                                                                                               
C                                                                               H2                                                  
CFCATS   - INVERSE FOURIER TRANSFORM TO RETREIVE TIME SERIES FILE.      1  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHORS - F. ZWIERS AND B. DUGAS                                               H3                                                  
C                                                                               H3                                                  
CPURPOSE -  COMPUTE ACTUAL TIME SERIES FROM TIME SERIES FOURIER                 H3                                                  
C           COEFFICIENTS.                                                       H3                                                  
C           NOTE - IF INPUT TYPE OF FCOEFS IS 4HSPEC THE DATA PACKING HAS       H3                                                  
C                  TO BE OFF OR THE FIRST LM IMAGINARY COMPONENTS WILL BE       H3                                                  
C                  LOST. IN THIS CASE THE PROGRAM WILL PRINT A WARNING...       H3                                                  
C                                                                               H3                                                  
CINPUT FILE...                                                                  H3                                                  
C                                                                               H3                                                  
C      FCOEFS = A FILE OF TIME SERIES FOURIER COEFFICIENTS (TYPE FOUR).         H3                                                  
C                                                                               H3                                                  
COUTPUT FILE...                                                                 H3                                                  
C                                                                               H3                                                  
C      TSER   = THE FINITE FOURIER TRANSFORMS OF FILE FCOEFS WITH ITS           H3                                                  
C               RECORDS HAVING THE LABEL TIME:                                  H3                                                  
C                                                                               H3                                                  
C        IBUF(1) =   4HTIME                                                     H3                                                  
C        IBUF(2) =  IF IKIND = 1, 2 OR 9, (GRID, ZONL OR SUBA) THEN IT'S THE    H3                                                  
C                                           LOCATION POINT IN LATLON FORMAT.    H3                                                  
C                   IF IKIND = 3 (COEF) THEN IT'S THE LOCATION OF COEFFICIENT   H3                                                  
C                                         IN THE REAL LA*2 LENGTH VECTOR.       H3                                                  
C        IBUF(3) =  NAME .                                                      H3                                                  
C        IBUF(4) =  LEVEL .                                                     H3                                                  
C        IBUF(5) =  2*(LENGTH-1)  WHERE LENGTH IS THE LENGTH OF THE FOURIER     H3                                                  
C                   COEFFICIENTS SERIES.  THE LONGEST TIME  SERIES THAT CAN     H3                                                  
C                   BE PROCESSED IS $TSL$ POINTS LONG.                          H3                                                  
C        IBUF(6) =  1                                                           H3                                                  
C        IBUF(7) =  DIMENSION OF GRID, KHEM AND IKIND (IN FORMAT CCCRRRKQ) OR   H3                                                  
C                    LRLMTQ, WHERE Q IS IKIND.                                  H3                                                  
C        IBUF(8) =  PACKING DENSITY                                             H3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               H5                                                  
C        LENGTH = THE ACTUAL  LENGTH USED IN  THE  TRANSFORM. THIS NUMBER       H5                                                  
C                 SHOULD  BE EVEN.  IF IT IS'NT  THE PROGRAM WILL DISCARD       H5                                                  
C                 ONE POINT. THE REST (IF ANY) OF THE SERIES IS DICARDED.       H5                                                  
C                                                                               H5                                                  
C          NOTE - IF THE INPUT PARAMETER LENGTH IS EXPRESSIBLE AS A PRODUCT     H5                                                  
C          ====   OF  2, 3, 4, 5, 6  AND/OR 8'S, THIS PROGRAM USES SUBROUTINE   H5                                                  
C                 VFFT TO COMPUTE THE FFT'S OF THE SERIES.                      H5                                                  
C                                                                               H5                                                  
C               - IF ANY OTHER PRIME NUMBER COME INTO THE DECOMPOSITION OF      H5                                                  
C                 LENGTH, ROUTINE FOURT IS USED TO DO THE FFT'S.                H5                                                  
C                                                                               H5                                                  
C               - IF THE INPUT VALUE OF LENGTH IS ZERO, THE PROGRAM WILL USE    H5                                                  
C                 THE ACTUAL LENGTH OF THE FIRST INPUT SERIES TO DETERMINE      H5                                                  
C                 WHICH OF THE TWO PRECEEDING CASES HOLD.                       H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARD...                                                       H5                                                  
C                                                                               H5                                                  
C*FCATS    9732                                                                 H5                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
                                                                                                                                    
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      DIMENSION X(129696),IX(512),WORK(129697),TRIGS(18528),                                                                        
     1          IFAX(100),IBUF(8),IDAT(18528),LFAX(7)                                                                               
      INTEGER   ADD,BASE,RRR                                                                                                        
      LOGICAL   OK                                                                                                                  
      COMMON/ICOM/IBUF,IDAT                                                                                                         
      DATA      MAXX  /18528/, MAXL /18528/, MAXW /129696/,                                                                         
     1          MAXS /64/, LFAX / 6,8,5,4,3,2,1 /                                                                                   
C     DATA      NFOUR/4HFOUR/,NTIME/4HTIME/                                                                                         
                                                                                                                                    
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     ******************************************************************                                                            
C     ****  IF INPUT TYPE OF FCOEFS IS 4HSPEC THE DATA PACKING HAS  ****                                                            
C     ****  TO BE OFF OR THE FIRST LM IMAGINARY COMPONENTS WILL BE  ****                                                            
C     ****  LOST. IN THIS CASE THE PROGRAM WILL PRINT A WARNING ... ****                                                            
C     ******************************************************************                                                            
                                                                                                                                    
C-----------------------------------------------------------------------------                                                      
      NFOUR=NC4TO8("FOUR")                                                                                                          
      NTIME=NC4TO8("TIME")                                                                                                          
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.0)                                         THEN                                                                 
          WRITE(6,6005) LENGTH                                                                                                      
          CALL                                     XIT('FCATS',-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                                                                                                                      
                                                                                                                                    
C     * THE UNPACKING OF SPECTRAL DATA RESETS CERTAIN                                                                               
C     * IMAGINARY COEFFICIENTS TO ZERO. KIND = NSPEC SHOULD                                                                         
C     * NOT BE PACKED. PROGRAM ADVISES USER OF THIS FACT...                                                                         
                                                                                                                                    
      IF (IBUF(8).NE.1 .AND. IBUF(1).EQ.NFOUR) WRITE(6,6090)                                                                        
                                                                                                                                    
C     * DOES THE INPUT HAVE THE RIGHT TYPE OF LABEL ?                                                                               
                                                                                                                                    
      IF (IBUF(6).NE.1                          .OR.                                                                                
     1   IBUF(1).NE.NFOUR )                                    THEN                                                                 
          WRITE(6,6030)                                                                                                             
          CALL                                     XIT('FCATS',-2)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      KIND = NTIME                                                                                                                  
                                                                                                                                    
C     * LENGTH HAS TO BE EVEN. IF EQUAL TO ZERO, RESET TO LEN.                                                                      
C     * IF LENGTH IS TOO SMALL FOR IBUF(5), TRUNCATE INPUT FIELD                                                                    
C     * TO ACCOMODATE IT.                                                                                                           
                                                                                                                                    
      LENGTH  = (LENGTH/2)*2                                                                                                        
                                                                                                                                    
      LENB2P1 = IBUF(5)                                                                                                             
      LEN     = 2*(LENB2P1-1)                                                                                                       
      IF (LENGTH.LT.LEN .AND. LENGTH.NE.0)                     THEN                                                                 
          LEN    = LENGTH                                                                                                           
          WRITE(6,6045) LEN/2+1                                                                                                     
      ELSE IF (LENGTH.EQ.0)                                    THEN                                                                 
          LENGTH = LEN                                                                                                              
          WRITE(6,6035) LENGTH                                                                                                      
      ENDIF                                                                                                                         
                                                                                                                                    
C     * IS LENGTH TOO BIG ?                                                                                                         
                                                                                                                                    
      IF (LENGTH.GT.MAXL)                                      THEN                                                                 
          WRITE(6,6040)                                                                                                             
          CALL                                     XIT('FCATS',-5)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      JUMP    = LENGTH+2                                                                                                            
      NSET    = -1                                                                                                                  
                                                                                                                                    
C     * TIME SERIES WILL BE PROCESSED IN GROUPS OF NBLOCK AT A                                                                      
C     * TIME TO TAKE ADVANTAGE OF THE VECTORIZATION OF THE FFT.                                                                     
                                                                                                                                    
      NBLOCK = MIN( (MAXW-2*LENB2P1)/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.NFOUR )                            THEN                                                                 
                 WRITE(6,6050) NSET*NBLOCK + N                                                                                      
                 WRITE(6,6020) IBUF                                                                                                 
                 WRITE(6,6060)                                                                                                      
                 CALL                              XIT('FCATS',-6)                                                                  
             ENDIF                                                                                                                  
                                                                                                                                    
             IF (IBUF(5).NE.LENB2P1)                           THEN                                                                 
                 WRITE(6,6050) NSET*NBLOCK + N                                                                                      
                 WRITE(6,6020) IBUF                                                                                                 
                 WRITE(6,6070)                                                                                                      
                 CALL                              XIT('FCATS',-7)                                                                  
             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 THE COEFFICIENTS WITH ZEROS.                                                                              
                                                                                                                                    
         DO 1150 N=1,NS                                                                                                             
             ADD = (N-1)*JUMP                                                                                                       
             DO 1150 I=LEN+3,JUMP                                                                                                   
                 X(ADD+I) = 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,-1,WORK,LENGTH)                                                                      
 1170        CONTINUE                                                                                                               
         ENDIF                                                                                                                      
                                                                                                                                    
C        * OUTPUT THE BLOCK                                                                                                         
                                                                                                                                    
         IBUF(1) = KIND                                                                                                             
         DO 1200 N=1,NS                                                                                                             
             BASE = (N-1)*8                                                                                                         
             DO 1190 J=2,8                                                                                                          
                 IBUF(J) = IX(BASE+J)                                                                                               
 1190        CONTINUE                                                                                                               
             IBUF(5) = LENGTH                                                                                                       
             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.                                                                                                                   
                                                                                                                                    
      IF (NS.EQ.0) NS = NBLOCK                                                                                                      
      BASE            = (NS-1)*8                                                                                                    
      IX(BASE+1)      = KIND                                                                                                        
      IX(BASE+5)      = LENGTH                                                                                                      
      WRITE(6,6020) (IX(BASE+I),I=1,8)                                                                                              
      WRITE(6,6080) NSET*NBLOCK+NS                                                                                                  
      CALL                                         XIT('FCATS',0)                                                                   
                                                                                                                                    
C     * PREMATURE EOF ON UNIT 1.                                                                                                    
                                                                                                                                    
 9000 WRITE(6,6000)                                                                                                                 
      CALL                                         XIT('FCATS',-8)                                                                  
                                                                                                                                    
C     * PREMATURE EOF ON UNIT 1.                                                                                                    
                                                                                                                                    
 9010 WRITE(6,6010)                                                                                                                 
      CALL                                         XIT('FCATS',-9)                                                                  
C-----------------------------------------------------------------------------                                                      
                                                                                                                                    
 5000 FORMAT(10X,I5)                                                            H4                                                  
 6000 FORMAT('0INPUT FILE IS EMPTY.')                                                                                               
 6005 FORMAT('0LENGTH IS ',I6)                                                                                                      
 6010 FORMAT('0INPUT FILE OF FOURIER COEFFICIENTS IS EMPTY.')                                                                       
 6020 FORMAT(1X,A4,I10,1X,A4,5I10)                                                                                                  
 6030 FORMAT('0FIRST RECORD IN INPUT FILE OF FOURIER COEFFICIENTS',                                                                 
     1       ' DOES NOT HAVE A FOURIER LIKE LABEL.')                                                                                
 6035 FORMAT('0LENGTH SET TO 2*(IBUF(5)-1) =',I5)                                                                                   
 6040 FORMAT('0A RECORD IN THE INPUT FILE IS TOO LONG.')                                                                            
 6045 FORMAT('0ONLY ',I5,' COEFFICIENTS WILL BE USED IN EXPANSION.')                                                                
 6050 FORMAT('0LABEL FOR RECORD ',I5,':')                                                                                           
 6060 FORMAT(11X,'IS NOT A FOURIER COEFFICIENTS LABEL.')                                                                            
 6070 FORMAT('0CHANGE OF LENGTH DETECTED IN THE INPUT SERIES.')                                                                     
 6080 FORMAT('0PROCESSED ',I10,' RECORDS.')                                                                                         
 6090 FORMAT('0INPUT SPECTRAL FIELD IS PACKED AND WILL LOOSE DATA')                                                                 
      END                                                                                                                           
