! Tinh tuong quan tre giua NTC hoac TCday voi cac chi so khi hau
    Use NetCDF
    Real, Parameter :: Miss=-99.0
    !Integer, Parameter :: Yr1=1970, Yr2=2017, Nyear=Yr2-Yr1+1, Nmon=12, Nidx=16, NTC_Cha=6
    Integer, Parameter :: Yr1=1979, Yr2=2018, Nyear=Yr2-Yr1+1, Nmon=12, Nidx=30, NTC_Cha=6
    !Character (Len=8), Parameter :: Clim_nm(Nidx)=(/"ONI","MEI","REQSOI","SOI","RINDO","Nino34", &
    !                                "CPAC850","BEST","WPC850","PDO","PNA","DMI","WHWP","TNI","WP","QBO"/)

    Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "oni", "tni", &
                                                     "mei", "pna", "censo", "ammsst", "CAR_ersst", "NTA_ersst", &
                                                     "soi", "reqsoi", "whwp", "amonus", "np", "repac_slpa", &
                                                     "rindo_slpa", "tsa", "tna","zwnd20", "wpac85", "epac85", &
                                                     "cpac85", "qbo50", "qbo30", "wp", "ao", "nao"/)
    Character (Len=5),  Parameter :: TC_nm(NTC_Cha)=(/"tc","tc3m","tc6m","ntc","ntc3m","ntc6m"/)

    Real, Allocatable :: Corr_lag(:,:,:,:), Prob(:,:,:,:) ! (NTC_Cha, Nmon, Nidx, 0:Nmon)
    Real, Allocatable :: TC_Cha(:,:,:), ClimIDX(:,:,:)    ! NTC_Cha/Nidx,Nyear,Nmon
    Real, Allocatable :: Clim_Dat(:,:,:)                  ! Nyear-1,Nlag*Nmon,Nmon (Nlag=10)

    Character (Len=200) :: FnameTC, FnameIDX, Fname
    !--------------------------------------------------------------------------------------
    ! Tinh cac EOF cua cac ClimIDX va tuong quan giua cac EOF voi TC_Char
    Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:)
    Real*8, Allocatable :: EOFs(:,:), PCs(:,:)
    Real*8, Allocatable   :: Eigenvalues(:), ECumulate(:)
    Integer :: Nt, Nx   ! Kich thuoc thuc cua mang X
    Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K
    Real*8 :: Factor
    Integer, Parameter :: Nidx1=Nidx   ! =16-2

    Real, Allocatable :: EOF_Clim(:,:,:), ClimIDX1(:,:,:)    ! Loai bo CPAC850 va WPC850 vi SL ngan
    Integer :: imon, iyr, iidx, imod, Nlag
    Real, Allocatable :: Corr_all_lag(:,:,:),Prob_all_lag(:,:,:)  ! Nmon,Nmode,NTC_Cha
    !--------------------------------------------------------------------------------------
    FnameTC="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc"   ! Contains: ntc, and tcday for 1,3,6 mon
    FnameIDX="/home/tanpv/Data/ClimIDX_2020/Clim_IDX_2020.nc"

    Call Read_TC (FnameTC)
    Call Read_ClimIDX (FnameIDX)

    Fname="TQ_tre_TC_ClimIDX_2020.nc"
    Call Calc_Corre_lag (Fname)

      !---------------------------------------------------------------------------------
      switch=0; Standard=1; PC_Stard=0; Nmode=6
      Nt=Nyear; Nx=Nidx1

      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Factor = Sqrt(Real(Nt))

      Allocate (xdat(Nt,Nx), ydat(Nt,Nx))
      Allocate (Xtb(Nx), Sx(Nx))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (Eigenvalues(Mcount))          ! Tri rieng va Vecto rieng (1D)
      Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
      Allocate (EOF_Clim(Nmode,Nt,Nmon))
      !---------------------------------------------------------------------------------
      ! ClimIDX1(Nidx1,Nyear,Nmon)
    Do imon=1,Nmon
       Do iyr=1,Nt  !=Nyear
          Do iidx=1,Nx ! =Nidx1
             xdat(iyr,iidx)=ClimIDX1(iidx,iyr,imon)
          Enddo
       Enddo
       Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                     Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode)
       Do imod=1,Nmode
          Do iyr=1,Nt
             EOF_Clim(imod,iyr,imon)=EOFs(iyr,imod)
          Enddo
       Enddo
    Enddo
    Fname="TQ_tre_TC_EOF_Clim_2020.nc"
    Call Calc_Corre_lag_eof (Fname)
    Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate,EOF_Clim)
    !--------------------------------------------------------------------------------------
    ! Clim_Dat(1:Nyear-1,Nidx1*Nlag,imon)

      Call Multi_lag   ! Return: Clim_Dat
      print*,"Done"; pause

      switch=0; Standard=1; PC_Stard=0; Nmode=6
      Nt=Nyear-1; Nlag=Nmon-2; Nx=Nidx1*Nlag

      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Factor = Sqrt(Real(Nt))

      Allocate (xdat(Nt,Nx), ydat(Nt,Nx))
      Allocate (Xtb(Nx), Sx(Nx))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (Eigenvalues(Mcount))          ! Tri rieng va Vecto rieng (1D)
      Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
      Allocate (EOF_Clim(Nmode,Nt,Nmon))
      
    Do imon=1,Nmon
       print*,"=============imon==>",imon
       Do iyr=1,Nt  !=Nyear
          Do iidx=1,Nx ! =Nidx1
             xdat(iyr,iidx)=Clim_Dat(iyr,iidx,imon)
          Enddo
       Enddo
       Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                     Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode)
       Do imod=1,Nmode
          Do iyr=1,Nt
             EOF_Clim(imod,iyr,imon)=EOFs(iyr,imod)
          Enddo
       Enddo
    Enddo
    Fname="TQ_TC_All_EOF_Clim_2020.nc"
    Call Corre_EOF_All_TC (Fname)      
    !--------------------------------------------------------------------------------------
    CONTAINS
    !--------------------------------------------------------------------------------------
    Subroutine Corre_EOF_All_TC (Fname)
    Integer :: ireg, imon, k, lag, imon1, imon11, Nxy
    Real :: Corr, Probe
    Character (Len=*) :: Fname
    Character (Len=100) :: var_nm, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
    Real, Allocatable :: X(:), Y(:)
    Real :: xmon(Nmon), xtc(NTC_Cha), xidx(Nmode)

    Allocate (Corr_all_lag(Nmon,Nmode,NTC_Cha),Prob_all_lag(Nmon,Nmode,NTC_Cha))
    Allocate(X(Nyear), Y(Nyear))

    Do ist=1, NTC_Cha
       Do iidx=1, Nmode
          Do imon=1, Nmon
             X(1:Nt)=EOF_Clim(iidx,1:Nt,imon); Y(1:Nt)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1
             Call Correlation (X, Y, Nxy, Corr, Probe, Miss)
             Corr_all_lag(imon,iidx,ist)=Corr; Prob_all_lag(imon,iidx,ist)=Probe
          Enddo
       Enddo
    Enddo

    Do imon=1,Nmon; xmon(imon)=imon; Enddo
    Do itc=1,NTC_Cha; xtc(itc)=itc; Enddo
    Do iidx=1,Nmode; xidx(iidx)=iidx; Enddo
    var_nm="corr"; var_ln="Lag Correlation between (n)tc,(n)tc3m,(n)tc6m and EOFs of Climate indices"
    var_unit="-"; t_unit="months since 1971-01-01 "; lon_nm="mon"; lat_nm="pc"; t_nm="tc"

    Call wr_ncdf_3d (Fname,xmon,xidx,xtc,Nmon,Nmode,NTC_Cha,Corr_all_lag,var_nm,t_unit, var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
    var_nm="prob"; var_ln="P-value"
    Call wr_ncdf_3d_apd (Fname, Prob_all_lag,var_nm,Nmon,Nmode,NTC_Cha, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

    End Subroutine Corre_EOF_All_TC
    !--------------------------------------------------------------------------------------
    Subroutine Read_TC (Fname)
       Character (Len=*) :: Fname
       Real, Allocatable :: Year(:),Month(:), tmp(:,:)
       Integer :: Nyr, Nmn, icha, iyr, i,j

       Allocate (TC_Cha(NTC_Cha,Nyear,Nmon))
       Do icha=1,NTC_Cha
          Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(TC_nm(icha)))
          j=0
          Do iyr=Yr1,Yr2
             Do i=1,Nyr
                if (Year(i)==iyr) then
                   j=j+1; TC_Cha(icha,j,:)=tmp(i,:)
                endif
             Enddo
          Enddo
          Deallocate (Year,Month,tmp)
       Enddo
    End Subroutine Read_TC
    !--------------------------------------------------------------------------------------
    Subroutine Read_ClimIDX (Fname)
       Character (Len=*) :: Fname
       Real, Allocatable :: Year(:),Month(:), tmp(:,:)
       Integer :: Nyr, Nmn, iidx, iyr, i,j, k

       Allocate (ClimIDX(Nidx,Nyear,Nmon),ClimIDX1(Nidx1,Nyear,Nmon))
       k=0   ! De chon cac IDX cho ClimIDX1
       Do iidx=1,Nidx
          Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx)))
          j=0
          Do iyr=Yr1,Yr2
             Do i=1,Nyr
                if (Year(i)==iyr) then
                   j=j+1; ClimIDX(iidx,j,:)=tmp(i,:)
          !         print "(F5.0,12F8.3)",Year(i),ClimIDX(iidx,j,:)
                endif
             Enddo
          Enddo
          if (trim(Clim_nm(iidx))/="CPAC850".AND.trim(Clim_nm(iidx))/="WPC850") then
             k=k+1;ClimIDX1(k,:,:)=ClimIDX(iidx,:,:)
          endif
          Deallocate (Year,Month,tmp)
          !print*,Clim_nm(iidx);pause
       Enddo
    End Subroutine Read_ClimIDX
    !--------------------------------------------------------------------------------------
    Subroutine Calc_Corre_lag_eof (Fname)

    Integer :: ireg, imon, k, lag, imon1, imon11, Nxy
    Real :: Corr, Probe
    Character (Len=*) :: Fname
    Character (Len=100) :: var_nm, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
    Real, Allocatable :: X(:), Y(:)
    Real :: xmon(Nmon), xlag(Nmon+1), xtc(NTC_Cha), xidx(Nmode)

    !Allocate(Corr_lag(NTC_Cha, Nmon, Nidx, 0:Nmon), Prob(NTC_Cha, Nmon, Nidx, 0:Nmon))
    Allocate (Corr_lag(Nmon,Nmon+1,NTC_Cha,Nmode),Prob(Nmon,Nmon+1,NTC_Cha,Nmode))
    Allocate(X(Nyear), Y(Nyear))

    Do ist=1, NTC_Cha
       Do iidx=1, Nmode
          Do imon=1, Nmon
             Do lag=0, Nmon
                if (lag==0) then
                   X(1:Nyear)=EOF_Clim(iidx,1:Nyear,imon); Y(1:Nyear)=TC_Cha(ist,1:Nyear,imon); Nxy=Nyear
                else
                   imon1=imon-lag
                   if (imon1>0) then
                      X(1:Nyear)=EOF_Clim(iidx,1:Nyear,imon1); Y(1:Nyear)=TC_Cha(ist,1:Nyear,imon); Nxy=Nyear
                   else
                      imon11 = 12 - ABS(imon1)
                      !X(1:Nyear-1)=EOF_Clim(iidx,1:Nyear-1,imon11); Y(1:Nyear-1)=TC_Cha(ist,1:Nyear-1,imon); Nxy=Nyear-1
                      X(1:Nyear-1)=EOF_Clim(iidx,1:Nyear-1,imon11); Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1
                   endif
                endif
                Call Correlation (X, Y, Nxy, Corr, Probe, Miss)
                !Corr_lag(ist,imon,iidx,lag)=Corr; Prob(ist,imon,iidx,lag)=Probe
                Corr_lag(imon,lag+1,ist,iidx)=Corr; Prob(imon,lag+1,ist,iidx)=Probe

    !            Call Corre (X, Y, Nxy, Corr_lag(ist,imon,iidx,lag), Prob(ist,imon,iidx,lag), Miss)
    !            print*, ist, imon,iidx,lag, Corr(ist,imon,iidx,lag), Prob(ist,imon,iidx,lag)
             Enddo
          Enddo
       Enddo
    Enddo
    Deallocate(X, Y)

    Do imon=1,Nmon; xmon(imon)=imon; Enddo
    Do lag=1,Nmon+1; xlag(lag)=lag-1; Enddo
    Do itc=1,NTC_Cha; xtc(itc)=itc; Enddo
    Do iidx=1,Nmode; xidx(iidx)=iidx; Enddo
    var_nm="corr"; var_ln="Lag Correlation between (n)tc,(n)tc3m,(n)tc6m and EOFs of Climate indices"
    var_unit="-"; t_unit="months since 1971-01-01 "; lon_nm="mon"; lat_nm="lag"; z_nm="tc"; t_nm="idx"
    Call wr_ncdf_4d (Fname,xmon,xlag,xtc,xidx,Nmon,Nmon+1,NTC_Cha,Nmode, Corr_lag, var_nm, &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
    var_nm="prob"; var_ln="P-value"
    Call wr_ncdf_4d_apd (Fname, Prob, var_nm, Nmon,Nmon+1,NTC_Cha,Nmode, &
                           var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
    Deallocate (Corr_lag,Prob)
    !do lag=0, Nmon
    !   do iidx=1, Nidx
    !      print*,lag,iidx
    !      print "(12F7.3)", ((Corr_av(ist,imon,iidx,lag),imon=1,Nmon),ist=1,Nreg)
    !   enddo
    !   pause
    !enddo
    End Subroutine Calc_Corre_lag_eof
    !--------------------------------------------------------------------------------------
    Subroutine Calc_Corre_lag (Fname)

    Integer :: ireg, imon, k, lag, imon1, imon11, Nxy
    Real :: Corr, Probe
    Character (Len=*) :: Fname
    Character (Len=100) :: var_nm, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
    Real, Allocatable :: X(:), Y(:)
    Real :: xmon(Nmon), xlag(Nmon+1), xtc(NTC_Cha), xidx(Nidx)

    !Allocate(Corr_lag(NTC_Cha, Nmon, Nidx, 0:Nmon), Prob(NTC_Cha, Nmon, Nidx, 0:Nmon))
    Allocate (Corr_lag(Nmon,Nmon+1,NTC_Cha,Nidx),Prob(Nmon,Nmon+1,NTC_Cha,Nidx))
    Allocate(X(Nyear), Y(Nyear))

    Do ist=1, NTC_Cha
       Do iidx=1, Nidx
          Do imon=1, Nmon
             Do lag=0, Nmon
                if (lag==0) then
                   X(1:Nyear)=ClimIDX(iidx,1:Nyear,imon); Y(1:Nyear)=TC_Cha(ist,1:Nyear,imon); Nxy=Nyear
                else
                   imon1=imon-lag
                   if (imon1>0) then
                      X(1:Nyear)=ClimIDX(iidx,1:Nyear,imon1); Y(1:Nyear)=TC_Cha(ist,1:Nyear,imon); Nxy=Nyear
                   else
                      imon11 = 12 - ABS(imon1)
                      !X(1:Nyear-1)=ClimIDX(iidx,1:Nyear-1,imon11); Y(1:Nyear-1)=TC_Cha(ist,1:Nyear-1,imon); Nxy=Nyear-1
                      X(1:Nyear-1)=ClimIDX(iidx,1:Nyear-1,imon11); Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1
                   endif
                endif
                Call Correlation (X, Y, Nxy, Corr, Probe, Miss)
                !Corr_lag(ist,imon,iidx,lag)=Corr; Prob(ist,imon,iidx,lag)=Probe
                Corr_lag(imon,lag+1,ist,iidx)=Corr; Prob(imon,lag+1,ist,iidx)=Probe

    !            Call Corre (X, Y, Nxy, Corr_lag(ist,imon,iidx,lag), Prob(ist,imon,iidx,lag), Miss)
    !            print*, ist, imon,iidx,lag, Corr(ist,imon,iidx,lag), Prob(ist,imon,iidx,lag)
             Enddo
          Enddo
       Enddo
    Enddo
    Deallocate(X, Y)

    Do imon=1,Nmon; xmon(imon)=imon; Enddo
    Do lag=1,Nmon+1; xlag(lag)=lag-1; Enddo
    Do itc=1,NTC_Cha; xtc(itc)=itc; Enddo
    Do iidx=1,Nidx; xidx(iidx)=iidx; Enddo
    var_nm="corr"; var_ln="Lag Correlation between (n)tc,(n)tc3m,(n)tc6m and Climate indices"
    var_unit="-"; t_unit="months since 1971-01-01 "; lon_nm="mon"; lat_nm="lag"; z_nm="tc"; t_nm="idx"
    Call wr_ncdf_4d (Fname,xmon,xlag,xtc,xidx,Nmon,Nmon+1,NTC_Cha,Nidx, Corr_lag, var_nm, &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
    var_nm="prob"; var_ln="P-value"
    Call wr_ncdf_4d_apd (Fname, Prob, var_nm, Nmon,Nmon+1,NTC_Cha,Nidx, &
                           var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
    Deallocate (Corr_lag,Prob)
    End Subroutine Calc_Corre_lag 
    !--------------------------------------------------------------------------------------
    Subroutine Multi_lag
       Integer :: imon, iidx, lag, imon1, imon11, Nlag, k
       Real, Allocatable :: X(:), Y(:)

    Nlag=Nmon-2  
    Allocate (Clim_Dat(Nyear-1,Nidx1*Nlag,Nmon))
    Allocate(X(Nyear-1), Y(Nyear-1))
    !Do ist=1, NTC_Cha
       Do imon=1, Nmon    ! iidx=1, Nidx
          k=0
          Do iidx=1, Nidx1 ! imon=1, Nmon
             Do lag=3,Nmon ! 0, Nmon -> Lay tu lag=3-lag=12 -> 10 lags
                !if (lag==0) then
                !   X(1:Nyear)=ClimIDX(iidx,1:Nyear,imon); Y(1:Nyear)=TC_Cha(ist,1:Nyear,imon); Nxy=Nyear
                !else
                   imon1=imon-lag
                   if (imon1>0) then
                      X(1:Nyear-1)=ClimIDX1(iidx,2:Nyear,imon1); !Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1
                   else
                      imon11 = 12 - ABS(imon1)
                      !X(1:Nyear-1)=ClimIDX(iidx,1:Nyear-1,imon11); Y(1:Nyear-1)=TC_Cha(ist,1:Nyear-1,imon); Nxy=Nyear-1
                      X(1:Nyear-1)=ClimIDX1(iidx,1:Nyear-1,imon11); !Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1
                   endif
                !endif
                k=k+1; Clim_Dat(1:Nyear-1,k,imon)=X(1:Nyear-1)
                print*,"k=",k,Nidx1*Nlag
             Enddo
          Enddo
       Enddo
    !Enddo
    Deallocate(X, Y)

    End Subroutine Multi_lag
    !--------------------------------------------------------------------------------------

    Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
    Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
    Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90"
    Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"

    !--------------------------------------------------------------------------------------
    END

    Include "/home/tanpv/Mylibs/Student_Prob.f90"
    !--------------------------------------------------------------------------------------