! Maximum Covariance Analysis of Drought Inde:wqx

      Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Yr1=1960, Yr2=2019, Nidx=16
      Integer, Parameter :: Nmode=6, Nreg=6
      !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(:)
      Real, Allocatable :: Lon_sst(:), Lat_sst(:), Mask_sst(:,:)
      Real, Allocatable :: Mask(:,:), xlon(:), xlat(:), ylon(:), ylat(:)
      Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j
      Integer :: SeaPoint, Nsst, Nlon_sst, Nlat_sst
      Integer :: Yr_sta
      Real :: Miss_sst
      !---------------------------------------------------------------- 
      Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ &
                "nino3", "nino34", "nino4", "oni", "censo", & ! "tni", &   ! 6
                "rindo_slpa","soi", "reqsoi", "repac_slpa", & !"whwp", "ammsst", "pdo", "NTA_ersst", "dmi",  &   ! 5 Thay CAR_ersst bang pdo, thay tsa bang dmi
                "whwp", "ammsst", "NTA_ersst", "tna", "amonus", &
                "pdo", "dmi"/) ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 
                ! "rindo_slpa", "qbo30", "wp", "ao", "nao"/)
                !"rindo_slpa", "tsa", "wp", "ao", "nao"/)

!Nino3, Nino34, Nino4, ONI, BEST, SOI, REQsoi, REPAC_slp, RINDO_slp,

!WHWP, AMMsst, NTA_SST, TNA, AMONUS,

!PDO, DMI


      !Character (Len=15) :: varn="ori_pdsi" !"spi_12"
      !Character (Len=15) :: varn="spei_12"
      ! 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 
      Integer :: Mov_Ave_Step                    ! Buoc truot trong Moving average                
      Real :: Sumx,Corr,prob

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

      Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab,Dir_Out
      !Character (Len=15)  :: Dat_src="Ens" ! "APH" ! "Ens"
      Character (Len=15)  :: Dat_src="APH" ! "Ens"
      Character (Len=15)  :: varn="spei_12"
      !--------------------------------------------------------------------
      ! LAPACK:
      Integer :: M, N, LDA, LDU, LDVT, LWMAX, INFO, LWORK
      DOUBLE PRECISION, Allocatable :: AA(:,:), UU(:,:), VT(:,:), SS(:), WORK(:)
      EXTERNAL ::        DGESVD
      !--------------------------------------------------------------------

      Mov_Ave_Step = 0           ! Smooth with moving average
      Ntime=(Yr2-Yr1+1)*Nmon

      !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_"//trim(Dat_src)//".nc"
      !Fname="/home/tanpv/Progs/Dr_SEA/Res_APH/BL_"//trim(Dat_src)//".nc"
      !Call Read_drou_data (Fname,varn)

      Fname="/work/projects/G09/SEACLID/GPCC_CRU_SEA/"// &
            "SPEI_SPI_GPCC_CRU_59_19_New.nc"
      Call Read_SPEI (Fname)
      print*,"Reading dought index ... done ! Landpoint, Ntime",Landpoint, Ntime
 
      !Fname="/home/tanpv/Data/SST_globe/sst.mon.mean_COBE.nc" !HadISST_sst.nc"
      !Yr_sta=1891; Miss_sst=1.e+20                            ! Depend on sst source
      !Call Read_sst_data (Fname,varn)
      !print*,"Reading SST ... done !"

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

      Nt=Ntime; Nx=Landpoint; Ny=Nidx 
      print*,"Nx,Ny,Nt: ",Nx,Ny,Nt

      Allocate (Xc(Nx,Nt),Yc(Ny,Nt),Xbar(Nx),Ybar(Ny), Sx(Nx), Sy(Ny))
      !Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nt,Nx), B(Nt,Ny))
      Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nx,Nt), B(Ny,Nt))
      Allocate (WW(Ny,Ny))
!      print*,">>>>>>>>>>>"
      Allocate (Uk(Nx,Nmode), Vk(Ny,Nmode), WWk(Nmode,Nmode), Wk(Nmode))
      !Allocate (Ak(Nt,Nmode), Bk(Nt,Nmode), Corre(Nmode),Probs(Nmode))
      Allocate (Ak(Nmode,Nt), Bk(Nmode,Nt), Corre(Nmode),Probs(Nmode))
      Allocate (Accum(Nmode),Percent(Nmode))
!            print*,">>>>>>>>>>>"

      Allocate (U3D(Nlon,Nlat,Nmode)) !,V3D(Nlon_sst,Nlat_sst,Nmode))
      print*,">>>>>>>>>>>"
      Allocate (XX(Nt),YY(Nt))
      
      !--------------------------------------------------------------------

      print*,"Centering data ..."
      Normalize = .true.
      SDV_MatLab = .true.
      Calc_mtx=.true.
      !Calc_mtx=.false.
      !--------------------------------------------------------------------

      Call Centered_Data (Normalize)            ! Su dung de tinh Ak, Bk

      if (Calc_mtx) then            ! Chi tinh Ma tran Covariance

         ! Covariance matrix:
         print*,"Calculating Covariance matrix..."
         C = MatMul(Xc,Transpose(Yc))/Real(Nt-1)

         ! Write C to nc file:
         Allocate (xx1d(Nx), yy1d(Ny))
         Do i=1,Nx; xx1d(i)=i; Enddo
         Do i=1,Ny; yy1d(i)=i; Enddo
         Call wr_ncdf_2d ("C_Matrix.nc",xx1d,yy1d,Nx,Ny,C,"cor", &
              "-", "-", "x", "y", Miss)
         Deallocate (xx1d,yy1d)

         ! Calculating SVD:
         if (SDV_MatLab) then
            !isys = system ("./run_mlab.sh ")
            !Fnu="/home/tanpv/Progs/Dr_SEA/MCA/NC/U.nc"
            !Fnv="/home/tanpv/Progs/Dr_SEA/MCA/NC/V.nc"
            !Fnab="/home/tanpv/Progs/Dr_SEA/MCA/NC/S.nc"
            !Call Read_UVW_MatLab (Fnu,Fnv,Fnab)

            isys = system ("./run_mlab.sh ")
            Fnu="/home/tanpv/Progs/2021/NC/U.nc"
            Fnv="/home/tanpv/Progs/2021/NC/V.nc"
            Fnab="/home/tanpv/Progs/2021/NC/S.nc"
            Call Read_UVW_MatLab (Fnu,Fnv,Fnab)
         else
            !isys = system ("./run_r.sh")
            !Fnu="/home/tanpv/Progs/Dr_SEA/MCA/NC/u.nc"
            !Fnv="/home/tanpv/Progs/Dr_SEA/MCA/NC/v.nc"
            !Fnab="/home/tanpv/Progs/Dr_SEA/MCA/NC/w.nc"
            !Call Read_UVW (Fnu,Fnv,Fnab)

            isys = system ("./run_r.sh ")
            Fnu="/home/tanpv/Progs/2021/NC/u.nc"
            Fnv="/home/tanpv/Progs/2021/NC/v.nc"
            Fnab="/home/tanpv/Progs/2021/NC/w.nc"
            Call Read_UVW (Fnu,Fnv,Fnab)
         endif

         Call Calc_3D_AB

         !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      else                          ! Xu ky KQ sau khi thuc hien SVD bang R

         if (SDV_MatLab) then
            !Fnu="/home/tanpv/Progs/Dr_SEA/MCA/NC/U.nc"
            !Fnv="/home/tanpv/Progs/Dr_SEA/MCA/NC/V.nc"
            !Fnab="/home/tanpv/Progs/Dr_SEA/MCA/NC/S.nc"
            !Call Read_UVW_MatLab (Fnu,Fnv,Fnab)

            Fnu="/home/tanpv/Progs/2021/NC/U.nc"
            Fnv="/home/tanpv/Progs/2021/NC/V.nc"
            Fnab="/home/tanpv/Progs/2021/NC/S.nc"
            Call Read_UVW_MatLab (Fnu,Fnv,Fnab)
         else
            !Fnu="/home/tanpv/Progs/Dr_SEA/MCA/NC/u.nc"
            !Fnv="/home/tanpv/Progs/Dr_SEA/MCA/NC/v.nc"
            !Fnab="/home/tanpv/Progs/Dr_SEA/MCA/NC/w.nc"
            !Call Read_UVW (Fnu,Fnv,Fnab)

            Fnu="/home/tanpv/Progs/2021/NC/u.nc"
            Fnv="/home/tanpv/Progs/2021/NC/v.nc"
            Fnab="/home/tanpv/Progs/2021/NC/w.nc"
            Call Read_UVW (Fnu,Fnv,Fnab)
         endif

         Call Calc_3D_AB

         ! print test whether U and V are orthogonal:
         print"(6F10.2)", MatMul(Transpose(Vk),Vk)

      endif

      Dir_Out="/home/tanpv/Progs/2021/NC/"
      !--------------------------------------------------------
      ! Write to NC files:
      Fnu=trim(Dir_Out)//"U_Dr_Cli_Revise.nc"
      Fnv=trim(Dir_Out)//"V_Dr_Cli_Revise.nc"
      Fnab=trim(Dir_Out)//"AB_Dr_Cli_Revise.nc"
      Call Write_to_NC (Fnu,Fnv,Fnab)

      ! Calculate and write Homogenerous and Heterogenerous correlation maps: 
      Call Homo_And_Hetero_Map
      Fnu =trim(Dir_Out)//"Ho_He_Cli_Revise.nc"
      Call Write_Homo_Hete_Maps_to_NC (Fnu)
      !--------------------------------------------------------
      print*,"Wk", Wk
      print*,"Corr=", Corre
      print*,"Prob=",Probs
      print*,"Perc=",Percent
      print*,"Accum",Accum

      CONTAINS
      !--------------------------------------------------------------------
      Subroutine Read_UVW_MatLab (Fu,Fv,Fw)
        Character (Len=*) :: Fu,Fv,Fw
        Real, Allocatable :: UU(:,:)
        Real, Allocatable :: xlo(:), xla(:)
        Integer :: nlon,nlat,i,j

        Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S")
        print*, nlon,nlat,nx,ny

        !Do i=1,Ny; W(i) = SQRT(UU(i,i)); Enddo
        Do i=1,Ny; W(i) = UU(i,i); Enddo

        Deallocate (xlo,xla,UU)

        Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"U")
        print*, nlon,nlat,nx,ny
        Do i=1,Nx
           U(i,1:Ny) = UU(i,1:Ny) !*W(i)
        Enddo
        Deallocate (xlo,xla,UU)
        Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"V")
        print*, nlon,nlat,nx,ny

        Do i=1,Ny
           V(i,1:Ny) = UU(i,1:Ny) !*W(i)
        Enddo
        Deallocate (xlo,xla,UU)
        !Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S")
        !print*, nlon,nlat,nx,ny

        !Do i=1,Ny; W(i) = UU(i,i); Enddo
        !Deallocate (xlo,xla,UU)
      End Subroutine Read_UVW_MatLab
      !--------------------------------------------------------------------
      Subroutine Read_UVW (Fu,Fv,Fw)
        Character (Len=*) :: Fu,Fv,Fw
        Real, Allocatable :: UU(:,:)
        Real, Allocatable :: xlo(:), xla(:)
        Integer :: nlon,nlat,i,j

        Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"var")
        print*, nlon,nlat,nx,ny
        Do i=1,Nx
           U(i,1:Ny) = UU(i,1:Ny)
        Enddo
        Deallocate (xlo,xla,UU)
        Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"var")
        print*, nlon,nlat,nx,ny

        Do i=1,Ny
           V(i,1:Ny) = UU(1:Ny,i)
        Enddo
        Deallocate (xlo,xla,UU)
        Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"var")
        print*, nlon,nlat,nx,ny

        W(1:Ny) = UU(1:Ny,1)
        Deallocate (xlo,xla,UU)
      End Subroutine Read_UVW
      !--------------------------------------------------------------------
      Subroutine Calc_3D_AB
         Integer :: Nmin
         Real :: Sumx

         ! Sign inversion of U and V:
         !U = -U; V = -V

         Wk(1:Nmode) = W(1:Nmode)                  ! Singular values
         Do i=1,Nmode
            Uk(:,i)=U(:,i); Vk(:,i)=V(:,i)
         Enddo

         ! Singular matrix:
         WW=0;  Do j=1,Ny; WW(j,j)=W(j); Enddo      ! Full 
         WWk=0; Do j=1,Nmode; WWk(j,j)=W(j); Enddo  ! Nmode


         ! Timeseries of Nmode
         Ak = MatMul(Transpose(Uk),Xc)
         Bk = MatMul(Transpose(Vk),Yc)

         ! Correlations:
         Do i=1,Nmode
            XX(1:Nt)=Ak(i,1:Nt); YY(1:Nt)=Bk(i,1:Nt)
            Call Correlation (XX,YY,Nt,Corr,Prob, Miss)
            Corre(i)=Corr; Probs(i)=prob
         Enddo
      
         ! Accum (%) of singular values:
         Nmin=Min(Nx,Ny)
         Sumx=Sum(W(1:Nmin)); Accum(:)=0
         Do i=1,Nmode
            Percent(i)=W(i)/Sumx*100
            Accum(i)=Sum(Percent(1:i))
         Enddo

         ! Convert Uk to U3D:
         U3D=Miss
         Do i=1,Nmode
            j=0
            Do ilat=1,Nlat
               Do ilon=1,Nlon
                  if (Mask(ilon,ilat)==1) then
                     j=j+1
                     U3D(ilon,ilat,i) = Uk(j,i)
                  endif
               Enddo
            Enddo
         Enddo

         ! Convert Vk to V3D:  >>>>> Do not need for ClimIDX
         !V3D=Miss
         !Do i=1,Nmode
         !   j=0
         !   Do ilat=1,Nlat_sst
         !      Do ilon=1,Nlon_sst
         !         if (Mask_sst(ilon,ilat)==1) then
         !            j=j+1
         !            V3D(ilon,ilat,i) = Vk(j,i)
         !         endif
         !      Enddo
         !   Enddo
         !Enddo
  
      End Subroutine Calc_3D_AB
      !--------------------------------------------------------------------
      Subroutine Homo_And_Hetero_Map
        ! Calculate Homogenerous correlation map = r[X(t), a_k (t)]
        ! and     Heterogenerous correlation map = r[X(t), b_k (t)]

        Integer :: ilon,ilat,it,i,j
        Real :: ctmp, ptmp
        Real, Allocatable :: Ho(:,:),He(:,:),pHo(:,:),pHe(:,:),X1d(:), A1d(:), B1d(:)

        Allocate (X1d(Nt),A1d(Nt),B1d(Nt),Ho(Nx,Nmode),He(Nx,Nmode))
        Allocate (pHo(Nx,Nmode),pHe(Nx,Nmode))

        Do i=1,Nmode
           A1d(1:Nt)=Ak(i,1:Nt); B1d(1:Nt)=Bk(i,1:Nt)
           Do j=1,Nx
              X1d(1:Nt)=X(j,1:Nt)
              Call Correlation(X1d,A1d,Nt,ctmp,ptmp,Miss)
              Ho(j,i)=ctmp; pHo(j,i)=ptmp
              Call Correlation(X1d,B1d,Nt,ctmp,ptmp,Miss)
              He(j,i)=ctmp; pHe(j,i)=ptmp
           Enddo
        Enddo
        ! Convert to 3D:
        Allocate (Homo(Nlon,Nlat,Nmode),p_Ho(Nlon,Nlat,Nmode))
        Allocate (Hete(Nlon,Nlat,Nmode),p_He(Nlon,Nlat,Nmode))
         Homo=Miss; Hete=Miss; p_Ho=Miss; p_He=Miss
         Do i=1,Nmode
            j=0
            Do ilat=1,Nlat
               Do ilon=1,Nlon
                  if (Mask(ilon,ilat)==1) then
                     j=j+1
                     Homo(ilon,ilat,i)=Ho(j,i); p_Ho(ilon,ilat,i)=pHo(j,i)
                     Hete(ilon,ilat,i)=He(j,i); p_He(ilon,ilat,i)=pHe(j,i)
                  endif
               Enddo
            Enddo
         Enddo

      End Subroutine Homo_And_Hetero_Map
      !--------------------------------------------------------------------
      Subroutine Write_Homo_Hete_Maps_to_NC (Fnu)
        Character (Len=*) :: Fnu
        Real, Allocatable :: xmod(:)
        Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm
        Integer :: i,j,k

        Allocate (xmod(Nmode))
        Do i=1,Nmode; xmod(i)=i; Enddo
        t_unit="months since 1962-01-01 "
        var_ln="Homogenerous correlation map r[Dr(t), a_k (t)]"
        Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,Homo,"ho", &
             t_unit, var_ln, "-", "lon", "lat", "mod", Miss)
        var_ln="p-values of Homo"
        Call wr_ncdf_3d_apd (Fnu,p_Ho,"pho",nlon,nlat,Nmode, var_ln, &
             "-", "lon", "lat", "mod", Miss)
        var_ln="Heterogenerous correlation map r[Dr(t), b_k (t)]"
        Call wr_ncdf_3d_apd (Fnu,Hete,"he",nlon,nlat,Nmode, var_ln, &
             "-", "lon", "lat", "mod", Miss)
        var_ln="p-values of Hete"
        Call wr_ncdf_3d_apd (Fnu,p_He,"phe",nlon,nlat,Nmode, var_ln, &
             "-", "lon", "lat", "mod", Miss)
      End Subroutine Write_Homo_Hete_Maps_to_NC
      !--------------------------------------------------------------------
      Subroutine Write_to_NC (Fnu,Fnv,Fnab)
        Character (Len=*) :: Fnu, Fnv, Fnab  ! Filename of U, V and A,B
        Real, Allocatable :: xmod(:), X2D(:,:), xtime(:), climidx(:)
        Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm
        Integer :: i,j,k

        Allocate (xmod(Nmode),X2D(Nmode,Ny),xtime(Nt),climidx(Nidx))

        Do i=1,Nmode; xmod(i)=i; Enddo
        Do i=1,Nt; xtime(i)=i; Enddo
        Do i=1,Nidx; climidx(i)=i; Enddo

        ! U matrix:
        t_unit="months since 1961-01-01 "
        var_ln="MCA for drought index and SST - U-matrix"
        Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,U3D,"u", &
             t_unit, var_ln, "-", "lon", "lat", "mod", Miss)

        ! V matrix:
        !var_ln="MCA for drought index and SST - V-matrix"
        !Call wr_ncdf_3d (Fnv,Lon_sst,Lat_sst,xmod,nlon_sst,nlat_sst,Nmode,V3D,"v", &
        !     t_unit, var_ln, "-", "lon", "lat", "mod", Miss)

        ! V matrix:
        Do i=1,Nmode
           X2D(i,:)=Vk(:,i)
        Enddo
        var_ln="MCA for climate indices"
        Call wr_ncdf_2d (Fnv,xmod,climidx,Nmode,Nidx, X2D, "v", var_ln, &
                         "-", "mod", "idx", Miss)

        !var_ln="Percentage of contribution of singular values"
        !Call wr_ncdf_1d_apd (Fnv,Percent,"per", Nmode, var_ln,"-","mod", Miss)
        !var_ln="Accumulated Percent. of contribution of singular values"
        !Call wr_ncdf_1d_apd (Fnv,Accum,"cum", Nmode, var_ln,"-","mod", Miss)

        ! A & B matrix
        var_ln="Timeseries of MCA of drought index (Ak)"
        Call wr_ncdf_2d (Fnab,xmod,xtime,Nmode,Nt,Ak, "a", var_ln, &
                         "-", "mod", "mon", Miss)
        var_ln="Timeseries of MCA of SST (Bk)"
        Call wr_ncdf_2d_apd (Fnab,Bk,"b",Nmode,Nt, var_ln, "-", &
                         "mod", "mon", Miss)
        var_ln="Correlation coefficients between A & B"
        Call wr_ncdf_1d_apd (Fnab,Corre,"cor", Nmode, var_ln,"-","mod", Miss)
        var_ln="p-values of Correlation coefficients between A & B"
        Call wr_ncdf_1d_apd (Fnab,probs,"prob", Nmode, var_ln,"-","mod", Miss)

        var_ln="Percentage of contribution of singular values"
        Call wr_ncdf_1d_apd (Fnab,Percent,"per", Nmode, var_ln,"-","mod", Miss)
        var_ln="Accumulated Percent. of contribution of singular values"
        Call wr_ncdf_1d_apd (Fnab,Accum,"cum", Nmode, var_ln,"-","mod", Miss)

      End Subroutine Write_to_NC
      !--------------------------------------------------------------------
      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 Centered_Data (Normalize)
        Integer :: i,j
        Logical :: Normalize
        Real, Allocatable :: Std(:), X1d(:)
       
        ! Centering:
        Allocate (Std(Nx), X1d(Nt))
        Do j=1,Nx
           X1d(1:Nt)=X(j,1:Nt)
           Xbar(j)=Ave(X1d,Nt,Miss)
           Std (j)=Stdev_1(X1d,Nt,Miss)
           Do i=1,Nt
              Xc(j,i)=X(j,i)-Xbar(j)
              if (Normalize) Xc(j,i)=Xc(j,i)/Std(j) !/Stdev_1(X(j,1:Nt),Nt,Miss)
           Enddo
        Enddo
        Deallocate (Std); Allocate (Std(Ny))
        Do j=1,Ny
           X1d(1:Nt)=Y(j,1:Nt)
           Ybar(j)=Ave(X1d,Nt,Miss)
           Std (j)=Stdev_1(X1d,Nt,Miss)
           Do i=1,Nt
              Yc(j,i)=Y(j,i)-Ybar(j)
              if (Normalize) Yc(j,i)=Yc(j,i)/Std(j) ! Stdev_1(Y(j,1:Nt),Nt,Miss)
           Enddo
        Enddo
      End Subroutine Centered_Data
      !--------------------------------------------------------------------------------------
      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  Lay bang chuoi SPEI sau khi da tru nam dau tien
       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
          if (j/=Ntime) then; print*,"Check ClimIDX: ",j,Ntime; endif
          Deallocate (Year,Month,tmp)
       Enddo
      End Subroutine Read_ClimIDX
      !--------------------------------------------------------------------
      Subroutine Read_sst_data (Fname,varn)
        Character (Len=*) :: Fname,varn
        !Integer, Parameter :: Yr_sta=1891   ! Nam bat dau chuoi SL SST
        Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:)
        Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlon, nlat, imon, imon1
        Integer :: M, Md2, it1, it2, ilon1, ilon2,ilat1,ilat2
        Character (Len=255) :: Fn

        ! Lon_sst(:), Lat_sst(:), Mask_sst(:,:), SeaPoint, Nsst, Nlon_sst, Nlat_sst

        Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,XX,"sst") !trim(varn))
        Where (IsNaN(XX)); XX=Miss; Endwhere
        Where (XX==Miss_sst); XX=Miss; Endwhere
        i=1; Do while (tlon(i)<Slon1); i=i+1; Enddo; ilon1=i
        i=1; Do while (tlon(i)<Slon2); i=i+1; Enddo; ilon2=i
        j=1; Do while (tlat(j)>Slat1); j=j+1; Enddo; ilat1=j
        j=1; Do while (tlat(j)>Slat2); j=j+1; Enddo; ilat2=j
        it1=(Yr1-Yr_sta)*Nmon+1 + Nmon         ! Bo nam dau tien
        it2=it1+Ntime-1

        ! Lay them M/2 phan tu dau va M/2 phan tu cuoi cho Moving Average:
        Md2=M/2
        it1=it1-Md2; it2=it2+Md2
        Nt = Ntime                ! Luu gia tri Ntime
        Ntime = it2-it1+1         ! Kich thuoc chuoi moi dung trong Moving Average

        Nlon_sst=ilon2-ilon1+1; 
        Nlat_sst=ilat1-ilat2+1       ! Chu y: ilat1>ilat2 vi Latitude giam tu North to South

        !print*,Yr_sta+(it1-1)/Nmon, Yr_sta+it2/Nmon, ilon1,ilon2,ilat1, ilat2
        !print*,tlon(ilon1),tlon(ilon2),tlat(ilat1),tlat(ilat2); pause

        Allocate (YY(Nlon_sst,Nlat_sst,Ntime),Lon_sst(Nlon_sst), Lat_sst(Nlat_sst))
        Allocate (Mask_sst(Nlon_sst,Nlat_sst), X1d(Ntime))
        Mask_sst(:,:)=Miss
        X1d(1:Ntime) = xtime(it1:it2)
        i=0
        Do ilon=1,Nlon
           if (ilon>=ilon1.AND.ilon<=ilon2) then
              i=i+1; Lon_sst(i)=tlon(ilon)
              j=Nlat_sst+1
              Do ilat=1,Nlat
                 if (ilat>=ilat2.AND.ilat<=ilat1) then
                    j=j-1; Lat_sst(j)=tlat(ilat)
                    YY(i,j,1:Ntime) = XX(ilon,ilat,it1:it2)
                 endif
              Enddo
           endif
        Enddo
        ! Processing for MCA:
        SeaPoint=0
        Do i=1,Nlon_sst
           Do j=1,Nlat_sst
              X1d(1:Ntime) = YY(i,j,1:Ntime)
              if (.not.Ar1_Empty(X1d,Ntime,Miss)) then 
                 SeaPoint=SeaPoint+1; Mask_sst(i,j)=1
              endif
           Enddo
        Enddo

        ! Moving Average:
        M = Mov_Ave_Step               ! Buoc truot = 13 thang
        if (M>0) then
           Nt=Ntime-M+1
        else
           Nt=Ntime
        endif

        Allocate (Y(SeaPoint,Nt),Y1d(Nt)) 
        ipoint=0
        Do ilat=1, Nlat_sst
           Do ilon=1, Nlon_sst
              if (Mask_sst(ilon,ilat)/=Miss) then
                 ipoint=ipoint+1
                 X1d(1:Ntime)=YY(ilon,ilat,1:Ntime)
                 if (M>0) then
                    Call Moving_Ave (X1d,Ntime,M,Y1d)
                    Y(ipoint,1:Nt)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime)
                 else
                    Y(ipoint,1:Nt)=X1d(1:Ntime)   ! Nt=Ntime
                 endif
              endif
           Enddo
        Enddo
        Ntime = Nt

      End Subroutine Read_sst_data
      !---------------------------------------------------------------------------------
      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 Read_SPEI (Fname)
      Character (Len=*) :: Fname
      Integer :: ilon,ilat,it,it1, Nt, Nt1, ix
      Real, Allocatable :: spe3d(:,:,:), X3d(:,:,:), Y3d(:,:,:), xtime(:), x1d(:)


      !Lon_spe(:), Lat_spe(:)
      !Call rd_ncdf_3d (Fname,lon_spe,lat_spe,xtime,nlon_spe,nlat_spe,nt,X3d,"spei12")
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,Nt,X3d,"spei12")
      !Deallocate (lon,lat,xtime)
      !Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,Mask,"land")
      Where (IsNaN(X3d)) X3d=Miss
      it1=Nmon+1; Nt1=Nt-Nmon

      if (Nt1/=Ntime) then; print*,"Check: Nt,Ntime:",Nt1,Ntime; stop; endif

      Allocate (spe3d(nlon,nlat,Nt1),mask(nlon,nlat))
      Allocate (x1d(Nt1))
      spe=Miss; mask=Miss; Landpoint=0
      Do ilon=1,nlon
         Do ilat=1,nlat
            x1d(1:Nt1)=X3d(ilon,ilat,it1:Nt)
            if (.not.Ar1_Empty (X1d,Nt1,Miss)) then
               mask(ilon,ilat)=1
               Landpoint=Landpoint+1
               spe3d(ilon,ilat,1:Nt1)=X3d(ilon,ilat,it1:Nt)
            endif
         Enddo
      Enddo
      ! Reshape to 2D arr: Ntime=Nt1
      Allocate (X(Landpoint,Ntime))
      !Allocate (spe2d(LandSpe,Ntime))
      ix=0
      Do ilat=1,Nlat
         Do ilon=1,Nlon
            if (mask(ilon,ilat)/=Miss) then
               ix=ix+1; X(ix,1:Ntime)=spe3d(ilon,ilat,1:Ntime)
            endif
         Enddo
      Enddo
     End Subroutine Read_SPEI
      !--------------------------------------------------------------------
      Subroutine Read_drou_data (Fname,varn)
        Character (Len=*) :: Fname,varn
        Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:)
        Integer :: M, ist, ilon, ilat, ipoint, nt, it, nlo, nla, imon, imon1

        ! Read in drought index:
        print*,"Reading Drought..."
        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,nt,XX,trim(varn))

        Where (IsNaN(XX)); XX=Miss; Endwhere
        
        Ntime=Nt-Nmon; it=Nmon+1
        Allocate (PDSI(Nlon,Nlat,Ntime)) !,Dr_IDX(Nlon,Nlat,Ntime))

        PDSI(:,:,1:Ntime) = XX(:,:,it:Nt)     ! Bo Nmon thang dau tien

        Nst=Nlon*Nlat
        Allocate (Mask(Nlon,Nlat))
        Allocate (lon(Nst), lat(Nst), tmp(Ntime))

        ist=0; Landpoint=0; Mask(:,:)=Miss
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              ist=ist+1;
              Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat)
              tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (.not.Ar1_Empty(tmp,Ntime,Miss)) then !.AND.Reg2D(ilon,ilat)/=Miss) then
                 Landpoint=Landpoint+1; Mask(ilon,ilat)=1
              endif
           Enddo
        Enddo

        ! Moving Average:
        M = Mov_Ave_Step               ! Buoc truot = 13 thang
        if (M>0) then
           Nt=Ntime-M+1
        else
           Nt=Ntime
        endif

        Allocate (X(Landpoint,Nt),Dr_IDX(Nlon,Nlat,Nt),Y1d(Nt),X1d(Ntime)) !,Reg(Landpoint))
        print*,Ntime,Nt,Nlon,Nlat
        ipoint=0
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              if (Mask(ilon,ilat)/=Miss) then
                 ipoint=ipoint+1
                 X1d(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
                 !X(ipoint,1:Ntime)=PDSI(ilon,ilat,1:Ntime)
                 if (M>0) then
                    Call Moving_Ave (X1d,Ntime,M,Y1d)
                    X(ipoint,1:Nt)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime)
                    Dr_IDX(ilon,ilat,1:Nt)=Y1d(1:Nt)
                 else
                    X(ipoint,1:Nt)=X1d(1:Ntime)   ! Nt=Ntime
                    Dr_IDX(ilon,ilat,1:Nt)=X1d(1:Ntime)
                 endif
                 !print*,ipoint
                 !Dr_IDX(ilon,ilat,1:Ntime)=PDSI(ilon,ilat,1:Ntime)  ! For further use
              endif
           Enddo
        Enddo
        Ntime=Nt
        !print*,X; pause
        Deallocate (PDSI,XX,xtime)
      End Subroutine Read_drou_data
      !--------------------------------------------------------------------------------------
      Logical Function Ar1_Empty(X,N,Miss)
        Integer :: N, i
        Real :: X(N), Miss
        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"