! Chuong trinh nay dung de tinh cac he so tuong quan giua R (hoac T) voi cac yeu to ENSO
! Cac chuoi R (T) de tinh HSTQ co the lay tre thoi gian so voi ENSO 0, 1, 2,...12 thang
! Ky hieu Tau la do tre do
! Cac HSTQ duoc tinh cho tung thang (1..12)
! Tap hop cac HSTQ lap thanh ma tran Rxy(0:12, 1:12) voi so thu tu dong la do tre, so TT cot la thang
! 
! Cau truc file so lieu:
! - File so lieu Mua:
!   Dong 1: Ten tram
!   Dong 2: 3 so nguyen tuong ung la So nam co SL (N), Nam dau tien (y11), Nam cuoi cung (y12)
!   Tu dong thu 3 tro di, moi dong co 13 so, so dau tien la thu tu nam, 12 so tiep theo
!                         la Tong luong mua tung thang trong nam, tu thang 1 den thang 12
! - File so lieu ENSO:
!   Dong 1: 3 so nguyen tuong ung la So nam co SL (N), Nam dau tien (y21), Nam cuoi cung (y22)
!   Tu dong thu 2 tro di, moi dong co 13 so, so dau tien la thu tu nam, 12 so tiep theo
!                         la Tong luong mua tung thang trong nam, tu thang 1 den thang 12

Real, Allocatable :: SST(:,:), R(:,:) ! So lieu ban dau
Real, Allocatable :: X(:), Y(:)       ! Mang lam viec tam thoi
!Real :: X(100), Y(100)               ! Mang lam viec tam thoi
Real :: Rxy(0:12,1:12)                ! Mang cac he so tuong quan (0:12 - do tre; 1:12 - thang)
Integer :: y11, y12, y21, y22         ! Nam dau, nam cuoi cua cac chuoi SL
Integer :: y1, y2
Integer :: Tau, Mon, Mon1, yr, yr1, N, Nxy
Character (Len = 50) Fin1, Fin2, Fout ! Ten file
Real :: Rmiss_Value, Emiss_Value      ! Gia tri khuyet cua R va ENSO

Rmiss_Value = -99.0 ; Emiss_Value = -99.99   ! Co the thay the cho phu hop voi so lieu
Fin1 = 'R_Da_Nang.txt'                       ! Ten file so lieu Mua
Fin2 = 'SST_NINO34.txt'                      ! Ten file so lieu ENSO
Fout = 'HSTQ.txt'                            ! Ten file he so tuong quan

Open (1, file=Fin1); Open (2,file=Fin2)
Read (1,*)                                   ! Bo qua dong Ten tram
Read (1,*) N, y11, y12 ; Read (2,*) N, y21, y22

Allocate(R(y11:y12,1:12), SST(y21:y22, 1:12))  ! Khoi tao bo nho cho cac mang

Do i=y11, y12
   Read (1,*) N, (R(i,j), j=1,12)             ! So lieu mua
Enddo
Do i=y21, y22
   Read (2,*) N, (SST(i,j), j=1,12)           ! So lieu SST
Enddo

y1 = MAX (y11, y21) ;  y2 = MIN (y12, y22)    ! Xac lap chuoi thoi gian chung cho ca R va SST
N = y2 - y1 + 1                               ! Do dai chuoi

Allocate (X(N), Y(N))

Do Mon = 1,12
   Do Tau = 0,12
      Mon1 = Mon - Tau
      if (Mon1 > 0) then
         k = 0
         Do i=y1, y2   ! SST
            k = k + 1
            X(k) = SST(i,Mon1)
            Y(k) =   R(i,Mon)
         Enddo    
         Nxy = k
      Else     
         Mon11 = 12 - ABS(Mon1)
         k = 0
         Do i=y1, y2 - 1   ! SST
            k = k + 1
            X(k) = SST(i,Mon11)
            Y(k) = R(i+1,Mon)
         Enddo    
         Nxy = k
     endif         
      Rxy(Tau,Mon) = HSTQ
      print*, Mon, Tau, Nxy
   Enddo
Enddo

Open (3,file = Fout)
Do Tau = 0,12
   Write(3,'(I4,12F8.4)') Tau, (Rxy(Tau, Mon), Mon = 1, 12)
Enddo

CONTAINS
   Function HSTQ 
   Integer :: i, k, N1
   Real :: TBx, TBy, Sx, Sy, RRxy
   k = 0
   TBx = 0; TBy = 0; Sx = 0; Sy = 0; RRxy = 0;
   Do i = 1, Nxy
      if ((X(i) /= Emiss_Value).AND.(Y(i) /= Rmiss_Value)) then
         k = k + 1
         TBx = TBx + X(i); TBy = TBy + Y(i);
         Sx = Sx + X(i)*X(i) ; Sy = Sy + Y(i)*Y(i);
         RRxy = RRxy + X(i)*Y(i)
      endif     
   Enddo
   TBx = TBx/Real(k); TBy = TBy/Real(k); Sx = Sx/Real(k)-TBx*TBx; Sy = Sy/Real(k)-TBy*TBy; RRxy = RRxy/Real(k);
   HSTQ = (RRxy - TBx*TBy)/SQRT(Sx*Sy)
   End Function HSTQ

End