! Maximum Covariance Analysis of Drought Inde:wqx

      Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Yr1=1980, Yr2=2020, Nidx=20
      Integer, Parameter :: Nmode=6, Ntime = (Yr2-Yr1+1)*Nmon
      !Real, Parameter :: Slon1=70, Slon2=280, Slat1=-20, Slat2=40
      Real, Parameter :: Slon1=70, Slon2=290, Slat1=-20, Slat2=40

      Real, Allocatable :: Lon(:), Lat(:), xtime(:), Regs(:)
      Real, Allocatable :: Lon_sst(:), Lat_sst(:), Mask_sst(:,:)
      Real, Allocatable :: Mask(:,:), xlon(:), xlat(:), ylon(:), ylat(:)
      Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, i,j
      Integer :: SeaPoint, Nsst, Nlon_sst, Nlat_sst
      Integer :: Yr_sta
      Real :: Miss_sst
      !---------------------------------------------------------------- 
      Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ &
              "oni", "nino34", "nino3", "nino4", "mei", "rindo_slpa", &
              "cpac85", "soi", "reqsoi", "repac_slpa", &
              "wpac85", "pdo", "dmi", "tni", "pna", "whwp", "wp", &
              "qbo30","tna", "amonus"/)

      ! For MCA:
      Real, Allocatable :: Xbar(:), Sx(:), Ybar(:), Sy(:)
      Real, Allocatable :: X(:,:), Y(:,:), Xc(:,:), Yc(:,:)
      Real, Allocatable :: C(:,:), U(:,:), V(:,:), WW(:,:), W(:)
      Real, Allocatable :: A(:,:), B(:,:)   ! Timeseries of X & Y
      Real, Allocatable :: Uk(:,:), Vk(:,:), WWk(:,:), Wk(:)
      Real, Allocatable :: Ak(:,:), Bk(:,:), Corre(:), Probs(:)
      Real, Allocatable :: Accum(:), Percent(:)         ! Accumulated percentage of Singular values
      Real, Allocatable :: XX(:), YY(:), Dr_IDX(:,:,:)  ! (Nlon,Nlat,Ntime)
      Real, Allocatable :: U3D(:,:,:),V3D(:,:,:)        ! Uk,Vk in 3D array
      Real, Allocatable :: Homo(:,:),Hete(:,:)      ! Homo and Hetero maps (Nlon,Nlat,Nmode)
      Real, Allocatable :: p_Ho(:,:),p_He(:,:)      ! Pval

      Integer :: Nt, Nx, Ny, isys, ilen 
      Integer :: Mov_Ave_Step                    ! Buoc truot trong Moving average                
      Real :: Sumx,Corr,prob

      Real, Allocatable :: xx1d(:), yy1d(:)
      Logical :: Calc_mtx, Normalize, SVD_MatLab, SDV_MatLab, IsAnomaly
      !--------------------------------------------------------------------

      Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab, Dir_Out
      Character (Len=15)  :: Dat_src="Stn" ! "Ens"
      Character (Len=15)  :: varn="HWD"
      Character (Len=3)   :: Normalize_Opt
      Character (Len=20), Allocatable :: Stname(:)
      !--------------------------------------------------------------------
      ! LAPACK:
      Integer :: M, N, LDA, LDU, LDVT, LWMAX, INFO, LWORK
      DOUBLE PRECISION, Allocatable :: AA(:,:), UU(:,:), VT(:,:), SS(:), WORK(:)
      EXTERNAL ::        DGESVD
      !--------------------------------------------------------------------
      Real, Dimension (Nidx,Nidx) :: Corr_cli, prob_cli
      Real :: X1d(Ntime), Y1d(Ntime)


      Fname="/work/projects/G09/SEACLID/ClimIDX/Clim_IDX_2022.nc"
      Call Read_ClimIDX (Fname) !,varn)

      print"(A11,20(A11,A1))","Var,",((trim(Clim_nm(j)),","),j=1,Nidx)
      Do i=1,Nidx
         X1d(1:Ntime) = Y(i,1:Ntime)
         Do j=1,Nidx
            Y1d(1:Ntime) = Y(j,1:Ntime)
            Call Correlation(X1d,Y1d,Ntime,Corr,prob,Miss)
            Corr_cli(i,j) = Corr
            prob_cli(i,j) = prob
         Enddo
         print"(A11,20(F11.3,A1))",trim(Clim_nm(i))//",",((Corr_cli(i,j),","),j=1,Nidx)
      Enddo

      !
      Open (3,file="Clim.txt", Status="Unknown")
      write (3,"(A5,4A11)") "time",trim(Clim_nm(1)),trim(Clim_nm(2)),trim(Clim_nm(12)),trim(Clim_nm(13))
      Do i=1,Ntime
         write (3,"(I5,4)")
      Enddo
      !--------------------------------------------------------

      CONTAINS
      !-----------------------------------------------------------------
      Subroutine Correlation(X,Y,N,Corr,prob,Miss)
        Integer :: N, N1, i, j
        Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t

        N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0
        Do i=1,N
           if (X(i)/=Miss.AND.Y(i)/=Miss) then
              n1=n1+1
              Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i)
              Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i)
              Rxy=Rxy+X(i)*Y(i)
           endif
        Enddo
        if (N1>2) then
           Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1)
           Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb
           tmp=Rxy/Real(N1)
           Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy)
           t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) )
           Prob = T_Dist (t, Real(N1-2))
        else
           Rxy=Miss; prob=Miss
        endif
        Corr=Rxy
        
      End Subroutine Correlation
      !---------------------------------------------------------------------------
      Subroutine Read_ClimIDX (Fname) !,varn)
       Character (Len=*) :: Fname     !,varn
       Real, Allocatable :: Year(:),Month(:), tmp(:,:)
       Integer :: Yr11, Nyr, Nmn, iidx, iyr, i,j, k

       Allocate (Y(Nidx,Ntime))
       Yr11=Yr1 !+1
       Do iidx=1,Nidx
          print*,Ntime,iidx,trim(Clim_nm(iidx))
          Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx)))
          !
          j=0
          Do iyr=Yr11,Yr2
             Do i=1,Nyr
                if (Year(i)==iyr) then
                   Do k=1,Nmon
                      j=j+1; Y(iidx,j)=tmp(i,k)
                   Enddo
                endif
             Enddo
          Enddo
          Deallocate (Year,Month,tmp)
       Enddo
      End Subroutine Read_ClimIDX
      !---------------------------------------------------------------------------------
      Subroutine Moving_Ave (X,N,m,Y)
         Integer :: N,m, i,j,k   ! m must be an odd number
         Real :: X(N), Y(N-m+1), tmp

         Do j=1,N-m+1
            tmp=0
            Do i=j,j+m-1
               tmp=tmp+X(i)
            Enddo
            Y(j)=tmp/Real(m)
         Enddo
      End Subroutine Moving_Ave
      !--------------------------------------------------------------------------------------
      Subroutine Write_to_TXT (Fname,X,N,M,Varname)
        Character (Len=*) :: Fname,Varname
        Integer :: N,M,i,j
        Real :: X(N,M)

        print*,"Write_to_TXT X ",N,M,"--->", trim(Fname)
        Open (3,file=Fname,Status="Unknown")
        Write (3,"(A)") "Varname= "//trim(Varname)
        Do i=1,N
           Write (3,"(I4,50F10.4)") i,(X(i,j),j=1,M)
           print*,i
        Enddo
        Close (3)
        print*,"Writing ",trim(Varname)," ... Done !"
      End Subroutine Write_to_TXT
      !--------------------------------------------------------------------------------------
      Subroutine Anomaly (X,N)
        Integer :: N, it, j, imon, ic
        Real :: X(N), Xtb(Nmon), tmp

        Do imon=1,Nmon
           tmp=0; ic=0
           Do it=1,N
              j=mod(it,Nmon); if (j==0) j=Nmon
              if (j==imon) then; ic=ic+1; tmp=tmp+X(it); endif
           Enddo
           Xtb(imon)=tmp/ic
        Enddo
        Do it=1,N
           j=mod(it,Nmon); if (j==0) j=Nmon
           X(it)=X(it)-Xtb(j)
        Enddo
      End Subroutine Anomaly
      !--------------------------------------------------------------------------------------
      Logical Function Ar1_Empty(X,N)
        Integer :: N, i
        Real :: X(N)
        Logical :: Empty
        Empty=.false.
        Do i=1,N
           if (X(i)==Miss) then
              Empty=.true.; exit
           endif
        Enddo
        Ar1_Empty=Empty
      End Function Ar1_Empty
      !--------------------------------------------------------------------
!      Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90"
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/Ave.f90"
      Include "/home/tanpv/Mylibs/Stdev_1.f90"
      !Include "Sub_MCA.f90"
      !--------------------------------------------------------------------
      END
      Include "/home/tanpv/Mylibs/Student_Prob.f90"