      PROGRAM STD                                                                                                                   
C     PROGRAM STD (TSIN,       TSMEAN,       STDEV,       W2STAT,               H2                                                  
C    1                                       INPUT,       OUTPUT,       )       H2                                                  
C    2       TAPE1=TSIN, TAPE2=TSMEAN, TAPE3=STDEV, TAPE4=W2STAT,                                                                   
C    3                                 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     DEC 17/85 - B.DUGAS.                                                                                                          
C                                                                               H2                                                  
CSTD     - COMPUTES ST. DEV. FOR TIME SERIES WITH ANNUAL CYCLES         2  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - B.DUGAS                                                              H3                                                  
C                                                                               H3                                                  
CPURPOSE - CE PROGRAMME CALCULE L'ECART  TYPE DES CY CYCLES  ANNUELS DES        H3                                                  
C          SERIES TEMPORELLES  TSIN A CHAQUE  BOITE DE CE CYCLE  ANNUEL.        H3                                                  
C          TSIN ET TSMEAN SONT PRODUITS PAR LE PROGRAMME BINNING. AINSI,        H3                                                  
C          UNE BOITE DE TSIN N'EST SENSEE CONTENIR QU'UN SEUL POINT. UNE        H3                                                  
C          SERIE DE TSMEAN CONTIENT CY FOIS LE CYCLE ANNUEL DE LA  SERIE        H3                                                  
C          CORRESPONDANTE DE TSIN.                                              H3                                                  
C          ON CALCULE OPTIONNNELLEMENT, W2STAT QUI EST LA STATISTIQUE DE        H3                                                  
C          CRAMER-VON MISES  EVALUANT  LA NORMALITE  DE LA  DISTRIBUTION        H3                                                  
C          DES  CY CYCLES  ANNUELS. UNE VALEUR  EST PRODUITE POUR CHAQUE        H3                                                  
C          BOITE.  LE NOMBRE TOTAL DE CES VALEURS EST DONC  (LONGUEUR DE        H3                                                  
C          TSIN) / CY.                                                          H3                                                  
C          (VOIR " EDF STATISTICS FOR GOODNESS  OF FIT AND SOME COMPARI-        H3                                                  
C          SONS ", M.A. STEPHENS, JOURNAL  OF THE AMERICAN  STATISTICAL         H3                                                  
C          ASSOCIATION, SEPTEMBER 1974, PP 730-737.)                            H3                                                  
C                                                                               H3                                                  
CINPUT FILES...                                                                 H3                                                  
C                                                                               H3                                                  
C      TSIN   = TIME SERIES WITH 2 OR MORE CYCLES                               H3                                                  
C      TSMEAN = CONTAINS THE TIME SERIES OF AVERAGES                            H3                                                  
C                                                                               H3                                                  
COUTPUT FILES...                                                                H3                                                  
C                                                                               H3                                                  
C      STDEV  = CONTAINS STANDARD DEVIATION                                     H3                                                  
C      W2STAT = CRAMER-VON MISES TEST RESULT (OPTIONAL)                         H3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               H5                                                  
C      CY  = NUMBER OF CYCLES PER "TSIN" TIME SERIES                            H5                                                  
C            (GREATER THAN 1 AND LESS THAN "TSIN"TIME SERIES LENGTH)            H5                                                  
C      OPT = SWITCH TO APPLY THE CRAMER-VON MISES TEST OR NOT                   H5                                                  
C            DEFAULT IS ZERO THAT IS NOT TO APPLY THE TEST                      H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARD...                                                       H5                                                  
C                                                                               H5                                                  
C* STD        5    0                                                            H5                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK, OPTION                                                                                                            
                                                                                                                                    
      INTEGER  BC, CY, OPT, TYPE                                                                                                    
      REAL TSMEAN(18528), TSIN  (18528), STDEV(18528),                                                                              
     1     W2STAT(18528), VALEUR( 9264)                                                                                             
                                                                                                                                    
      COMMON /ICOM/ IBUF(8), IDAT(18528)                                                                                            
                                                                                                                                    
      DATA  MAXX   / 18528 /                                                                                                        
      DATA  ZERO  /   0.0E0   /, UN / 1.0E0 /, DEUX / 2.0E0 /                                                                       
      DATA OPTION / .FALSE. /                                                                                                       
                                                                                                                                    
C---------------------------------------------------------------------                                                              
      NF = 6                                                                                                                        
      CALL JCLPNT (NF,1,2,3,4,5,6)                                                                                                  
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
                                                                                                                                    
C     * LIRE LE PREMIER TAMPON D'INFORMATION DE TSIN.                                                                               
                                                                                                                                    
      CALL FBUFFIN(1,IBUF,-8,K,LEN)                                                                                                 
      IF (K.GE.0) GOTO 901                                                                                                          
      IF (IBUF(6).NE.1) CALL                       XIT('STD',-1)                                                                    
      REWIND 1                                                                                                                      
                                                                                                                                    
C     * SAUVER CERTAINES VALEURS A DES FINS DE COMPARAISONS.                                                                        
                                                                                                                                    
      TYPE    = IBUF(1)                                                                                                             
      NOM     = IBUF(3)                                                                                                             
      NIVEAU  = IBUF(4)                                                                                                             
      LONG    = IBUF(5)                                                                                                             
                                                                                                                                    
      IF (LONG.GT.MAXX) CALL                       XIT('STD',-2)                                                                    
                                                                                                                                    
C     * LIRE LA CARTE DE CONTROLE ET S'ASSURER QUE CY EST PLUS                                                                      
C     * GRAND QUE UN ET PLUS PETIT OU EGAL A LA LONGUEUR DE TSIN.                                                                   
                                                                                                                                    
      READ(5,5000,END=902)  CY, OPT                                             H4                                                  
      IF (CY.LE.1 .OR. CY.GT.LONG )  CALL          XIT('STD',-3)                                                                    
                                                                                                                                    
      WRITE(6,6000) CY                                                                                                              
      WRITE(6,6005) IBUF                                                                                                            
                                                                                                                                    
      IF (OPT.NE.0)                                            THEN                                                                 
          REWIND 4                                                                                                                  
          OPTION = .TRUE.                                                                                                           
          WRITE(6,6040)                                                                                                             
      END IF                                                                                                                        
                                                                                                                                    
C     * BC EST LE NOMBRE DE BOITES PAR CYCLES. VERIFIER QU'UN                                                                       
C     * CYCLE CORRESPOND BIEN A NOMBRE ENTIER DE POINTS.                                                                            
                                                                                                                                    
      BC    = LONG/CY                                                                                                               
      OCYM1 = UN/ FLOAT(CY-1)                                                                                                       
      O12CY = UN/(FLOAT(CY)*12.0E0)                                                                                                 
      O2CY  = UN/(FLOAT(CY)*DEUX)                                                                                                   
                                                                                                                                    
                                                                                                                                    
      IF (LONG-BC*CY.NE.0)                                     THEN                                                                 
          WRITE(6,6030)                                                                                                             
          CALL                                     XIT('STD',-4)                                                                    
      END IF                                                                                                                        
                                                                                                                                    
C---------------------------------------------------------------------                                                              
C     * ALLER CHERCHER LA PROCHAINE SERIE.                                                                                          
                                                                                                                                    
      NSERIE = 0                                                                                                                    
  200 CALL GETFLD2(1,TSIN,TYPE,-1,NOM,NIVEAU,IBUF,MAXX,OK)                                                                          
                                                                                                                                    
          IF (.NOT.OK)                                         THEN                                                                 
              WRITE(6,6010) NSERIE                                                                                                  
              IF (NSERIE.EQ.0) CALL                XIT('STD',-5)                                                                    
              CALL                                 XIT('STD',0)                                                                     
          END IF                                                                                                                    
                                                                                                                                    
          IF (IBUF(5).NE.LONG .OR. IBUF(6).NE.1)               THEN                                                                 
              WRITE(6,6020) NSERIE+1, IBUF                                                                                          
              CALL                                 XIT('STD',-6)                                                                    
          END IF                                                                                                                    
                                                                                                                                    
          CALL GETFLD2(2,TSMEAN,TYPE,-1,NOM,NIVEAU,IBUF,MAXX,OK)                                                                    
                                                                                                                                    
          IF (.NOT.OK) CALL                        XIT('STD',-7)                                                                    
                                                                                                                                    
          IF (IBUF(5).NE.LONG .OR. IBUF(6).NE.1)               THEN                                                                 
              WRITE(6,6020) NSERIE+1, IBUF                                                                                          
              CALL                                 XIT('STD',-8)                                                                    
          ENDIF                                                                                                                     
                                                                                                                                    
C         * INITIALISER LE PREMIER CYCLE DE STDEV A ZERO.                                                                           
                                                                                                                                    
          DO 300 L=1,BC                                                                                                             
              STDEV(L) = ZERO                                                                                                       
  300     CONTINUE                                                                                                                  
                                                                                                                                    
C         * TROUVER L'ECART TYPE DE CHAQUE BOITE DE TSIN.                                                                           
C         * LA METTRE DANS LE PREMIER CYCLE DE STDEV.                                                                               
                                                                                                                                    
          DO 400 I=1,CY                                                                                                             
              IJ = (I-1)*BC                                                                                                         
              DO 400 J=1,BC                                                                                                         
                  STDEV(J) = STDEV(J)+(TSIN(J+IJ)-TSMEAN(J+IJ))**2                                                                  
  400     CONTINUE                                                                                                                  
                                                                                                                                    
          DO 500 J=1,BC                                                                                                             
              STDEV(J) = SQRT(STDEV(J)*OCYM1)                                                                                       
  500     CONTINUE                                                                                                                  
                                                                                                                                    
C         * REPETER CE PREMIER CYCLE CY FOIS.                                                                                       
                                                                                                                                    
          DO 600 I=2,CY                                                                                                             
              IJ = (I-1)*BC                                                                                                         
                                                                                                                                    
              DO 600 J=1,BC                                                                                                         
                  STDEV(J+IJ) = STDEV(J)                                                                                            
  600     CONTINUE                                                                                                                  
                                                                                                                                    
C         * (OPTIONNELLEMENT,) CALCULER W2STAT POUR CHAQUE BOITE.                                                                   
                                                                                                                                    
          IF (OPTION)                                          THEN                                                                 
                                                                                                                                    
              DO 800 J=1,BC                                                                                                         
                                                                                                                                    
                  U  = TSMEAN(J)                                                                                                    
                  S  = STDEV(J)                                                                                                     
                  OS = UN/S                                                                                                         
                  W2 = ZERO                                                                                                         
                                                                                                                                    
C                 * EXTRAIRE LES CY VALEURS CORRESPONDANTS A                                                                        
C                 * CETTE BOITE.                                                                                                    
                                                                                                                                    
                  INDEX = J                                                                                                         
                  DO 650 I=1,CY                                                                                                     
                      VALEUR(I) = TSIN(INDEX)                                                                                       
                      INDEX     = INDEX+BC                                                                                          
  650             CONTINUE                                                                                                          
                                                                                                                                    
C                 * LES METTRE EN ORDRE CROISSANT ET LES NORMALISER.                                                                
                                                                                                                                    
                  DO 700 I=1,CY                                                                                                     
                      VALEUR(I) = (VALEUR(I)-U)*OS                                                                                  
  700             CONTINUE                                                                                                          
                                                                                                                                    
                  CALL VSRTA (VALEUR,CY)                                                                                            
                                                                                                                                    
C                 * CALCULER LES PROBABILITES NORMALES ASSOCIEES A                                                                  
C                 * VALEUR PAR  APPROXIMATIONS  RATIONNELLES.  LES                                                                  
C                 * METTRE DANS VALEUR.                                                                                             
                                                                                                                                    
                  CALL NORMAR (VALEUR,CY)                                                                                           
                                                                                                                                    
C                 * EVALUER W2 ET APPLIQUER LA CORRECTION DUE A LA                                                                  
C                 * PETITESSE DE L'ECHANTILLON.                                                                                     
                                                                                                                                    
                  DO 750 I=1,CY                                                                                                     
                      W2 = W2+(VALEUR(I)-(DEUX*FLOAT(I)-1)*O2CY)**2                                                                 
  750             CONTINUE                                                                                                          
                  W2STAT(J) = (W2+O12CY)*(1.0E0+O2CY)                                                                               
                                                                                                                                    
  800         CONTINUE                                                                                                              
                                                                                                                                    
          END IF                                                                                                                    
                                                                                                                                    
C         * SAUVER STDEV ET (OPTIONNELLEMENT) W2STAT.                                                                               
                                                                                                                                    
          CALL PUTFLD2 (3,STDEV,IBUF,MAXX)                                                                                          
                                                                                                                                    
          IF (OPTION)                                          THEN                                                                 
              IBUF(5) = BC                                                                                                          
              CALL PUTFLD2 (4,W2STAT,IBUF,MAXX)                                                                                     
          END IF                                                                                                                    
                                                                                                                                    
          NSERIE = NSERIE+1                                                                                                         
                                                                                                                                    
      GOTO 200                                                                                                                      
                                                                                                                                    
C     * FIN PREMATUREE DU FICHIER TSIN.                                                                                             
                                                                                                                                    
  901 CALL                                         XIT('STD',-9)                                                                    
                                                                                                                                    
C     * E.O.F. SUR INPUT (IN).                                                                                                      
                                                                                                                                    
  902 CALL                                         XIT('STD',-10)                                                                   
                                                                                                                                    
C---------------------------------------------------------------------                                                              
 5000 FORMAT(10X,2I5)                                                           H4                                                  
 6000 FORMAT('0INPUT TIME SERIES CONSIST OF ',I5,' CYCLES.')                                                                        
 6005 FORMAT('0STD ON : ',A4,1X,I10,1X,A4,I10,2I5,I10,I5)                                                                           
 6010 FORMAT('0STD READ ',I5,' SERIES.')                                                                                            
 6020 FORMAT('0RECORD NO,=',I5,', LABEL = ',A4,1X,I10,1X,                                                                           
     1                                      A4,I10,2I5,I10,I5)                                                                      
 6030 FORMAT('0THE NUMBER OF BOXES PER CYCLE IS NOT AN INTEGER.')                                                                   
 6040 FORMAT('0CRAMER-VON MISES STATISTICS TO BE CALCULATED.')                                                                      
      END                                                                                                                           
