      PROGRAM PERMUT                                                                                                                
C     PROGRAM PERMUT (X,           Y,         XBAR,         SX,                 H2                                                  
C    1                YBAR,        SY,        INPUT,        OUTPUT,     )       H2                                                  
C    2         TAPE1 =X,    TAPE2 =Y,  TAPE11=XBAR,  TAPE12=SX,                                                                     
C    3         TAPE13=YBAR, TAPE14=SY, TAPE5 =INPUT, TAPE6 =OUTPUT)                                                                 
C     -------------------------------------------------------------             H2                                                  
C                                                                               H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       H2                                                  
C     JUL 22/92 - E. CHAN  (REPLACE UNFORMATTED I/O WITH I/O ROUTINES)                                                              
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     MAR 28/85 - F. ZWIERS                                                                                                         
C                                                                               H2                                                  
CPERMUT  - FIRST STEP OF PERMUTATION TESTING PROCEDURE.                 2  4 C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - F. ZWIERS                                                            H3                                                  
C                                                                               H3                                                  
CPURPOSE - STEP ONE OF A PERMUTATION PROCEDURE FOR ASCRIBING A SIGNFICANCE      H3                                                  
C          LEVEL TO THE DIFFERENCES IN MEANS AND VARIABLITY FOR TWO FIELDS.     H3                                                  
C                                                                               H3                                                  
CINPUT FILES...                                                                 H3                                                  
C                                                                               H3                                                  
C                                                                               H3                                                  
C      X    = FILE OF INDEPENDENT REALIZATIONS OF MEAN FIELDS FOR X.            H3                                                  
C      Y    = FILE OF INDEPENDENT REALIZATIONS OF MEAN FIELDS FOR Y.            H3                                                  
C                                                                               H3                                                  
COUTPUT FILES...                                                                H3                                                  
C                                                                               H3                                                  
C      XBAR = FILE OF MEANS. THE FIRST RECORD IS THE MEAN OF XBAR. THE          H3                                                  
C             SUBSEQUENT RECORDS ARE MEANS OF NX RECORDS  SAMPLED FROM          H3                                                  
C             THE POOL OF RECORDS IN FILES X AND Y. THE NX IS THE NUM-          H3                                                  
C             BER OF RECORDS IN FILE X.                                         H3                                                  
C      SX   = FILE  OF  STANDARD DEVIATIONS.  THE FIRST  RECORD IS THE          H3                                                  
C             STANDARD DEVIATION OF X. THE SUBSEQUENT NPERM-1  RECORDS          H3                                                  
C             ARE STANDARD DEVIATIONS OF THE NPERM-1 SAMPLES  TAKEN TO          H3                                                  
C             COMPUTE THE NPERM-1 MEANS IN XBAR.                                H3                                                  
C      YBAR = FILE OF MEANS AS ABOVE  EXCEPT THAT FIRST RECORD  IS THE          H3                                                  
C             THE MEAN OF Y AND THE SAMPLES ARE OF SIZE NY,(THE NUMBER          H3                                                  
C             OF RECORDS IN Y).                                                 H3                                                  
C      SY   = THE FILE OF CORRESPONDING STANDARD DEVIATIONS                     H3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               H5                                                  
C      NPERM = NUMBER OF TIMES TO RESAMPLE THE POOL OF RECORDS IN FILES         H5                                                  
C              X AND Y.                                                         H5                                                  
C      DSEED = SEED NUMBER FOR THE (IMSL) RANDOM PERMUTATION GENERATOR.         H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARD...                                                       H5                                                  
C                                                                               H5                                                  
C*PERMUT.  1000    12345.                                                       H5                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      DIMENSION X(18528),POOL(741120),XBAR(18528),YBAR(18528),                                                                      
     1          SX(18528), SY(18528),IBUFX(8),IBUFY(8),J(40)                                                                        
      LOGICAL    OK                                                                                                                 
      DOUBLE PRECISION DSEED                                                                                                        
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      DATA     MAXX/18528/,MPOOL/40/                                                                                                
C     DATA     NSPEC/4HSPEC/                                                                                                        
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
C                                                                                                                                   
      NSPEC=NC4TO8("SPEC")                                                                                                          
      NFILES=8                                                                                                                      
      CALL JCLPNT(NFILES,1,2,11,12,13,14,5,6)                                                                                       
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
C                                                                                                                                   
C     * GET THE PERMUT CARD                                                                                                         
C                                                                                                                                   
      READ(5,5010,END=10) NPERM,DSEED                                           H4                                                  
      GOTO 20                                                                                                                       
   10 CONTINUE                                                                                                                      
         WRITE(6,6010)                                                                                                              
         CALL                                      XIT('PERMUT',-1)                                                                 
   20 CONTINUE                                                                                                                      
      WRITE(6,6015) NPERM,DSEED                                                                                                     
C                                                                                                                                   
C     * CHECK THAT THE INPUT FILES CONTAIN DATA AND COMPARE THE LABELS                                                              
C                                                                                                                                   
      CALL FBUFFIN(1,IBUFX,-8,K,LEN)                                                                                                
      IF (K.GE.0) GOTO 30                                                                                                           
      REWIND 1                                                                                                                      
      GOTO 40                                                                                                                       
   30 CONTINUE                                                                                                                      
         WRITE(6,6020)                                                                                                              
         CALL                                      XIT('PERMUT',-2)                                                                 
   40 CONTINUE                                                                                                                      
      CALL FBUFFIN(2,IBUFY,-8,K,LEN)                                                                                                
      IF (K.GE.0) GOTO 50                                                                                                           
      REWIND 2                                                                                                                      
      GOTO 60                                                                                                                       
   50 CONTINUE                                                                                                                      
         WRITE(6,6030)                                                                                                              
         CALL                                      XIT('PERMUT',-3)                                                                 
   60 CONTINUE                                                                                                                      
      IF( IBUFX(1).NE.IBUFY(1) .OR. IBUFX(3).NE.IBUFY(3) .OR.                                                                       
     1    IBUFX(4).NE.IBUFY(4) .OR. IBUFX(5).NE.IBUFY(5) .OR.                                                                       
     2                              IBUFX(6).NE.IBUFY(6) )THEN                                                                      
         WRITE(6,6040)                                                                                                              
         CALL                                      XIT('PERMUT',-4)                                                                 
      ENDIF                                                                                                                         
      NWDS=IBUFX(5)*IBUFX(6)                                                                                                        
      IF(IBUFX(1).EQ.NSPEC) NWDS=NWDS*2                                                                                             
C                                                                                                                                   
C     * LOAD THE DATA                                                                                                               
C                                                                                                                                   
      NX=0                                                                                                                          
  100 CONTINUE                                                                                                                      
         CALL GETFLD2(1,X,-1,0,0,0,IBUF,MAXX,OK)                                                                                    
         IF(OK)THEN                                                                                                                 
            NX=NX+1                                                                                                                 
            IF(NX.GT.MPOOL)THEN                                                                                                     
               WRITE(6,6050)                                                                                                        
               CALL                                XIT('PERMUT',-5)                                                                 
            ENDIF                                                                                                                   
            IBASE=(NX-1)*NWDS                                                                                                       
            DO 110 I=1,NWDS                                                                                                         
               POOL(IBASE+I)=X(I)                                                                                                   
  110       CONTINUE                                                                                                                
            GOTO 100                                                                                                                
         ENDIF                                                                                                                      
C                                                                                                                                   
      WRITE(6,6055) NX                                                                                                              
C                                                                                                                                   
      NY=0                                                                                                                          
  150 CONTINUE                                                                                                                      
         CALL GETFLD2(2,X,-1,0,0,0,IBUF,MAXX,OK)                                                                                    
         IF(OK)THEN                                                                                                                 
            NY=NY+1                                                                                                                 
            IF(NX+NY.GT.MPOOL)THEN                                                                                                  
               WRITE(6,6060)                                                                                                        
               CALL                                XIT('PERMUT',-6)                                                                 
            ENDIF                                                                                                                   
            IBASE=(NX+NY-1)*NWDS                                                                                                    
            DO 160 I=1,NWDS                                                                                                         
               POOL(IBASE+I)=X(I)                                                                                                   
  160       CONTINUE                                                                                                                
            GOTO 150                                                                                                                
         ENDIF                                                                                                                      
C                                                                                                                                   
      WRITE(6,6065) NY                                                                                                              
C                                                                                                                                   
C                                                                                                                                   
      IBUFX(2)=NX                                                                                                                   
      IBUFY(2)=NY                                                                                                                   
C                                                                                                                                   
C     * COMPUTE STATISTICS FOR SAMPLES OF SIZE NX AND NY TAKEN FROM THE                                                             
C     * POOL OF RECORDS STORED IN ARRAY POOL.  THE FIRST SET OF                                                                     
C     * STATISTICS ARE FOR FILES X AND Y.                                                                                           
C                                                                                                                                   
      NPOOL=NX+NY                                                                                                                   
      NXP1=NX+1                                                                                                                     
      FNX=FLOAT(NX)                                                                                                                 
      FNY=FLOAT(NY)                                                                                                                 
      FNXM1=FNX-1.0E0                                                                                                               
      FNYM1=FNY-1.0E0                                                                                                               
      DO 200 I=1,NPOOL                                                                                                              
         J(I)=I                                                                                                                     
 200  CONTINUE                                                                                                                      
C                                                                                                                                   
C                                                                                                                                   
      DO 2000 L=1,NPERM                                                                                                             
C                                                                                                                                   
C        * ZERO THE STATS ARRAYS                                                                                                    
C                                                                                                                                   
         DO 1010 I=1,NWDS                                                                                                           
            XBAR(I)=0.0E0                                                                                                           
            YBAR(I)=0.0E0                                                                                                           
            SX(I)=0.0E0                                                                                                             
            SY(I)=0.0E0                                                                                                             
 1010    CONTINUE                                                                                                                   
C                                                                                                                                   
C        * ACCUMULATE SUMS AND SUMS OF SQUARES                                                                                      
C                                                                                                                                   
         DO 1050 I=1,NX                                                                                                             
            KBASE=(J(I)-1)*NWDS                                                                                                     
            DO 1040 K=1,NWDS                                                                                                        
               XBAR(K)=XBAR(K)+POOL(KBASE+K)                                                                                        
               SX(K)=SX(K)+POOL(KBASE+K)*POOL(KBASE+K)                                                                              
 1040       CONTINUE                                                                                                                
 1050    CONTINUE                                                                                                                   
         DO 1100 I=NXP1,NPOOL                                                                                                       
            KBASE=(J(I)-1)*NWDS                                                                                                     
            DO 1090 K=1,NWDS                                                                                                        
               YBAR(K)=YBAR(K)+POOL(KBASE+K)                                                                                        
               SY(K)=SY(K)+POOL(KBASE+K)*POOL(KBASE+K)                                                                              
 1090       CONTINUE                                                                                                                
 1100    CONTINUE                                                                                                                   
C                                                                                                                                   
C        * COMPUTE MEANS AND STANDARD DEVIATIONS                                                                                    
C                                                                                                                                   
         DO 1150 I=1,NWDS                                                                                                           
            XBAR(I)=XBAR(I)/FNX                                                                                                     
            YBAR(I)=YBAR(I)/FNY                                                                                                     
 1150    CONTINUE                                                                                                                   
         DO 1160 I=1,NWDS                                                                                                           
            SX(I)=SQRT((SX(I)-FNX*XBAR(I)*XBAR(I))/FNXM1)                                                                           
            SY(I)=SQRT((SY(I)-FNY*YBAR(I)*YBAR(I))/FNYM1)                                                                           
 1160    CONTINUE                                                                                                                   
C                                                                                                                                   
C        * OUTPUT THE RESULTS                                                                                                       
C                                                                                                                                   
         DO 1200 I=1,8                                                                                                              
            IBUF(I)=IBUFX(I)                                                                                                        
 1200    CONTINUE                                                                                                                   
         CALL PUTFLD2(11,XBAR,IBUF,MAXX)                                                                                            
         CALL PUTFLD2(12,SX,IBUF,MAXX)                                                                                              
         DO 1210 I=1,8                                                                                                              
            IBUF(I)=IBUFY(I)                                                                                                        
 1210    CONTINUE                                                                                                                   
         CALL PUTFLD2(13,YBAR,IBUF,MAXX)                                                                                            
         CALL PUTFLD2(14,SY,IBUF,MAXX)                                                                                              
C                                                                                                                                   
C        * PERMUTE THE RECORD LABELS                                                                                                
C                                                                                                                                   
         CALL GGPER(DSEED,NPOOL,J)                                                                                                  
C                                                                                                                                   
C        * DO IT ALL AGAIN                                                                                                          
C                                                                                                                                   
 2000 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * ALL DONE                                                                                                                    
C                                                                                                                                   
      CALL                                         XIT('PERMUT',0)                                                                  
C                                                                                                                                   
C-----------------------------------------------------------------------------                                                      
C                                                                                                                                   
 5010 FORMAT(10X,I5,E10.0)                                                      H4                                                  
 6010 FORMAT('0THE PERMUT CARD IS MISSING.')                                                                                        
 6015 FORMAT('0PERMUT:  NPERM,DSEED=',I5,E15.5)                                                                                     
 6020 FORMAT('0FILE X IS EMPTY.')                                                                                                   
 6030 FORMAT('0FILE Y IS EMPTY.')                                                                                                   
 6040 FORMAT('0LABELS ON THE FIRST RECORDS OF FILES X AND Y DO',                                                                    
     1       ' NOT MATCH.')                                                                                                         
 6050 FORMAT('0FILE X CONTAINS TOO MANY RECORDS.')                                                                                  
 6055 FORMAT('0FOUND ',I2,' RECORDS IN FILE X.')                                                                                    
 6060 FORMAT('0THE POOL OF RECORDS IN FILES X AND Y IS TOO LARGE.')                                                                 
 6065 FORMAT('0FOUND ',I2,' RECORDS IN FILE Y.')                                                                                    
      END                                                                                                                           
