      PROGRAM MKEOF                                                                                                                 
C     PROGRAM MKEOF (X,       XEOFS,       INPUT,       OUTPUT,         )       H2                                                  
C    1         TAPE1=X, TAPE2=XEOFS, TAPE5=INPUT, TAPE6=OUTPUT)                                                                     
C     ---------------------------------------------------------                 H2                                                  
C                                                                               H2                                                  
C     AUG 20/08 - S.KHARIN (FIX A BUG IN CALCULATION OF THE COVARIANCE MATRIX   H2                                                  
C                           WHEN NTIME>LENGTH.                                  H2                                                  
C                           SET THE SIGN OF EOFS TO BE MOSTLY POSITIVE.)        H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     APR 03/03 - F.MAJAESS (EMPLOY CLOSE/OPEN STATEMENTS TO ACCOMPLISH PROPER                                                      
C                            FILE POINTER POSITIONING)                                                                              
C     NOV 12/02 - SLAVA KHARIN (USE "1" FOR TIMESTEP IN WRITING OUT "TAVG" AND                                                      
C                               "VARI" RECORDS INSTEAD OF "NSTEP")                                                                  
C     APR 12/99 - SLAVA KHARIN (LIMIT PACKING DENSITY WARNING MESSAGE TO 1)                                                         
C     MAR 09/99 - SLAVA KHARIN.                                                                                                     
C                                                                               H2                                                  
CMKEOF  - COMPUTES EOFS AND PCS.                                        1  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - SLAVA KHARIN                                                         H3                                                  
C                                                                               H3                                                  
CPURPOSE - COMPUTES EMPIRICAL ORTHOGONAL FUNCTIONS AND PRINCIPAL COMPONENTS.    H3                                                  
C                                                                               H3                                                  
CINPUT FILE...                                                                  H3                                                  
C                                                                               H3                                                  
C     X = INPUT FILE WITH TIME SERIES OF FIELDS IN CCRN FORMAT.                 H3                                                  
C         THE TIME MEAN WILL BE REMOVED FROM THE FIELDS.                        H3                                                  
C         IF THE LAST LONGITUDE IS THE SAME AS THE FIRST ONE IN ALL             H3                                                  
C         INPUT FIELDS, THE PROGRAM ASSUMES THAT THIS LONGITUDE IS ADDED        H3                                                  
C         FOR CYCLICITY AND IT IS IGNORED IN CALCULATION OF THE COVARIANCE/     H3                                                  
C         CORRELATION MATRIX.                                                   H3                                                  
C                                                                               H3                                                  
COUTPUT FILE...                                                                 H3                                                  
C                                                                               H3                                                  
C       XEOFS = OUPUT FILE WITH EVS, EOFS AND PCS ETC. IN CCRN FORMAT           H3                                                  
C               1ST RECORD = EIGENVALUES               (NAME = 4HEVAL)          H3                                                  
C               2ND RECORD = PERCENTAGES OF EXPL. VAR. (NAME = 4H PCT)          H3                                                  
C               3RD RECORD = CUMULATIVE PERCENTAGES    (NAME = 4HPCTC)          H3                                                  
C               4TH RECORD = THE TIME MEAN             (NAME = 4HTAVG)          H3                                                  
C               5TH RECORD = THE VARIANCE              (NAME = 4HVARI)          H3                                                  
C               6TH RECORD = 1ST EOF                   (NAME = 4HEVEC)          H3                                                  
C               7TH RECORD = 1ST PC                    (NAME = 4H  PC)          H3                                                  
C               8TH RECORD = 2ND EOF                   (NAME = 4HEVEC)          H3                                                  
C               9TH RECORD = 2ND PC                    (NAME = 4H  PC)          H3                                                  
C               ... ETC.                                                        H3                                                  
C                                                                               H3                                                  
C   NOTE - IF YOU USE AREA WEIGHTING OF INPUT FIELDS BEFORE RUNNING MKEOF,      H3                                                  
C          DO NOT FORGET TO REMOVE THIS WEIGHTING FROM THE EOFS (EVEC),         H3                                                  
C          TIME MEAN (TAVG) AND VARIANCE (VARI) AFTERWARDS.                     H3                                                  
C                                                                                                                                   
CINPUT PARAMETERS:                                                                                                                  
C                                                                               H5                                                  
C     NEOF  =   THE NUMBER OF REQUIERED EOFS.                                   H5                                                  
C               IF NEOF=0 ALL EOFS ARE WRITTEN OUT.                             H5                                                  
C     NORM  = 0 EOFS ARE NORMALIZED BY 1 (|EOF|^2=1) AND VAR(PC)=EVAL.          H5                                                  
C               THIS IS THE STANDARD DEFINITION OF EOFS AND PCS.                H5                                                  
C           = 1 THE PRINCIPAL COMPONENTS ARE STANDARDIZED (VAR(PC)=1) AND       H5                                                  
C               EOFS ARE NORMALIZED BY THE SQUARE ROOT OF THE EXPLAINED         H5                                                  
C               VARIANCE (|EOF|^2=EVAL).                                        H5                                                  
C     ICOV  = 0 EIGENVALUE PROBLEM IS SOLVED FOR THE COVARIANCE MATRIX          H5                                                  
C               (THIS IS STANDARD DEFINITION).                                  H5                                                  
C           = 1 EIGENVALUE PROBLEM IS SOLVED FOR THE CORRELATION MATRIX.        H5                                                  
C     ITRAN = 0 TRANSPOSE INPUT DATA IN TIME/SPACE ONLY IF THE NUMBER OF        H5                                                  
C               TIME STEPS IS SMALLER THAN THE LENGTH OF THE VECTORS.           H5                                                  
C               (RECOMMENDED).                                                  H5                                                  
C           = 1 DO NOT TRANSPOSE INPUT DATA.                                    H5                                                  
C           =-1 DO TRANSPOSE THE INPUT DATA.                                    H5                                                  
C               NOTE: THE TIME/SPACE TRANSPOSED DATA ARE WRITTEN TO             H5                                                  
C               FORTRAN UNIT 77                                                 H5                                                  
C     SPVAL  =  SPECIAL VALUE TO INDICATE MISSING VALUES IN INPUT FIELDS.       H5                                                  
C     PCTMIN =  MINIMAL PERCENTAGE (%) OF NON-MISSING VALUES AT A GRID BOX      H5                                                  
C               FOR ESTIMATING THE COVARIANCE MATRIX.                           H5                                                  
C               IF THE PERCENTAGE OF NON-MISSING VALUES AT ANY PARTICULAR GRID  H5                                                  
C               POINT IS LESS THAN PCTMIN, THIS GRID POINT IS IGNORED           H5                                                  
C               (I.E. IT IS ASSUMED THAT ALL VALUES ARE MISSING).               H5                                                  
C               PCTMIN=0 IS EQUIVALENT TO PCTMIN=100.                           H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARD...                                                       H5                                                  
C                                                                               H5                                                  
C*MKEOF       5    0    0    0    1.E+38 75.0                                   H5                                                  
C                                                                               H5                                                  
C (COMPUTE THE FIRST 5 EOFS WHICH ARE NORMALIZED ONE.                           H5                                                  
C  VAR(PCS) = EXPLAINED VARIANCE.                                               H5                                                  
C  COVARIANCE MATRIX IS USED.                                                   H5                                                  
C  MISSING VALUE IS 1.E+38.                                                     H5                                                  
C  IF THE PERCENTAGE OF NON-MISSING DATA AT A GRID POINT                        H5                                                  
C  IS LESS THAN 75% THAN THIS GRID POINT IS IGNORED.                            H5                                                  
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
C                                                                                                                                   
C     * LENMAX = MAX LENGTH OF A DATA RECORD.                                                                                       
C     * NTMAX  = MAX NUMBER OF TIME STEPS.                                                                                          
C     * NMAX   = MAX DIMENSION OF THE COVARIANCE/CORRELATION MATRIX.                                                                
C                EITHER THE NUMBER OF INPUT FIELDS OR THE LENGTH OF FIELDS                                                          
C                MUST BE SMALLER THAN OR EQUAL TO NMAX.                                                                             
C                12012=1001 YRS * 12 MONTHS                                                                                         
C     * NEOFOUT= NUMBER OF EOFS FOR WHICH EXPLAINED VARIANCES ARE PRINTED OUT.                                                      
C     * NBUF   = SIZE OF BUFFER USED FOR TIME/SPACE TRANSFORMATION                                                                  
C                                                                                                                                   
      PARAMETER (NTMAX=18528,LENMAX=65341,NMAX=12012,NEOFOUT=50,                                                                    
     1     NBUF=NMAX*(NMAX+4)+NTMAX+LENMAX)                                                                                         
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
      COMMON/JCOM/JBUF(8),JDAT(18528)                                                                                               
C                                                                                                                                   
C     * ARRAY FOR INPUT FIELDS                                                                                                      
C                                                                                                                                   
      REAL X(LENMAX)                                                                                                                
C                                                                                                                                   
C     * ARRAYS FOR OUTPUT FIELDS                                                                                                    
C                                                                                                                                   
      REAL  EVAL(NMAX),EVAL0(NMAX),PCT(NMAX),PCTC(NMAX)                                                                             
     1     ,AVG(LENMAX),VAR(LENMAX),EVEC(LENMAX),PC(NTMAX)                                                                          
C                                                                                                                                   
C     * INTERNAL ARRAYS:                                                                                                            
C     * COV IS COVARIANCE MATRIX                                                                                                    
C     * DATBUF IS USED FOR TIME/SPACE TRANSFORMATION OF INPUT FIELDS                                                                
C                                                                                                                                   
      REAL  COV(NMAX,NMAX),DATBUF(NBUF)                                                                                             
C                                                                                                                                   
C     * SAVE SPACE BY MAKING EQUIVALENCE                                                                                            
C                                                                                                                                   
      COMMON /MKEOFBUF/ COV,EVAL,EVAL0,PCT,PCTC,EVEC,PC                                                                             
      EQUIVALENCE (DATBUF(1),COV(1,1))                                                                                              
C                                                                                                                                   
      INTEGER MISSL(LENMAX),MISST(NTMAX),IVAL(LENMAX),INDX(NMAX)                                                                    
C                                                                                                                                   
      DATA MAXX/65341/,MAXT/18528/                                                                                                  
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
C     * ASSIGN FILES TO FORTRAN UNITS                                                                                               
C                                                                                                                                   
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,1,2,5,6)                                                                                                      
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      READ(5,5010,END=900)NEOF,NORM,ICOV,ITRAN,SPVAL,PCTMIN                     H4                                                  
      SPVALT=1.E-6*MAX(1.E0,ABS(SPVAL))                                                                                             
C                                                                                                                                   
C     * PRINT THE INPUT PARAMETERS TO STANDARD OUTPUT                                                                               
C                                                                                                                                   
      WRITE (6,6010) NEOF,NORM,ICOV,ITRAN,SPVAL,PCTMIN                                                                              
                                                                                                                                    
      IF (NEOF.EQ.0) THEN                                                                                                           
        WRITE(6,6020)                                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (NORM.EQ.0) THEN                                                                                                           
        WRITE(6,6030)                                                                                                               
      ELSE                                                                                                                          
        WRITE(6,6035)                                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (ICOV.EQ.0) THEN                                                                                                           
        WRITE(6,6040)                                                                                                               
      ELSE                                                                                                                          
        WRITE(6,6045)                                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (ITRAN.EQ.0) THEN                                                                                                          
        WRITE(6,6050)                                                                                                               
      ELSE IF (ITRAN.EQ.1) THEN                                                                                                     
        WRITE(6,6055)                                                                                                               
      ELSE IF (ITRAN.EQ.-1) THEN                                                                                                    
        WRITE(6,6060)                                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (PCTMIN.EQ.0.E0) PCTMIN=100.E0                                                                                             
C                                                                                                                                   
C     * INITIALIZE SOME ARRAYS                                                                                                      
C                                                                                                                                   
      DO IJ=1,LENMAX                                                                                                                
        AVG(IJ)=0.E0                                                                                                                
        VAR(IJ)=0.E0                                                                                                                
C                                                                                                                                   
C       * MISST=NUMBER OF MISSING VALUES AT EACH GRID POINT                                                                         
C       * IVAL=GRID POINTS THAT SATISFY THE PCTMIN CRITERIUM.                                                                       
C                                                                                                                                   
        MISSL(IJ)=0                                                                                                                 
        IVAL(IJ)=0                                                                                                                  
      ENDDO                                                                                                                         
      DO NT=1,NTMAX                                                                                                                 
C                                                                                                                                   
C       * MISST=NUMBER OF MISSING VALUES AT EACH TIME STEP                                                                          
C                                                                                                                                   
        MISST(NT)=0                                                                                                                 
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * COUNT THE NUMBER OF TIME STEPS AND VALID GRID POINTS IN THE                                                                 
C     * INPUT FILE                                                                                                                  
C                                                                                                                                   
      NSTEP=0                                                                                                                       
      ICYCLE=1                                                                                                                      
 100  CONTINUE                                                                                                                      
      CALL GETFLD2(1,X,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                    
      IF(.NOT.OK) THEN                                                                                                              
        IF(NSTEP.EQ.0) CALL                        XIT('MKEOF',-1)                                                                  
        GO TO 300                                                                                                                   
      ENDIF                                                                                                                         
      NSTEP=NSTEP+1                                                                                                                 
C                                                                                                                                   
C     * CHECK WHETHER THE LAST MERIDIAN IS THE SAME AS THE FIRST ONE                                                                
C     * AND SET ICYCLE TO 1 IF IT IS THE CASE.                                                                                      
C     * (DO NOT CHECK IF ICYCLE IS ALREADY EQUAL TO 0)                                                                              
C                                                                                                                                   
      KIND=IBUF(1)                                                                                                                  
      LEV=IBUF(4)                                                                                                                   
      NLG=IBUF(5)                                                                                                                   
      NLAT=IBUF(6)                                                                                                                  
      KHEM=IBUF(7)                                                                                                                  
      NPCK=IBUF(8)                                                                                                                  
      NLON=NLG                                                                                                                      
C                                                                                                                                   
C     * GIVE A WARNING IF NPCK > 1                                                                                                  
C                                                                                                                                   
      IF (NPCK.GT.1.AND.NSTEP.EQ.1) THEN                                                                                            
        WRITE(6,'(3A)')                                                                                                             
     1       ' *** WARNING: PACKING DENSITY > 1.',                                                                                  
     2       ' MAKE SURE INPUT DATA DOES NOT CONTAIN VERY BIG SPECIAL',                                                             
     3       ' VALUE, E.G. 1.E38 ***'                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (KIND.EQ.NC4TO8("SPEC").OR.                                                                                                
     +    KIND.EQ.NC4TO8("FOUR")) NLON=NLG*2                                                                                        
      IF (NLON*NLAT.GT.LENMAX) THEN                                                                                                 
        WRITE(6,'(A)')                                                                                                              
     1        ' ***ERROR: NLEN>LENMAX. INCREASE LENMAX IN MKEOF.F ***'                                                              
        CALL                                       XIT('MKEOF',-2)                                                                  
      ENDIF                                                                                                                         
      IF (KIND.NE.NC4TO8("GRID")) ICYCLE=0                                                                                          
      IF (ICYCLE.EQ.0) GO TO 220                                                                                                    
C                                                                                                                                   
C     * TEST FOR CYCLIC MERIDIAN                                                                                                    
C                                                                                                                                   
      SUM=0.E0                                                                                                                      
      DO J=1,NLAT                                                                                                                   
        SUM=SUM+ABS(X((J-1)*NLON+1)-X(J*NLON))                                                                                      
      ENDDO                                                                                                                         
      IF (SUM/FLOAT(NLAT) .GT. 1.E-12) ICYCLE=0                                                                                     
 220  CONTINUE                                                                                                                      
      NLON1=NLG-ICYCLE                                                                                                              
      IF (KIND.EQ.NC4TO8("SPEC").OR.                                                                                                
     +    KIND.EQ.NC4TO8("FOUR")) NLON1=2*NLON1                                                                                     
      NWDS=NLON1*NLAT                                                                                                               
C                                                                                                                                   
C     * COUNT MISSING VALUES AND CALCULATE TIME MEAN AND VARIANCE                                                                   
C                                                                                                                                   
      DO J=1,NLAT                                                                                                                   
        DO I=1,NLON1                                                                                                                
          IJ=(J-1)*NLON+I                                                                                                           
          IF (ABS(X(IJ)-SPVAL).LT.SPVALT) THEN                                                                                      
            MISSL(IJ)=MISSL(IJ)+1                                                                                                   
            MISST(NSTEP)=MISST(NSTEP)+1                                                                                             
          ELSE                                                                                                                      
            AVG(IJ)=AVG(IJ)+X(IJ)                                                                                                   
            VAR(IJ)=VAR(IJ)+X(IJ)**2                                                                                                
          ENDIF                                                                                                                     
        ENDDO                                                                                                                       
      ENDDO                                                                                                                         
      GO TO 100                                                                                                                     
                                                                                                                                    
 300  CONTINUE                                                                                                                      
      IF (ICYCLE.EQ.1) THEN                                                                                                         
        WRITE(6,6065)                                                                                                               
      ENDIF                                                                                                                         
      WRITE(6,6070) NSTEP                                                                                                           
      IF (NSTEP.GT.NTMAX) THEN                                                                                                      
        WRITE(6,'(A)')                                                                                                              
     1       ' *** ERROR: NSTEP>NTMAX. INCREASE NTMAX IN MKEOF.F ***'                                                               
        CALL                                       XIT('MKEOF',-3)                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * APPLY CRITERIA FOR VALID POINTS                                                                                             
C                                                                                                                                   
      NLEN=0                                                                                                                        
      NMISS=0                                                                                                                       
      DO J=1,NLAT                                                                                                                   
        DO I=1,NLON1                                                                                                                
          IJ=(J-1)*NLON+I                                                                                                           
          IF ((1.E0-FLOAT(MISSL(IJ))/FLOAT(NSTEP))*100.E0                                                                           
     &         .LT.PCTMIN) THEN                                                                                                     
            NMISS=NMISS+1                                                                                                           
          ELSE                                                                                                                      
            NLEN=NLEN+1                                                                                                             
            IVAL(NLEN)=IJ                                                                                                           
          ENDIF                                                                                                                     
        ENDDO                                                                                                                       
      ENDDO                                                                                                                         
      WRITE(6,6075) NLEN                                                                                                            
      WRITE(6,6080) NMISS                                                                                                           
C                                                                                                                                   
C     * CHECK THE DIMENSION OF MATRIX                                                                                               
C                                                                                                                                   
      IF (ITRAN.EQ.0) THEN                                                                                                          
        NMIN=MIN(NLEN,NSTEP)                                                                                                        
      ELSE IF (ITRAN.EQ.1) THEN                                                                                                     
        NMIN=NLEN                                                                                                                   
      ELSE IF (ITRAN.EQ.-1) THEN                                                                                                    
        NMIN=NSTEP                                                                                                                  
      ENDIF                                                                                                                         
      IF (NMAX.LT.NMIN) THEN                                                                                                        
        WRITE (6,'(2A,I6)')                                                                                                         
     1       ' ***ERROR: THE COVARIANCE MATRIX IS TOO LARGE.',                                                                      
     2       ' INCREASE NMAX TO AT LEAST',NMIN                                                                                      
        CALL                                       XIT('MKEOF',-4)                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPLETE CALCULATION OF THE TIME MEAN AND VARIANCE                                                                          
C     * (DO IT ONLY FOR VALID GRID POINTS)                                                                                          
C                                                                                                                                   
      DO IV=1,NLEN                                                                                                                  
        AVG(IVAL(IV))=AVG(IVAL(IV))/FLOAT(NSTEP-MISSL(IVAL(IV)))                                                                    
        VAR(IVAL(IV))=MAX(0.E0,                                                                                                     
     1                VAR(IVAL(IV))/FLOAT(NSTEP-MISSL(IVAL(IV)))-                                                                   
     2                AVG(IVAL(IV))**2)                                                                                             
      ENDDO                                                                                                                         
                                                                                                                                    
      IF (ITRAN.EQ.1.OR.(ITRAN.EQ.0.AND.NSTEP.GE.NLEN)) THEN                                                                        
C                                                                                                                                   
C       * THE INPUT FIELDS ARE NOT TIME/SPACE TRANSPOSED.                                                                           
C                                                                                                                                   
        WRITE(6,'(A)')' CALCULATE THE DIRECT COVARIANCE MATRIX.'                                                                    
        LTRANS=0                                                                                                                    
        NMAT=NLEN                                                                                                                   
C                                                                                                                                   
C       * INITIALIZE THE COVARIANCE/CORRELATION MATRIX                                                                              
C                                                                                                                                   
        DO N2=1,NLEN                                                                                                                
          DO N1=1,NLEN                                                                                                              
            COV(N1,N2)=0.E0                                                                                                         
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * ACCUMULATE COVARIANCE/CORRELATION MATRIX COEFFICIENTS                                                                     
C                                                                                                                                   
        REWIND 1                                                                                                                    
        DO N=1,NSTEP                                                                                                                
          CALL GETFLD2(1,X,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                
          IF(.NOT.OK)THEN                                                                                                           
            CALL                                   XIT('MKEOF',-5)                                                                  
          ENDIF                                                                                                                     
C                                                                                                                                   
C         * FIND VALID VALUES AND SUBTRACT TIME MEAN                                                                                
C         * SET MISSING VALUES TO 0.                                                                                                
C                                                                                                                                   
          DO IV=1,NLEN                                                                                                              
            IF(ABS(X(IVAL(IV))-SPVAL).GT.SPVALT) THEN                                                                               
              X(IVAL(IV))=X(IVAL(IV))-AVG(IVAL(IV))                                                                                 
            ELSE                                                                                                                    
              X(IVAL(IV))=0.E0                                                                                                      
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * UPPER TRIANGLE                                                                                                          
C                                                                                                                                   
          DO IV2=1,NLEN                                                                                                             
            DO IV1=IV2,NLEN                                                                                                         
              COV(IV1,IV2)=COV(IV1,IV2)+X(IVAL(IV1))*X(IVAL(IV2))                                                                   
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * CALCULATION OF CORRELATION MATRIX                                                                                         
C                                                                                                                                   
        IF (ICOV.NE.0) THEN                                                                                                         
          WRITE(6,'(A)')' CALCULATION OF THE CORRELATION MATRIX.'                                                                   
          DO IV2=1,NLEN                                                                                                             
            DO IV1=IV2,NLEN                                                                                                         
              COV(IV1,IV2)=COV(IV1,IV2)/                                                                                            
     1             SQRT(VAR(IVAL(IV1))*VAR(IVAL(IV2)))                                                                              
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * COMPLETE CALCULATION OF COVARIANCE/CORRELATION MATRIX                                                                     
C                                                                                                                                   
        DO IV2=1,NLEN                                                                                                               
          DO IV1=IV2,NLEN                                                                                                           
            COV(IV1,IV2)=COV(IV1,IV2)/                                                                                              
     1           FLOAT(NSTEP-MAX(MISSL(IVAL(IV1)),MISSL(IVAL(IV2))))                                                                
            COV(IV2,IV1)=COV(IV1,IV2)                                                                                               
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
      ELSE                                                                                                                          
C                                                                                                                                   
C       * THE COVARIANCE/CORRELATION MATRIX WILL BE CALCULATED FOR                                                                  
C       * TRANSPOSED INPUT DATA.                                                                                                    
C                                                                                                                                   
C       * CALCULATE THE PORTION OF THE INPUT FIELDS (ALL TIME STEPS)                                                                
C       * FITTING IN THE BUFFER DATBUF.                                                                                             
C                                                                                                                                   
        WRITE(6,'(A)')' CALCULATE THE TRANSPOSED COVARIANCE MATRIX.'                                                                
        LTRANS=1                                                                                                                    
        NMAT=NSTEP                                                                                                                  
        NPOINTS=NBUF/NSTEP                                                                                                          
        NREAD=(NLEN-1)/NPOINTS+1                                                                                                    
C                                                                                                                                   
C       * READ THE INPUT FILE NREAD TIMES AND WRITE OUT THE TRANSPOSED FILE                                                         
C       * TO FORTRAN UNIT 77                                                                                                        
C                                                                                                                                   
        OPEN(77,IOSTAT=IRTC,FILE='XXXMKEOF',FORM='UNFORMATTED')                                                                     
        IF (IRTC.NE.0) CALL                        XIT('MKEOF',-6)                                                                  
        REWIND 77                                                                                                                   
        DO NR=1,NREAD                                                                                                               
          REWIND 1                                                                                                                  
          IJ=0                                                                                                                      
          DO N=1,NSTEP                                                                                                              
            CALL GETFLD2(1,X,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                              
            IF(.NOT.OK)THEN                                                                                                         
              CALL                                 XIT('MKEOF',-7)                                                                  
            ENDIF                                                                                                                   
            DO IV=(NR-1)*NPOINTS+1,MIN(NR*NPOINTS,NLEN)                                                                             
              IJ=IJ+1                                                                                                               
              DATBUF(IJ)=X(IVAL(IV))                                                                                                
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
          NP=IJ/NSTEP                                                                                                               
          DO I=1,NP                                                                                                                 
            WRITE(77)(DATBUF((N-1)*NP+I),N=1,NSTEP)                                                                                 
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * INITIALIZE THE COVARIANCE/CORRELATION MATRIX                                                                              
C                                                                                                                                   
        DO N2=1,NSTEP                                                                                                               
          DO N1=1,NSTEP                                                                                                             
            COV(N1,N2)=0.E0                                                                                                         
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * CALCULATE THE COVARIANCE/CORRELATION MATRIX FOR TIME/SPACE                                                                
C       * TRANSPOSED INPUT DATA BY ACCUMULATING THE MATRIX COEFFICIENTS                                                             
C                                                                                                                                   
        CLOSE (77)                                                                                                                  
        OPEN(77,IOSTAT=IRTC,FILE='XXXMKEOF',FORM='UNFORMATTED')                                                                     
        REWIND 77                                                                                                                   
        DO IV=1,NLEN                                                                                                                
          READ (77)(X(N),N=1,NSTEP)                                                                                                 
C                                                                                                                                   
C         * FIND VALID VALUES AND SUBTRACT MEAN                                                                                     
C         * SET MISSING VALUES TO 0.                                                                                                
C                                                                                                                                   
          DO N=1,NSTEP                                                                                                              
            IF (ABS(X(N)-SPVAL).GT.SPVALT) THEN                                                                                     
              X(N)=X(N)-AVG(IVAL(IV))                                                                                               
            ELSE                                                                                                                    
              X(N)=0.E0                                                                                                             
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * UPPER TRIANGLE                                                                                                          
C                                                                                                                                   
          IF (ICOV.EQ.0) THEN                                                                                                       
C                                                                                                                                   
C           * COVARIANCE MATRIX                                                                                                     
C                                                                                                                                   
            DO N2=1,NSTEP                                                                                                           
              DO N1=N2,NSTEP                                                                                                        
                COV(N1,N2)=COV(N1,N2)+X(N1)*X(N2)                                                                                   
              ENDDO                                                                                                                 
            ENDDO                                                                                                                   
          ELSE                                                                                                                      
C                                                                                                                                   
C           * CORRELATION MATRIX                                                                                                    
C                                                                                                                                   
            DO N2=1,NSTEP                                                                                                           
              DO N1=N2+1,NSTEP                                                                                                      
                COV(N1,N2)=COV(N1,N2)+X(N1)*X(N2)/VAR(IVAL(IV))                                                                     
              ENDDO                                                                                                                 
            ENDDO                                                                                                                   
          ENDIF                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * COMPLETE  CALCULATION OF COVARIANCE/CORRELATION MATRIX                                                                    
C                                                                                                                                   
        DO N2=1,NSTEP                                                                                                               
          COV(N2,N2)=COV(N2,N2)/FLOAT(NWDS-MISST(N2))                                                                               
          DO N1=N2+1,NSTEP                                                                                                          
            COV(N1,N2)=COV(N1,N2)/                                                                                                  
     1           FLOAT(NWDS-MAX(MISST(N1),MISST(N2)))                                                                               
            COV(N2,N1)=COV(N1,N2)                                                                                                   
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
      WRITE(6,'(A)')' COVARIANCE MATRIX IS COMPUTED.'                                                                               
C                                                                                                                                   
C     * THE COVARIANCE/CORRELATION MATRIX IS READY NOW.                                                                             
C     * WE USE SUBROUTINES TRED2 AND TQLI1 FROM NUMERICAL RECIPES TO SOLVE                                                          
C     * THE EIGENVALUE PROBLEM.                                                                                                     
C                                                                                                                                   
      CALL TRED2(COV,NMAT,NMAX,EVAL,PCT)                                                                                            
      WRITE(6,'(A)')' HOUSEHOLDER REDUCTION IS COMPLETED.'                                                                          
      CALL TQLI1(EVAL,PCT,NMAT,NMAX,COV)                                                                                            
      WRITE(6,'(A)')' QL DECOMPOSITION IS COMPLETED.'                                                                               
C                                                                                                                                   
C     * THE EIGENVALUES ARE IN EVAL AND EIGENVECTORS ARE IN COV                                                                     
C                                                                                                                                   
      IF (NEOF.EQ.0) NEOF=NMAT                                                                                                      
      NEOF=MIN(NEOF,NLEN,NSTEP-1)                                                                                                   
      WRITE(6,6090) NEOF                                                                                                            
C                                                                                                                                   
C     * RENORMALIZE EIGENVALUES IF THEY ARE COMPUTED FOR TRANSPOSED DATA                                                            
C                                                                                                                                   
      IF (LTRANS.EQ.1) THEN                                                                                                         
        DO NE=1,NMAT                                                                                                                
          EVAL(NE)=EVAL(NE)*FLOAT(NLEN)/FLOAT(NSTEP)                                                                                
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * INDEX THE EIGENVALUES FROM MIN TO MAX                                                                                       
C                                                                                                                                   
      CALL INDEXX(NMAT,EVAL,INDX)                                                                                                   
C                                                                                                                                   
C     * REVERT FROM MAX TO MIN                                                                                                      
C                                                                                                                                   
      DO NE=1,NMAT/2                                                                                                                
        NE1=INDX(NE)                                                                                                                
        INDX(NE)=INDX(NMAT-NE+1)                                                                                                    
        INDX(NMAT-NE+1)=NE1                                                                                                         
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * COMPUTE PERCENTAGE OF EXPLAINED VARIANCE                                                                                    
C                                                                                                                                   
      SUM=0.E0                                                                                                                      
      DO NE=1,NMAT                                                                                                                  
C                                                                                                                                   
C     * SOME EIGENVALUES ARE NEGATIVE DUE TO NUMERICAL ROUND-OFF.                                                                   
C     * REPLACE WITH ZEROS.                                                                                                         
C                                                                                                                                   
        IF(EVAL(NE).LT.0.E0)EVAL(NE)=0.E0                                                                                           
        SUM=SUM+EVAL(NE)                                                                                                            
      ENDDO                                                                                                                         
      DO NE=1,NEOF                                                                                                                  
        PCT(NE)=EVAL(INDX(NE))/SUM                                                                                                  
      ENDDO                                                                                                                         
      PCTC(1)=PCT(1)                                                                                                                
      DO NE=2,NEOF                                                                                                                  
        PCTC(NE)=PCTC(NE-1)+PCT(NE)                                                                                                 
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * PRINT EIGENVALUES AND EXPLAINED VARIANCES                                                                                   
C                                                                                                                                   
      WRITE(6,6100) SUM, NEOF                                                                                                       
      DO NE=1,MIN(NEOF,NEOFOUT)                                                                                                     
        WRITE(6,6105)NE,EVAL(INDX(NE)),100.E0*PCT(NE),100.E0*PCTC(NE)                                                               
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * WE PRINT OUT EVERY 10TH EV FOR HIGHER EOF NUMBERS                                                                           
C                                                                                                                                   
      DO NE=NEOFOUT+10,NEOF,10                                                                                                      
        WRITE(6,6105)NE,EVAL(INDX(NE)),100.E0*PCT(NE),100.E0*PCTC(NE)                                                               
      ENDDO                                                                                                                         
      IF (ICOV.EQ.1) WRITE(6,6110)                                                                                                  
C                                                                                                                                   
C     * SAVE THE RESULTS                                                                                                            
C                                                                                                                                   
C     * THE FIRST RECORD IN THE OUTPUT FILE CONTAINS THE EIGENVALUES                                                                
C     * (NAME=4HEVAL)                                                                                                               
C                                                                                                                                   
      DO NE=1,NEOF                                                                                                                  
        EVAL0(NE)=EVAL(INDX(NE))                                                                                                    
      ENDDO                                                                                                                         
      CALL SETLAB(IBUF,NC4TO8("TIME"),1,NC4TO8("EVAL"),LEV,                                                                         
     +                                       NEOF,1,KHEM,1)                                                                         
      CALL PUTFLD2(2,EVAL0,IBUF,MAXX)                                                                                               
C                                                                                                                                   
C     * THE SECOND RECORD IS THE PERCENTAGES OF THE TOTAL VARIANCE                                                                  
C     * EXPLAINED BY EACH COMPONENT (NAME=4H PCT)                                                                                   
C                                                                                                                                   
      CALL SETLAB(IBUF,NC4TO8("TIME"),1,NC4TO8(" PCT"),LEV,                                                                         
     +                                       NEOF,1,KHEM,1)                                                                         
      CALL PUTFLD2(2,PCT,IBUF,MAXX)                                                                                                 
C                                                                                                                                   
C     * THE THIRD RECORD IS THE CUMULATIVE PERCENTAGES OF THE TOTAL                                                                 
C     * VARIANCE (NAME=4HPCTC)                                                                                                      
C                                                                                                                                   
      CALL SETLAB(IBUF,NC4TO8("TIME"),1,NC4TO8("PCTC"),LEV,                                                                         
     +                                       NEOF,1,KHEM,1)                                                                         
      CALL PUTFLD2(2,PCTC,IBUF,MAXX)                                                                                                
C                                                                                                                                   
C     * THE 4TH RECORD CONTAINS THE TIME MEAN  (NAME=4HTAVG)                                                                        
C                                                                                                                                   
      DO IJ=1,NLAT*NLON                                                                                                             
        X(IJ)=SPVAL                                                                                                                 
      ENDDO                                                                                                                         
      DO IV=1,NLEN                                                                                                                  
        X(IVAL(IV))=AVG(IVAL(IV))                                                                                                   
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * ADD A MERIDIAN FOR CYCLICITY (IF REQUIRED)                                                                                  
C                                                                                                                                   
      IF (ICYCLE.EQ.1) THEN                                                                                                         
        DO J=1,NLAT                                                                                                                 
          X(J*NLON)=X((J-1)*NLON+1)                                                                                                 
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
      CALL SETLAB(IBUF,KIND,1,NC4TO8("TAVG"),LEV,NLG,NLAT,KHEM,NPCK)                                                                
      CALL PUTFLD2(2,X,IBUF,MAXX)                                                                                                   
C                                                                                                                                   
C     * THE 5TH RECORD CONTAINS THE VARIANCE  (NAME=4HVARI)                                                                         
C                                                                                                                                   
      DO IJ=1,NLAT*NLON                                                                                                             
        X(IJ)=SPVAL                                                                                                                 
      ENDDO                                                                                                                         
      DO IV=1,NLEN                                                                                                                  
        X(IVAL(IV))=VAR(IVAL(IV))                                                                                                   
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * ADD A MERIDIAN FOR CYCLICITY (IF REQUIRED)                                                                                  
C                                                                                                                                   
      IF (ICYCLE.EQ.1) THEN                                                                                                         
        DO J=1,NLAT                                                                                                                 
          X(J*NLON)=X((J-1)*NLON+1)                                                                                                 
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
      CALL SETLAB(IBUF,KIND,1,NC4TO8("VARI"),LEV,NLG,NLAT,KHEM,NPCK)                                                                
      CALL PUTFLD2(2,X,IBUF,MAXX)                                                                                                   
C                                                                                                                                   
C     * THE NEXT 2*NEOF RECORDS ARE PAIRS OF EOFS AND PCS                                                                           
C                                                                                                                                   
      DO NE=1,NEOF                                                                                                                  
        DO IJ=1,NLAT*NLON                                                                                                           
          EVEC(IJ)=SPVAL                                                                                                            
        ENDDO                                                                                                                       
        IF (LTRANS.EQ.0) THEN                                                                                                       
          DO IV=1,NLEN                                                                                                              
            EVEC(IVAL(IV))=COV(IV,INDX(NE))                                                                                         
          ENDDO                                                                                                                     
        ELSE                                                                                                                        
C                                                                                                                                   
C         * IF THE EIGENVECTORS ARE CALCULATED FOR THE TRANSPOSED INPUT DATA                                                        
C         * THEN WE HAVE TO TRANSFORM THEM BACK.                                                                                    
C                                                                                                                                   
          CLOSE (77)                                                                                                                
          OPEN(77,IOSTAT=IRTC,FILE='XXXMKEOF',FORM='UNFORMATTED')                                                                   
          REWIND 77                                                                                                                 
          DO IV=1,NLEN                                                                                                              
C                                                                                                                                   
C           * READ IN THE TRANSPOSED INPUT FIELDS NEEDED FOR TRANSFORMATION                                                         
C                                                                                                                                   
            READ(77) (X(NT),NT=1,NSTEP)                                                                                             
C                                                                                                                                   
C           * FIND VALID VALUES AND SUBTRACT MEAN                                                                                   
C           * SET MISSING VALUES TO 0.                                                                                              
C                                                                                                                                   
            DO NT=1,NSTEP                                                                                                           
              IF(ABS(X(NT)-SPVAL).GT.SPVALT) THEN                                                                                   
                X(NT)=X(NT)-AVG(IVAL(IV))                                                                                           
              ELSE                                                                                                                  
                X(NT)=0.E0                                                                                                          
              ENDIF                                                                                                                 
            ENDDO                                                                                                                   
C                                                                                                                                   
C           * TRANSFORM                                                                                                             
C                                                                                                                                   
            EVEC(IVAL(IV))=0.E0                                                                                                     
            IF (ICOV.EQ.0) THEN                                                                                                     
              DO NT=1,NSTEP                                                                                                         
                EVEC(IVAL(IV))=EVEC(IVAL(IV))+                                                                                      
     1               COV(NT,INDX(NE))*X(NT)                                                                                         
              ENDDO                                                                                                                 
            ELSE                                                                                                                    
              DO NT=1,NSTEP                                                                                                         
                EVEC(IVAL(IV))=EVEC(IVAL(IV))+                                                                                      
     1               COV(NT,INDX(NE))*X(NT)/SQRT(VAR(IVAL(IV)))                                                                     
              ENDDO                                                                                                                 
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * CHOOSE THE SIGN                                                                                                           
C                                                                                                                                   
        NP=0                                                                                                                        
        DO IV=1,NLEN                                                                                                                
          IF (EVEC(IVAL(IV)).GT.0.)NP=NP+1                                                                                          
          IF (EVEC(IVAL(IV)).LT.0.)NP=NP-1                                                                                          
        ENDDO                                                                                                                       
        IF(NP.LT.0)THEN                                                                                                             
          DO IV=1,NLEN                                                                                                              
            EVEC(IVAL(IV))=-EVEC(IVAL(IV))                                                                                          
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * CALCULATE PRINCIPAL COMPONENTS                                                                                            
C                                                                                                                                   
        REWIND 1                                                                                                                    
        DO NT=1,NSTEP                                                                                                               
          CALL GETFLD2(1,X,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                
          IF(.NOT.OK)THEN                                                                                                           
            CALL                                   XIT('MKEOF',-8)                                                                  
          ENDIF                                                                                                                     
C                                                                                                                                   
C         * FIND VALID VALUES AND SUBTRACT MEAN                                                                                     
C         * SET MISSING VALUES TO 0.                                                                                                
C                                                                                                                                   
          DO IV=1,NLEN                                                                                                              
            IF(ABS(X(IVAL(IV))-SPVAL).GT.SPVALT)THEN                                                                                
              X(IVAL(IV))=X(IVAL(IV))-AVG(IVAL(IV))                                                                                 
            ELSE                                                                                                                    
              X(IVAL(IV))=0.E0                                                                                                      
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
          PC(NT)=0.E0                                                                                                               
          DO IV=1,NLEN                                                                                                              
            PC(NT)=PC(NT)+EVEC(IVAL(IV))*X(IVAL(IV))                                                                                
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * NORMALIZE EOFS AND PCS PROPERLY AND CHOOSE THE RIGHT SIGN                                                                 
C       * SO THAT THE MAXIMAL EIGENVECTOR ELEMENT IS POSITIVE.                                                                      
C                                                                                                                                   
        EMAX=0.E0                                                                                                                   
        DO IV=1,NLEN                                                                                                                
          IF (ABS(EVEC(IVAL(IV))).GT.ABS(EMAX)) EMAX=EVEC(IVAL(IV))                                                                 
        ENDDO                                                                                                                       
        SUM=0.E0                                                                                                                    
        DO IV=1,NLEN                                                                                                                
          SUM=SUM+EVEC(IVAL(IV))**2                                                                                                 
        ENDDO                                                                                                                       
        FACT1=SIGN(1.E0,EMAX)/SQRT(SUM)                                                                                             
        FACT2=FACT1                                                                                                                 
        IF (NORM.EQ.1) THEN                                                                                                         
          FACT1=FACT1*SQRT(EVAL(INDX(NE)))                                                                                          
          FACT2=FACT2/SQRT(EVAL(INDX(NE)))                                                                                          
        ENDIF                                                                                                                       
        DO IV=1,NLEN                                                                                                                
          EVEC(IVAL(IV))=EVEC(IVAL(IV))*FACT1                                                                                       
        ENDDO                                                                                                                       
        DO NT=1,NSTEP                                                                                                               
          PC(NT)=PC(NT)*FACT2                                                                                                       
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * ADD LATITUDE FOR CYCLICITY (IF REQUIRED)                                                                                  
C                                                                                                                                   
        IF (ICYCLE.EQ.1) THEN                                                                                                       
          DO J=1,NLAT                                                                                                               
            EVEC(J*NLON)=EVEC((J-1)*NLON+1)                                                                                         
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * SAVE EOFS (NAME=4HEVEC)                                                                                                   
C                                                                                                                                   
        CALL SETLAB(IBUF,KIND,NE,NC4TO8("EVEC"),LEV,NLG,NLAT,KHEM,NPCK)                                                             
        CALL PUTFLD2(2,EVEC,IBUF,MAXX)                                                                                              
C                                                                                                                                   
C       * SAVE PCS (NAME=4H  PC)                                                                                                    
C                                                                                                                                   
        CALL SETLAB(JBUF,NC4TO8("TIME"),NE,NC4TO8("  PC"),LEV,                                                                      
     +                                       NSTEP,1,KHEM,NPCK)                                                                     
        CALL PUTFLD2(2,PC,JBUF,MAXT)                                                                                                
      ENDDO                                                                                                                         
                                                                                                                                    
      CALL                                         XIT('MKEOF',0)                                                                   
 900  CALL                                         XIT('MKEOF',-10)                                                                 
C-------------------------------------------------------------------------                                                          
 5010 FORMAT (10X,4I5,E10.0,F5.1)                                               H4                                                  
 6010 FORMAT (' INPUT PARAMETERS ARE:'/                                                                                             
     1     ' NEOF=',I5,' NORM=',I5,' ICOV=',I5,' ITRAN=',I5,                                                                        
     2     ' SPVAL=',1P1E12.5,' PCTMIN=',0P1F5.1/)                                                                                  
 6020 FORMAT (' ALL (POSSIBLE) EOFS ARE CALCULATED.')                                                                               
 6030 FORMAT (' EOFS ARE NORMALIZED BY UNITY.')                                                                                     
 6035 FORMAT (' PCS  ARE NORMALIZED BY UNITY.')                                                                                     
 6040 FORMAT (' EOFS ARE CALCULATED FOR COVARIANCE MATRIX.')                                                                        
 6045 FORMAT (' EOFS ARE CALCULATED FOR CORRELATION MATRIX.')                                                                       
 6050 FORMAT (' INPUT DATA ARE TRANSPOSED ONLY IF NTIME<LENGTH.')                                                                   
 6055 FORMAT (' INPUT DATA WILL NOT BE TRANSPOSED.')                                                                                
 6060 FORMAT (' INPUT DATA ARE TRANSPOSED')                                                                                         
 6065 FORMAT (' INPUT FILE SEEMS TO HAVE CYCLIC MERIDIAN. EXCLUDE IT.')                                                             
 6070 FORMAT (' NTIME =',I10,' RECORDS PROCESSED.')                                                                                 
 6075 FORMAT (' LENGTH=',I10,' VALID POINTS FOUND.')                                                                                
 6080 FORMAT (' NMISS =',I10,' POINTS EXCLUDED.')                                                                                   
 6090 FORMAT (' FIRST',I5,' EOFS ARE CALCULATED.')                                                                                  
 6100 FORMAT (/' SUM OF ALL EIGENVALUES IS',1P1E12.5,                                                                               
     1        /' FIRST ',I5,' EIGENVALUES ARE:',                                                                                    
     2        /'    NUMBER        EVAL    PCT(%) CUM.PCT(%)')                                                                       
 6105 FORMAT (I10,1P1E12.4,0P2F10.4)                                                                                                
 6110 FORMAT (                                                                                                                      
     1 ' NOTE: EXPLAINED VARIANCES ARE VALID FOR NORMALIZED DATA!',                                                                 
     2 '       BUT PCS ARE CALCULATED FOR THE ORIGINAL DATA.')                                                                      
      END                                                                                                                           
