      PROGRAM VARIMAX                                                                                                               
C     PROGRAM VARIMAX(X,       R,        INPUT,       OUTPUT,           )       H2                                                  
C    1          TAPE1=X, TAPE5=R,  TAPE5=INPUT, TAPE6=OUTPUT)                                                                       
C     -------------------------------------------------------                   H2                                                  
C                                                                               H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       H2                                                  
C     APR 25/03 - SLAVA KHARIN (PROTECT AGAINST POTENTIAL DIVISION BY ZERO)                                                         
C     APR 12/99 - SLAVA KHARIN (FIX A BUG)                                                                                          
C     MAR 09/99 - SLAVA KHARIN.                                                                                                     
C                                                                               H2                                                  
CVARIMAX  - VARIMAX ORTHOGONAL ROTATION.                                1  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - SLAVA KHARIN                                                         H3                                                  
C                                                                               H3                                                  
CPURPOSE - ROTATE ORTHOGONALLY A SET OF GRIDDED FIELDS BY THE VARIMAX METHOD.   H3                                                  
C                                                                               H3                                                  
CINPUT FILE...                                                                  H3                                                  
C                                                                               H3                                                  
C     X = INPUT FILE WITH VECTORS TO BE ROTATED IN CCRN FORMAT.                 H3                                                  
C                                                                               H3                                                  
COUTPUT FILE...                                                                 H3                                                  
C                                                                               H3                                                  
C     R = OUTPUT FILE WITH ROTATED VECTORS.                                     H3                                                  
C                                                                                                                                   
C INPUT PARAMETERS:                                                                                                                 
C                                                                               H5                                                  
C     NVEC  =   NUMBER OF VECTORS TO BE ROTATED.                                H5                                                  
C           = 0 ROTATE ALL VECTORS IN INPUT FILE.                               H5                                                  
C     NORM  = 0 DO NOT PRE-PROCESS INPUT VECTORS.                               H5                                                  
C               OTHERWISE, PERFORM KAISER NORMALIZATION.                        H5                                                  
C               IN KAISER NORMALIZATION, THE I-TH ELEMENT OF EACH VECTOR IS     H5                                                  
C               DIVIDED BY THE SQUARE ROOT OF THE SUM OF THE SQUARED I-TH       H5                                                  
C               ELEMENTS.                                                       H5                                                  
C     SPVAL =   SPECIAL VALUE IN INPUT VECTORS TO INDICATE MISSING VALUES.      H5                                                  
C                                                                               H5                                                  
CEXAMPLES OF INPUT CARD...                                                      H5                                                  
C                                                                               H5                                                  
C*VARIMAX     5    1    1.E+38                                                  H5                                                  
C                                                                               H5                                                  
C (VARIMAX ROTATION OF THE FIRST 5 VECTORS WITH KAISER NORMALIZATION)           H5                                                  
C-------------------------------------------------------------------------------                                                    
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
C                                                                                                                                   
C     NMAX    = MAXIMAL NUMBER OF EOFS.                                                                                             
C     LENMAX  = MAX LENGTH OF VECTORS.                                                                                              
C     ITERMAX = MAXIMAL NUMBER OF ITERATIONS.                                                                                       
C     EPS     = SMALL NUMBER FOR THE CONVERGENCE CRITERION.                                                                         
C                                                                                                                                   
      PARAMETER (LENMAX=65341,NMAX=100,ITERMAX=1000,                                                                                
     1           EPS=1.E-12)                                                                                                        
                                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
C                                                                                                                                   
C     * ARRAY FOR INPUT FIELDS                                                                                                      
C                                                                                                                                   
      REAL X(LENMAX,NMAX), H(LENMAX)                                                                                                
C                                                                                                                                   
      DATA MAXX/65341/                                                                                                              
      DATA PI/3.14159265358979E0/,SMALL/1.E-12/,BIG/1.E+39/                                                                         
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
C     * ASSIGN FILES TO FORTRAN UNITS                                                                                               
C                                                                                                                                   
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,1,2,5,6)                                                                                                      
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
C                                                                                                                                   
C     * READ INPUT PARAMETERS FROM CARD                                                                                             
C                                                                                                                                   
      READ(5,5010,END=900) NVEC,NORM,SPVAL                                      H4                                                  
C                                                                                                                                   
C     * PRINT THE INPUT PARAMETERS TO STANDARD OUTPUT                                                                               
C                                                                                                                                   
      WRITE (6,6010) NVEC,NORM,SPVAL                                                                                                
                                                                                                                                    
      IF (NVEC .EQ. 0)                                                                                                              
     1     WRITE (6,'(A)') ' ROTATE ALL VECTORS IN INPUT FILE.'                                                                     
                                                                                                                                    
C                                                                                                                                   
C     * READ VECTORS                                                                                                                
C                                                                                                                                   
      NR=0                                                                                                                          
 100  CONTINUE                                                                                                                      
      CALL GETFLD2(1,H,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                    
                                                                                                                                    
      IF(.NOT.OK) THEN                                                                                                              
         IF(NR.EQ.0) CALL                          XIT('VARIMAX',-1)                                                                
         IF(NR.LT.NVEC) THEN                                                                                                        
            WRITE(6,'(2A)')                                                                                                         
     1           ' ***WARNING! NUMBER OF VECTORS IN INPUT FILE IS',                                                                 
     2           ' SMALLER THAN NVEC IN INPUT CARD.'                                                                                
            NVEC=NR                                                                                                                 
            GOTO 200                                                                                                                
         ENDIF                                                                                                                      
         IF(NVEC.EQ.0) THEN                                                                                                         
            NVEC=NR                                                                                                                 
            GOTO 200                                                                                                                
         ENDIF                                                                                                                      
      ENDIF                                                                                                                         
                                                                                                                                    
      IF (NR.EQ.0) THEN                                                                                                             
         IF (IBUF(1).EQ.NC4TO8("SPEC") .OR. IBUF(1).EQ.NC4TO8("FOUR"))                                                              
     1        CALL                                 XIT('VARIMAX',-2)                                                                
         NWDS=IBUF(5)*IBUF(6)                                                                                                       
         NPCK=IBUF(8)                                                                                                               
C                                                                                                                                   
C     * GIVE A WARNING IF NPCK > 1                                                                                                  
C                                                                                                                                   
         IF (NPCK.GT.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                                                                                                                      
      ENDIF                                                                                                                         
                                                                                                                                    
      NR=NR+1                                                                                                                       
      DO I=1,NWDS                                                                                                                   
         X(I,NR)=H(I)                                                                                                               
      ENDDO                                                                                                                         
                                                                                                                                    
      IF (NR .EQ. NVEC) GOTO 200                                                                                                    
      GOTO 100                                                                                                                      
                                                                                                                                    
 200  CONTINUE                                                                                                                      
      WRITE(6,'(A,I5,A)')                                                                                                           
     1     ' ROTATE ', NVEC, ' INPUT VECTORS.'                                                                                      
C                                                                                                                                   
C     * NORMALIZE IF NEEDED                                                                                                         
C                                                                                                                                   
      IF (NORM.NE.0) THEN                                                                                                           
         WRITE(6,'(A)')                                                                                                             
     1        ' PERFORM KAISER NORMALIZATION OF UNROTATED VECTORS...'                                                               
C                                                                                                                                   
C     * CALCULATE NORMALIZING CONSTANTS                                                                                             
C                                                                                                                                   
         DO I=1,NWDS                                                                                                                
            H(I) = 0.E0                                                                                                             
            DO N = 1, NVEC                                                                                                          
               IF (ABS(X(I,N)-SPVAL).GT.SMALL)                                                                                      
     1              H(I) = H(I)+X(I,N)*X(I,N)                                                                                       
            ENDDO                                                                                                                   
            IF(H(I).GT.0.) THEN                                                                                                     
              H(I)=SQRT(H(I))                                                                                                       
              DO N = 1, NVEC                                                                                                        
                IF (ABS(X(I,N)-SPVAL).GT.SMALL)                                                                                     
     1              X(I,N) = X(I,N)/H(I)                                                                                            
              ENDDO                                                                                                                 
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
      ENDIF                                                                                                                         
                                                                                                                                    
C                                                                                                                                   
C     * CALCULATE THE VARIMAX CRITERION FOR UNROTATED VECTORS                                                                       
C     * (USE IT AS A NORMALIZATION FACTOR)                                                                                          
C                                                                                                                                   
      ITER=0                                                                                                                        
      VARIM0=0.E0                                                                                                                   
      DO N = 1, NVEC                                                                                                                
         SX2 = 0.E0                                                                                                                 
         SX4 = 0.E0                                                                                                                 
         NVAL = 0                                                                                                                   
         DO I = 1, NWDS                                                                                                             
            IF (ABS(X(I,N)-SPVAL).GT.SMALL) THEN                                                                                    
               SX2 = SX2+X(I,N)**2                                                                                                  
               SX4 = SX4+X(I,N)**4                                                                                                  
               NVAL = NVAL+1                                                                                                        
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         VARIM0 = VARIM0+(FLOAT(NVAL)*SX4-SX2*SX2)                                                                                  
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * ITERATION LOOP                                                                                                              
C                                                                                                                                   
      ITER=1                                                                                                                        
      VARIM1=BIG                                                                                                                    
 300  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * ROTATE PAIRS OF VECTORS                                                                                                     
C                                                                                                                                   
      DO N2=1, NVEC                                                                                                                 
         DO N1=N2+1, NVEC                                                                                                           
C                                                                                                                                   
C     * CALCULATE THE ROTATION ANGLE                                                                                                
C                                                                                                                                   
            A = 0.E0                                                                                                                
            B = 0.E0                                                                                                                
            C = 0.E0                                                                                                                
            D = 0.E0                                                                                                                
            DO I = 1, NWDS                                                                                                          
               IF ( ABS(X(I,N1)-SPVAL).GT.SMALL .AND.                                                                               
     1              ABS(X(I,N2)-SPVAL).GT.SMALL) THEN                                                                               
                  U = X(I,N1)*X(I,N1)-X(I,N2)*X(I,N2)                                                                               
                  V = 2.E0*X(I,N1)*X(I,N2)                                                                                          
                  A = A+U                                                                                                           
                  B = B+V                                                                                                           
                  C = C+U*U-V*V                                                                                                     
                  D = D+2.E0*U*V                                                                                                    
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
            PHI=ATAN2(D*FLOAT(NVAL)-2.E0*A*B,C*FLOAT(NVAL)-A*A+B*B)/4.E0                                                            
C                                                                                                                                   
C     * ROTATE THE PAIR                                                                                                             
C                                                                                                                                   
            CS = COS(PHI)                                                                                                           
            SN = SIN(PHI)                                                                                                           
            DO I = 1, NWDS                                                                                                          
               IF ( ABS(X(I,N1)-SPVAL).GT.SMALL .AND.                                                                               
     1              ABS(X(I,N2)-SPVAL).GT.SMALL ) THEN                                                                              
                  WORK   = X(I,N1)*CS+X(I,N2)*SN                                                                                    
                  X(I,N2)=-X(I,N1)*SN+X(I,N2)*CS                                                                                    
                  X(I,N1)= WORK                                                                                                     
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
C                                                                                                                                   
C     * END OF LOOP                                                                                                                 
C                                                                                                                                   
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE THE VARIMAX CRITERION FOR ROTATED VECTORS                                                                         
C                                                                                                                                   
      VARIM=0.E0                                                                                                                    
      DO N = 1, NVEC                                                                                                                
         SX2 = 0.E0                                                                                                                 
         SX4 = 0.E0                                                                                                                 
         NVAL = 0                                                                                                                   
         DO I = 1, NWDS                                                                                                             
            IF (ABS(X(I,N)-SPVAL).GT.SMALL) THEN                                                                                    
               SX2 = SX2+X(I,N)**2                                                                                                  
               SX4 = SX4+X(I,N)**4                                                                                                  
               NVAL = NVAL+1                                                                                                        
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         VARIM = VARIM+(FLOAT(NVAL)*SX4-SX2**2)                                                                                     
      ENDDO                                                                                                                         
      VARIM = VARIM/VARIM0                                                                                                          
C      WRITE (6,*) ' ITER=', ITER, ' VARIM = ', VARIM,                                                                              
C     1     ' DELTA=',ABS(VARIM-VARIM1)                                                                                             
C                                                                                                                                   
C     * CHECK CONVERGENCY                                                                                                           
C                                                                                                                                   
      IF (ABS(VARIM-VARIM1).LE.EPS .OR. ITER.GT.ITERMAX) GOTO 400                                                                   
                                                                                                                                    
      ITER = ITER+1                                                                                                                 
      VARIM1=VARIM                                                                                                                  
      GOTO 300                                                                                                                      
                                                                                                                                    
 400  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * REMOVE KAISER NORMALIZATION                                                                                                 
C                                                                                                                                   
      IF (NORM.NE.0) THEN                                                                                                           
         WRITE(6,'(A)')                                                                                                             
     1        ' REMOVE KAISER NORMALIZATION FROM THE ROTATED VECTORS.'                                                              
         DO I=1,NWDS                                                                                                                
            DO N = 1, NVEC                                                                                                          
               IF (ABS(X(I,N)-SPVAL).GT.SMALL)                                                                                      
     1              X(I,N)=X(I,N)*H(I)                                                                                              
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * SAVE                                                                                                                        
C                                                                                                                                   
      DO N=1,NVEC                                                                                                                   
         IBUF(2)=N                                                                                                                  
         CALL PUTFLD2(2,X(1,N),IBUF,MAXX)                                                                                           
      ENDDO                                                                                                                         
                                                                                                                                    
      CALL                                         XIT('VARIMAX',0)                                                                 
 900  CALL                                         XIT('VARIMAX',-3)                                                                
C-------------------------------------------------------------------------                                                          
 5010 FORMAT (10X,2I5,1E10.0)                                                   H4                                                  
 6010 FORMAT (' NVEC=',I5, ' NORM=',I5, ' SPVAL=',1P1E12.5)                                                                         
      END                                                                                                                           
