C*********************************************************************
C     1950年以降について、北極振動指数と各地の気温の相関を調べる。
C*********************************************************************
      PARAMETER(N=58,M=12, MP=13) 
C                          ! N:1950-2007年までの58年間,M:1-12月
C                          ! データ個数はN*M
      CHARACTER FN*3
      DIMENSION AO(N,M),AT(N,MP) ! AO:北極振動指数, AT:気温
      DIMENSION X(N),Y(N),IY(N) ! 特定の期間のデータのみピックアップ
C                                                              して格納
C ---------------------------------------------------------------------
      WRITE(6,*) 'INPUT FILE NAME : ' ! 画面に '' の内容を書き出す
      READ(5,500) FN    ! FN は3文字の文字 --> CHARACTER FN*3
  500 FORMAT(A3)
C
      OPEN(11,FILE='AO.TXT') ! AO データファイルをオープン
      OPEN(12,FILE=FN//'.TXT') ! AT データファイルをオープン
C
C                        !  // は文字を結合させる操作
C
                             
      OPEN(13,FILE='AO'//FN//'.TXT') ! 特定期間のAO,ATデータを出力する
C                                    ! ここでは、11-12月のデータ
C
C =====================================================================
C     北極振動指数と気温データの読み込み（月平均データ）
C ---------------------------------------------------------------------
      DO 100 I=1,N
        DO 10 J=1,M
          READ(11,*) IY(I), IM, AO(I,J)
   10   CONTINUE
  100 CONTINUE
C
      DO 110 I=1,N
          READ(12,*) IDUM,AT(I,1),AT(I,2),AT(I,3),AT(I,4),AT(I,5)
     &              ,AT(I,6),AT(I,7),AT(I,8),AT(I,9),AT(I,10)
     &              ,AT(I,11),AT(I,12),AT(I,13)
  110 CONTINUE

C =====================================================================
C     11月と12月の北極振動指数及び気温をピックアップし、
C     11-12月の平均値（2ヶ月平均値）を求める
C ---------------------------------------------------------------------
C
C*****************************ここに注意！*****************************
C >>>>>>>> 12/2に配布したプログラムから修正すべきところです <<<<<<<<
      DO 12 I=1,N     !配列X(I),Y(I)にゼロを与えておく
        X(I)=0.0      !　これはDO 20内のX(I)=AOTMP+X(I)の計算で 
        Y(I)=0.0      !　Jが11のとき（つまり）最初にIFでYESと
   12 CONTINUE        !判定されたとき、X(I)は定まらないから
C
C**********************************************************************
C
      DO 200 I=1,N
        DO 20 J=1,M
          IF(J.GE.11) THEN                               ! ※に補足
            AOTMP=AO(I,J)
            ATTMP=AT(I,J)
            X(I)=AOTMP+X(I)
            Y(I)=ATTMP+Y(I)
          END IF
   20 CONTINUE
        X(I)=0.5*X(I) 
        Y(I)=0.5*Y(I) 
        WRITE(13,*) IY(I),X(I),Y(I)
  200 CONTINUE
C
C　           (※) J=1,Mとしているが、Mは12でありAT(I,13）は使わない。
C　            実は、AT(I,13)は年平均気温の値が入っているから。
C
C =====================================================================
C     相関係数の計算副プログラムを呼び出す
C ---------------------------------------------------------------------
      CALL PEARSN(X,Y,N,R)
C
C =====================================================================
C     相関係数を画面に出力
C ---------------------------------------------------------------------
      WRITE(6,*) 'R = ', R
C
C =====================================================================
C     装置（ファイル）をクローズ
C ---------------------------------------------------------------------
      CLOSE(11)
      CLOSE(12)
      CLOSE(13)
C =====================================================================
C     出力形式のフォーマット
C     ここでは、
C     整数4桁、1空白、実数8桁で小数点以下4桁、1空白、
C                                                実数6桁で小数点以下2桁
C ---------------------------------------------------------------------
  501 FORMAT(I4,1X,F8.4,1X,F6.2)
C =====================================================================
      STOP
      END
C*********************************************************************
C    相関係数を計算する副プログラム
C
C     相関係数(R)
C
C              Σ[ {X-(Xの平均)} * {Y-(Yの平均)} ]
C  R= -----------------------------------------------------------------
C    SQRT [ Σ{(X-Xの平均)*(X-Xの平均)} * Σ{(X-Xの平均)*(X-Xの平均)} ]
C
C*********************************************************************
      SUBROUTINE PEARSN(X,Y,N,R)
C      PARAMETER (TINY=1.E-20)
      DIMENSION X(N),Y(N)
C =====================================================================
      AX=0.
      AY=0.
C =====================================================================
C     配列X, Y の平均の計算
C ---------------------------------------------------------------------
      DO 10 J=1,N
        AX=AX+X(J)
        AY=AY+Y(J)
   10 CONTINUE
      AX=AX/N
      AY=AY/N
C =====================================================================
      SXX=0.
      SYY=0.
      SXY=0.
      DO 20 J=1,N
        XT=X(J)-AX   ! 平均からのズレ(XT)の計算(ここでは北極振動指数)
C                    ! XT=X-(Xの平均)
        YT=Y(J)-AY   ! 平均からのズレ(YT)の計算(ここでは気温)
C                    ! YT=Y-(Yの平均)
C
        SXX=SXX+XT**2  ! Σ[ {(X-(Xの平均)} * {(X-(Xの平均)} ] の計算
        SYY=SYY+YT**2  ! Σ[ {(Y-(Yの平均)} * {(Y-(Yの平均)} ] の計算
        SXY=SXY+XT*YT  ! Σ[ {(X-(Xの平均)} * {(Y-(Yの平均)} ] の計算
C
   20 CONTINUE
      R=SXY/SQRT(SXX*SYY)  ! 相関係数の値
C
      RETURN
      END
