! Khao sat moi quan he giua cac chi so khi hau va cac chi so han
use netcdf
Implicit None
Real, Parameter :: Miss = -99.
Integer, Parameter :: Nreg = 7, Nmon = 12, Yr1 = 1980, Yr2 = 2014, Nyear = Yr2-Yr1+1
Integer, Parameter :: NdrIDX = 10   ! PDSI, scPDSI, SPI01-12, SPEI01-12
Integer, Parameter :: NClimIDX = 16 
Real, Allocatable, Dimension (:,:,:) :: DrIDX   ! (Nst, Ntime, NdrIDX)
Real, Allocatable, Dimension (:,:,:) :: DrIDX_Reg   ! (Nreg, Ntime, NdrIDX)
Real, Allocatable, Dimension (:,:)   :: ClimIDX ! (Ntime, NClimIDX)
Real, Allocatable, Dimension (:)     :: ENSO_Yrs ! (Ntime)
Real, Allocatable, Dimension (:,:)   :: PCA_Cli  ! (Ntime, Npca)
Real, Allocatable, Dimension (:)     :: Xlon, Xlat, Reg, Xtime
Logical :: status

Character (Len=6), Parameter :: DroIDX_Name(NdrIDX)  =(/"PDSI","scPDSI","SPI01","SPI03","SPI06","SPI12","SPEI01","SPEI03","SPEI06","SPEI12"/)
Character (Len=6), Parameter :: CliIDX_Name(NClimIDX)=(/"DMI","BEST",  "MEI",   "Nino34","ONI","SOI","WP",    "WPC850", &
                                                       "QBO","RINSLP","CPA850","PDO",   "TNI","PNA","REQSOI","WHWP"/)


Integer :: ivar, ist, it, Ntime, Nst, Npca
Character (Len=255) :: FName

Character (Len=50), Parameter :: Fname_idx(NClimIDX) = &
                   (/"DMI.txt","BEST.txt","MEI.txt","Nino3.4.txt","ONI.txt","SOI.txt","WP.txt","WPC850.txt", &
                     "QBO.txt","RINDO_SLPA.txt","CPAC850.txt","PDO.txt","TNI.txt","PNA.txt","REQSOI.txt","WHWP.txt"/)
!------------------------------------------------------------------------------
Real, Allocatable, Dimension (:,:,:) :: Corr, Prob   ! (Nst, NdrIDX, NClimIDX) 
Real, Allocatable, Dimension (:,:,:) :: Corr_Reg, Prob_Reg   ! (Nreg, NdrIDX, NClimIDX)
Real, Allocatable, Dimension (:,:,:) :: Corr_PCA, Prob_PCA   ! (Nst, NdrIDX, Npca) 
Real, Allocatable, Dimension (:,:,:) :: Corr_Reg_PCA, Prob_Reg_PCA   ! (Nreg, NdrIDX, Npca)

!------------------------------------------------------------------------------
!Allocate (ClimIDX(Nyear*Nmon,NClimIDX))
! Reading PDSI
   Fname="/home/tanpv/SEACLID_New/RES/PDSI_Res/PDSI_Stn_No_Rain_Center.nc";   Call Read_PDSI_nc (Fname)
! Reading SPI (SPI01, 12)
   Fname="/home/tanpv/SEACLID_New/RES/SPI_Res/SPI_Stn_No_Rain_Center.nc";     Call Read_SPI_nc (FName)
! Reading SPEI (SPEI06,12)
   Fname="/home/tanpv/SEACLID_New/RES/SPI_Res/SPEI_Stn_No_Rain_Center.nc";    Call Read_SPEI_nc (FName)
! Reading ClimIDX
   Call Read_Clim_IDX

   Call Cal_Drou_IDX_Reg
   Call Cal_Corr_Matrix
   Call Write_Corre_To_TXT
! print test
!   Call Print_test (ClimIDX, Ntime, NClimIDX)
CONTAINS
!------------------------------------------------------------------------------
Subroutine Print_test (X, N, M)
Integer :: i,j,N,M
Real :: X(N,M)
Do i=1, N
   print "(200F7.2)", (X(i,j), j=1,M)
Enddo
End Subroutine Print_test
!------------------------------------------------------------------------------
Subroutine Cal_Drou_IDX_Reg
Integer :: ist, idr, ireg, it, ic

Allocate (DrIDX_Reg(Nreg, Ntime, NdrIDX))
Do idr=1, NdrIDX
   Do it=1, Ntime
      Do ireg=1, Nreg
         DrIDX_Reg(ireg,it,idr)=0; ic=0
         Do ist=1, Nst
            if (Int(Reg(ist))==ireg) then
               if (DrIDX(ist,it,idr)/=Miss) then
                  DrIDX_Reg(ireg,it,idr)=DrIDX_Reg(ireg,it,idr)+DrIDX(ist,it,idr); ic=ic+1
         !         print*,Int(Reg(ist)), DrIDX(ist,it,idr), ic
               endif
            endif
         Enddo
         if (ic>0) then; DrIDX_Reg(ireg,it,idr)=DrIDX_Reg(ireg,it,idr)/Real(ic); else;DrIDX_Reg(ireg,it,idr)=Miss; endif
         print*,idr,it,ireg,DrIDX_Reg(ireg,it,idr)
      Enddo
   Enddo
Enddo
End Subroutine Cal_Drou_IDX_Reg
!------------------------------------------------------------------------------
Subroutine Write_Corre_To_TXT
Integer :: ist, ireg, idr, icl
Character (Len=6), Parameter :: PCA_Cli_Name(3)=(/"PCA1","PCA2","PCA3"/)
! Station
Open (3, file="Corr_Dr_IDX_Cli_IDX.txt", status="Unknown")
write(3, "(A)") "Correlation Coeffcients between Drought Indices and Climate Indices"
write(3,*)
Do icl=1, NClimIDX
   write(3,"(2A)") "Correlation Coeffcients between Drought Indices and ",CliIDX_Name(icl)
   write(3,"(A3,1x,2A8,20A8)") "Reg","Lon","Lat",(DroIDX_Name(idr),idr=1,NdrIDX)
   Do ist=1, Nst
      write(3,"(I3,1x,2F8.3,20F8.3)") Int(Reg(ist)), XLon(ist), XLat(ist), (Corr(ist,idr,icl),idr=1,NdrIDX)
!      write(3,"(12x,A8,20F8.3)") "p-value", (prob(ist,idr,icl),idr=1,NdrIDX)
   Enddo
   write(3,"(A)") "---------------------------------------------------"
Enddo
Close(3)
! Sub-region
Open (3, file="Corr_Dr_IDX_Cli_IDX_Reg.txt", status="Unknown")
write(3, "(A)") "Correlation Coeffcients between Drought Indices and Climate Indices"
write(3,*)
Do icl=1, NClimIDX
   write(3,"(2A)") "Correlation Coeffcients between Drought Indices and ",CliIDX_Name(icl)
   write(3,"(A3,1x,2A8,20A8)") "Reg",(DroIDX_Name(idr),idr=1,NdrIDX)
   Do ireg=1, Nreg
      write(3,"(I3,1x,20F8.3)") ireg, (Corr_Reg(ireg,idr,icl),idr=1,NdrIDX)
!      write(3,"(A3,20F8.3)") "p-v", (prob_Reg(ist,idr,icl),idr=1,NdrIDX)
   Enddo
   write(3,"(A)") "---------------------------------------------------"
Enddo
Close(3)
! Station
Open (3, file="Corr_Dr_IDX_Cli_IDX_PCA.txt", status="Unknown")
write(3, "(A)") "Correlation Coeffcients between Drought Indices and Climate Indices"
write(3,*)
Do icl=1, Npca
   write(3,"(2A)") "Correlation Coeffcients between Drought Indices and ",PCA_Cli_Name(icl)
   write(3,"(A3,1x,2A8,20A8)") "Reg","Lon","Lat",(DroIDX_Name(idr),idr=1,NdrIDX)
   Do ist=1, Nst
      write(3,"(I3,1x,2F8.3,20F8.3)") Int(Reg(ist)), XLon(ist), XLat(ist), (Corr_PCA(ist,idr,icl),idr=1,NdrIDX)
!      write(3,"(12x,A8,20F8.3)") "p-value", (prob_PCA(ist,idr,icl),idr=1,NdrIDX)
   Enddo
   write(3,"(A)") "---------------------------------------------------"
Enddo
Close(3)
! Sub-region
Open (3, file="Corr_Dr_IDX_Cli_IDX_Reg_PCA.txt", status="Unknown")
write(3, "(A)") "Correlation Coeffcients between Drought Indices and Climate Indices"
write(3,*)
Do icl=1, Npca
   write(3,"(2A)") "Correlation Coeffcients between Drought Indices and ",PCA_Cli_Name(icl)
   write(3,"(A3,1x,2A8,20A8)") "Reg",(DroIDX_Name(idr),idr=1,NdrIDX)
   Do ireg=1, Nreg
      write(3,"(I3,1x,20F8.3)") ireg, (Corr_Reg_PCA(ireg,idr,icl),idr=1,NdrIDX)
!      write(3,"(A3,20F8.3)") "p-v", (prob_Reg_PCA(ist,idr,icl),idr=1,NdrIDX)
   Enddo
   write(3,"(A)") "---------------------------------------------------"
Enddo
Close(3)
End Subroutine Write_Corre_To_TXT
!------------------------------------------------------------------------------
Subroutine Cal_Corr_Matrix
Integer :: ist, ireg, idr, icl
Real, Allocatable, Dimension(:) :: X, Y
Real :: Cor, Pro

! Correlation Drou_IDX & Cli_IDX
Allocate (Corr(Nst, NdrIDX, NClimIDX), Prob(Nst, NdrIDX, NClimIDX), X(Ntime), Y(Ntime))
Allocate (Corr_Reg(Nreg, NdrIDX, NClimIDX), Prob_Reg(Nreg, NdrIDX, NClimIDX))
Do icl=1, NClimIDX
   X(1:Ntime) = ClimIDX(1:Ntime, icl)
   Do idr=1, NdrIDX
      Do ist=1, Nst   ! For each station
         Y(1:Ntime) = DrIDX(ist,1:Ntime, idr); Call Corre (X, Y, Ntime, Cor, Pro, Miss)
         Corr(ist,idr,icl) = Cor; Prob(ist,idr,icl) = Pro
      Enddo
      Do ireg=1, Nreg  ! For each sub-region
         Y(1:Ntime) = DrIDX_Reg(ireg,1:Ntime, idr); Call Corre (X, Y, Ntime, Cor, Pro, Miss)
         Corr_Reg(ireg,idr,icl) = Cor; Prob_Reg(ireg,idr,icl) = Pro
      Enddo
   Enddo
Enddo
! Correlation Drou_IDX & PCA_Cli
Allocate (Corr_PCA(Nst, NdrIDX, Npca), Prob_PCA(Nst, NdrIDX, Npca))
Allocate (Corr_Reg_PCA(Nreg, NdrIDX, Npca), Prob_Reg_PCA(Nreg, NdrIDX, Npca))
Do icl=1, Npca
   X(1:Ntime) = PCA_Cli(1:Ntime, icl)
   Do idr=1, NdrIDX
      Do ist=1, Nst   ! For each station
         Y(1:Ntime) = DrIDX(ist,1:Ntime, idr); Call Corre (X, Y, Ntime, Cor, Pro, Miss)
         Corr_PCA(ist,idr,icl) = Cor; Prob_PCA(ist,idr,icl) = Pro
      Enddo
      Do ireg=1, Nreg  ! For each sub-region
         Y(1:Ntime) = DrIDX_Reg(ireg,1:Ntime, idr); Call Corre (X, Y, Ntime, Cor, Pro, Miss)
!         print*,"DrIDX_Reg=",Y; print*,"Corr=",Cor, Pro
         Corr_Reg_PCA(ireg,idr,icl) = Cor; Prob_Reg_PCA(ireg,idr,icl) = Pro
      Enddo
   Enddo
Enddo
End Subroutine Cal_Corr_Matrix
!------------------------------------------------------------------------------
Subroutine Corre (X, Y, N, Corr, Prob, Miss)
Integer :: N, i, j, N1
Real :: Corr, Prob, T_Dist
Real :: X(N), Y(N), tbx, tby, sx, sy, rxy, Miss, t
tbx=0.; tby=0.; sx=0.; sy=0.;  rxy=0.; n1=0
Do i=1, N
   if (x(i)/=Miss.And.y(i)/=Miss) then
      n1=n1+1
      tbx=tbx+x(i); tby=tby+y(i); sx=sx+x(i)*x(i); sy=sy+y(i)*y(i); rxy=rxy+x(i)*y(i)
   endif
Enddo
tbx=tbx/real(n1); tby=tby/real(n1); sx=sx/real(n1); sy=sy/real(n1); rxy=rxy/real(n1)
sx=sqrt(sx-tbx*tbx); sy=sqrt(sy-tby*tby); rxy=rxy-tbx*tby
Corr = rxy/(sx*sy)
if (N1<=2) then
   Corr = Miss; Prob = Miss
else
   t = ABS(Corr/ ( Sqrt(1-Corr**2) / Sqrt(Real(N1-2)) ) )
   Prob = T_Dist (t, Real(N1-2))
!   if (Prob<0.OR.Prob>0.5) print*,"T_Dist:", t, N1-2, Prob
endif
End Subroutine Corre
!------------------------------------------------------------------------------
Subroutine Read_PDSI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: tmplon, tmplat, tmpTime, XStn
Real, Allocatable, Dimension (:,:) :: PDSI, scPDSI
Integer :: Nlatmp, Nlotmp, Ntmp
Character (Len=50) :: VarName

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, PDSI, "ori_pdsi"); Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, scPDSI, "sel_pdsi"); Deallocate (XStn, tmpTime)
!
Allocate (DrIDX(Nst, Ntime, NdrIDX))
DrIDX(1:Nst,1:Ntime,1) = PDSI(1:Nst,1:Ntime); DrIDX(1:Nst,1:Ntime,2) = scPDSI(1:Nst,1:Ntime)
Deallocate (PDSI, scPDSI)

End Subroutine Read_PDSI_nc
!--------------------------------------------------------------------------------
Subroutine Read_SPI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: XStn, tmpTime
Real, Allocatable, Dimension (:,:) :: SPI
Integer ::  Nlon1, Nlat1, Nt

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_01")
DrIDX(1:Nst,1:Ntime,3) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI) 

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_03")
DrIDX(1:Nst,1:Ntime,4) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI) 

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_06")
DrIDX(1:Nst,1:Ntime,5) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI)

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_12")
DrIDX(1:Nst,1:Ntime,6) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI)

!print*,SPI12; pause
!Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlon, "lon"); Deallocate (XStn)
!Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlat, "lat"); Deallocate (XStn)
!Call rd_ncdf_1d (trim(FName), XStn, Nst, Reg,  "reg"); Deallocate (XStn)

End Subroutine Read_SPI_nc
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
Subroutine Read_SPEI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: XStn, tmpTime
Real, Allocatable, Dimension (:,:) :: SPI
Integer ::  Nlon1, Nlat1, Nt

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_01")
DrIDX(1:Nst,1:Ntime,7) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI)

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_03")
DrIDX(1:Nst,1:Ntime,8) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI)

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_06")
DrIDX(1:Nst,1:Ntime,9) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI)

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI, "SPI_12")
DrIDX(1:Nst,1:Ntime,10) = SPI(1:Nst,1:Ntime); Deallocate (XStn, tmpTime, SPI)

Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlon, "lon"); Deallocate (XStn)
Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlat, "lat"); Deallocate (XStn)
Call rd_ncdf_1d (trim(FName), XStn, Nst, Reg,  "reg"); Deallocate (XStn)

End Subroutine Read_SPEI_nc
!----------------------------------------------------------------------------------!
Subroutine Read_Clim_IDX
Integer :: iidx, iyr, imon, year, it1, it2, IO, icl, tt, ym
Character (Len=50) :: Fname
Real :: X(Nmon)
Logical :: Exists
!Allocate(ClimIDX(Nidx, Nyear, Nmon))
Allocate (ClimIDX(Nyear*Nmon,NClimIDX))
Do iidx=1, NClimIDX
   Fname="../../Nafosted/Data/ClimIDX/"//trim(Fname_idx(iidx))
   Inquire(file=trim(Fname), Exist=Exists)
   if (.Not.Exists) then
      print*,"File ",trim(Fname)," does not exist"
      cycle
   endif
   Open(1, file=trim(Fname), Status="Old")
   read(1,*)         ! Bo dong dau
   Do
      read(1,*,IOStat=IO) year, (X(imon),imon=1,Nmon)
      if (IO /= 0) exit
      Do iyr=Yr1, Yr2
         if (year==iyr) then
            it1=(iyr-Yr1)*Nmon+1; it2=it1+Nmon-1; ClimIDX(it1:it2,iidx)=X(1:Nmon)
         endif
      Enddo
   Enddo
   Close (1)
Enddo
! Read ENSO years
Allocate(ENSO_Yrs(Nyear*Nmon))
Fname="../../Nafosted/Data/ClimIDX/ENSO_Years.txt";
Inquire(file=trim(Fname), Exist=Exists)
   if (.Not.Exists) then
      print*,"File ",trim(Fname)," does not exist"
   else
      print*,"Reading ENSO..."
   endif
Open(1, file=trim(Fname), Status="Old")
   read(1,*)         ! Bo dong dau
   Do
      read(1,*,IOStat=IO) year, (X(imon),imon=1,Nmon)
!print "(A,12F6.1)","=====>",(X(imon),imon=1,Nmon)
      if (IO /= 0) exit
      Do iyr=Yr1, Yr2
         if (year==iyr) then
             it1=(iyr-Yr1)*Nmon+1; it2=it1+Nmon-1; ENSO_Yrs(it1:it2)=X(1:Nmon)
         endif
      Enddo
   Enddo
   Close (1)
! Reading PCA_Clim
Npca=3   ! 3 thanh phan chinh
Allocate(PCA_Cli(Ntime, Npca)); Open (1, file="PCA_Cli.txt")
Do icl=1, NClimIDX+2; read(1,*); Enddo ! Bo NClimIDX+2 dong dau
Do it1=1, Ntime
   read(1,*) tt, ym, (PCA_Cli(it1,icl),icl=1, Npca)
Enddo
Close(1)
End Subroutine Read_Clim_IDX
!------------------------------------------------------------------------------

Include "../Libs/sub_rd_ncdf_Tan.F90"
!Include "../Libs/sub_rd_ncdf.F90"
Include "../Libs/sub_wr_ncdf_Tan.F90"
!Include "../Sub_Trend_MK_Linear.f90"
!Include "Num_Dr_Spells.f90"
Include "../Libs/Sub_Write_GrADS_Stn.f90"
!----------------------------------------------------------------------------------!
END
!----------------------------------------------------------------------------------!
Include "/home/tanpv/SEACLID/Libs/Student_Prob.f90"
!Include "ENSO_FA.f90"
!----------------------------------------------------------------------------------!