Program TQ_ClimIDX_EOF_of_PDSI
! Chuong trinh tinh tuong quan giua cac EOF cua PDSI va ENSO

Use netcdf
Integer, Parameter :: Nreg=7, Yr1=1980, Yr2=2014, Nyear=Yr2-Yr1+1, Nmon=12, Ntime=Nyear*Nmon
Integer, Parameter :: N_ClimIDX=16
Real, Parameter    :: Miss=-99.0

!------------------------------------------
Character (Len = 11), Dimension(N_ClimIDX) :: CL_Name

Real, Allocatable                  :: X(:), Y(:), Corr(:,:), Prob(:,:), Lon(:), Lat(:)
Real, Allocatable                  :: Corr1(:,:), Prob1(:,:)
Real, Allocatable                  :: Corr2(:,:), Prob2(:,:)

Real, Allocatable                  :: EOF_PDSI(:,:), EOF_ClimIDX(:,:)
Real, Allocatable, Dimension (:,:) :: PDSI, Reg_PDSI, T2m, Prec   ! Monthly Temperature and Precipitation 
Real, Allocatable, Dimension (:)   :: Reg

Real, Dimension (N_ClimIDX, Ntime) :: ClimIDX
Real :: YY(Nyear)
Character (Len=255) :: Res_Dir, Fname, Varn, Var_ln, Var_unit, lon_n, Clat_n
Integer :: ist, Nst, NPCA, NPAC1, Nt, it, ipca, iidx
Real :: Corre, Proba


Res_Dir="/home/tanpv/SEACLID_New/RES_Stn/"

  Fname="/home/tanpv/SEACLID_New/RES/PDSI_Res/PDSI_Stn_No_Rain_Center.nc"
  Call Read_PDSI_Stn (Fname)     ! Read PDSI, T2m, Prec(Nst,Ntime), and Reg(Nst)
  Nst=131
  Where (PDSI<-20.0) PDSI=-20.0

  Allocate (Reg_PDSI(1:Ntime,1:Nreg), X(Nst), Y(Nreg))
  Do it=1, Ntime
     X(1:Nst)=PDSI(1:Nst,it); Call Regional_Mean (X, Y)
     Reg_PDSI(it,1:Nreg) = Y(1:Nreg)
     print "(I4,10F7.3)", it,Reg_PDSI(it,1:Nreg)
  Enddo
  Deallocate (X, Y)  

pause
Call Read_Clim_IDX              ! Return: ClimIDX, CL_Name (N_ClimIDX)
!print*,ClimIDX; pause

! Reading EOF of PDSI from file PCA_PDSI.nc
Call rd_ncdf_2d (trim(Res_Dir)//"PCA_PDSI.nc", Lon, Lat, Nt, NPCA, EOF_PDSI, "PCA")
Deallocate (Lon, Lat);

! Reading EOF of ClimIDX from file PCA_PDSI.nc
Call rd_ncdf_2d (trim(Res_Dir)//"PCA_ClimIDX.nc", Lon, Lat, Nt, NPCA1, EOF_ClimIDX,"PCA")

print*, "NPCA, Nt=", NPCA, Nt
pause

!------------------------------------------
Open (3, file="TQ_EOF_PDSI_Clim.txt", Status="Unknown")

Allocate (Corr(NPCA,N_ClimIDX), Prob(NPCA,N_ClimIDX), X(Nt), Y(Nt)) 

write (3,*)" He so tuong quan (va P-value tuong ung) giua cac EOF cua PDSI (cot) voi cac ClimIDX (hang)"
write (3,"(A11,5(2I7,2x))") "Cl_Name", ((ipca,ipca), ipca=1, NPCA)
Do iidx=1, N_ClimIDX
   X(1:Nt)=ClimIDX(iidx,1:Nt)
   Do ipca=1, NPCA
      Y(1:Nt)=EOF_PDSI(1:Nt,ipca)
      Call Correlation (X, Y, Nt, Corre, Proba, Miss)
      Corr(ipca,iidx)=Corre; Prob(ipca,iidx)=Proba
   Enddo
   write (3,"(A11,5(2F7.3,2x))") CL_Name(iidx), (Corr(ipca,iidx),Prob(ipca,iidx), ipca=1, NPCA)
   print "(A,5(2F7.3, 2x))", CL_Name(iidx), (Corr(ipca,iidx),Prob(ipca,iidx), ipca=1, NPCA)
Enddo
!------------------------------------------

Allocate (Corr1(NPCA,NPCA1), Prob1(NPCA,NPCA1))

write (3,*)" He so tuong quan (va P-value tuong ung) giua cac EOF cua PDSI (cot) voi cac EOF cua ClimIDX (hang)"

write (3,"(A11,5(2I7,2x))") "EOF_CL", ((ipca,ipca), ipca=1, NPCA)
Do iidx=1, NPCA1
   X(1:Nt)=EOF_ClimIDX(1:Nt,iidx)
   Do ipca=1, NPCA
      Y(1:Nt)=EOF_PDSI(1:Nt,ipca)
      Call Correlation (X, Y, Nt, Corre, Proba, Miss)
      Corr1(ipca,iidx)=Corre; Prob1(ipca,iidx)=Proba
   Enddo
   write (3,"(A9,I2.2,5(2F7.3,2x))") "EOF",iidx, (Corr1(ipca,iidx),Prob1(ipca,iidx), ipca=1, NPCA)
   print "(A,5(2F7.3, 2x))", "EOF-EOF", (Corr1(ipca,iidx),Prob1(ipca,iidx),ipca=1, NPCA)
Enddo
!------------------------------------------
  Deallocate (X, Y)
  Allocate (Corr2(N_ClimIDX,Nreg), Prob2(N_ClimIDX,Nreg), X(Ntime), Y(Ntime))
  write (3,*)" He so tuong quan (va P-value tuong ung) giua cac ClimIDX (hang) voi PDSI trung binh vung (cot)"
  write (3,"(A11,10(2I7,2x))") "IDX\Reg", ((ipca,ipca), ipca=1, Nreg)
  Do iidx=1, N_ClimIDX
     X(1:Ntime)=ClimIDX(iidx,1:Ntime)
     Do ipca=1, Nreg
        Y(1:Ntime)=Reg_PDSI(1:Ntime,ipca)
        Call Correlation (X, Y, Ntime, Corre, Proba, Miss)
        Corr2(ipca,iidx)=Corre; Prob2(ipca,iidx)=Proba
     Enddo
     write (3,"(A11,10(2F7.3,2x))") CL_Name(iidx), (Corr2(ipca,iidx),Prob2(ipca,iidx), ipca=1, Nreg)
  Enddo

  Close (3)
!------------------------------------------------------------------------------------------
CONTAINS
!------------------------------------------
Subroutine Regional_Mean (X, Y)
Integer :: ifile, ist, ireg, Nc
Real, Dimension (Nst)  :: X
Real, dimension (Nreg) :: Y
Real :: tmp
   Do ireg=1, Nreg ! +1
     ! if (ireg>=1.AND.ireg<=Nreg) then
         tmp=0; Nc=0
         Do ist=1, Nst
               if (Int(Reg(ist))==ireg) then
                  if (X(ist)/=Miss.AND.X(ist)/=0.) then; Nc=Nc+1; tmp=tmp+X(ist); endif
               endif
         Enddo
         if (Nc>1) then; Y(ireg)=tmp/Real(Nc); else; Y(ireg)=Miss; endif
     ! else if (ireg==Nreg+1) then
     !    tmp=0; Nc=0
     !    Do ist=1, Nst
     !          if (Int(Reg(ist))>=1.AND.Int(Reg(ist))<=Nreg) then
     !             if (X(ist)/=Miss.AND.X(ist)/=0.) then; Nc=Nc+1; tmp=tmp+X(ist); endif
     !          endif
     !    Enddo
     !    if (Nc>1) then; Y(ireg)=tmp/Real(Nc); else; Y(ireg)=Miss; endif
     ! endif
   Enddo
End Subroutine Regional_Mean
!------------------------------------------

Include "SUB/Read_Clim_IDX.f90"
Include "SUB/Read_PDSI_Stn.f90"
Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"

Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"
!Include "../Libs/Student_Prob.f90"

!------------------------------------------------------------------------------------------
END
Include "/home/tanpv/Mylibs/Student_Test.f90"
Include "/home/tanpv/Mylibs/Special_Function.f90"
!Include "../Libs/Student_Prob.f90"
!------------------------------------------------------------------------------------------