! Maximum Covariance Analysis of Drought Inde:wqx

      Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Yr1=1961, Yr2=2005, Nidx=24
      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)=(/ &
                "nino12", "nino3", "nino34", "nino4", "oni", "tni", &   ! 6
                "whwp", "ammsst", "pdo", "NTA_ersst", "dmi",  &   ! 5 Thay CAR_ersst bang pdo, thay tsa bang dmi
                "tna", "amonus", &
                "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 
                "rindo_slpa", "qbo30", "wp", "ao", "nao"/)
                !"rindo_slpa", "tsa", "wp", "ao", "nao"/)

      !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
      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
      !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

      !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)
      print*,"Reading dought index ... done !"
 
      !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="/home/tanpv/Data/ClimIDX/Clim_IDX_2020.nc"
      Call Read_ClimIDX (Fname) !,varn)

      Nt=Ntime; Nx=Landpoint; Ny=Nidx 

      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))
      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))
      Allocate (U3D(Nlon,Nlat,Nmode),V3D(Nlon_sst,Nlat_sst,Nmode))
      Allocate (XX(Nt),YY(Nt))
      
      !--------------------------------------------------------------------

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

      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)
         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)
         endif

         Call Calc_3D_AB

         ! Testing:
         !C = MatMul(MatMul(U,WW),Transpose(V))
         ! 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_Rec.nc",xx1d,yy1d,Nx,Ny,C,"cor", &
         !     "-", "-", "x", "y", Miss)
         !Deallocate (xx1d,yy1d)

         ! Nothing to do more !!!!!!!!!!!
         !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      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)
         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)
         endif

         Call Calc_3D_AB

         ! print test whether U and V are orthogonal:
         print"(6F10.2)", MatMul(Transpose(Vk),Vk)
         ! Testing:
         !C = MatMul(MatMul(U,WW),Transpose(V))
         ! 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_Rec.nc",xx1d,yy1d,Nx,Ny,C,"cor", &
         !     "-", "-", "x", "y", Miss)
         !Deallocate (xx1d,yy1d)

      endif

      ! SVD:
      !print*,"Calling svdcmp ..."
      !Call svdcmp (C, Nx, Ny, U, W, V)
      !print*,"Calling svdcmp ... Done!"

      !Call SVDsrt (U, Nx, Ny, V, W)

      !--------------------------------------------------------------------
      ! Calculating SVD using Lapack:
      !M = Nx; N = Ny; LDA = M; LDU = M; LDVT = N
      !LWMAX = Max(M,N)

      !Allocate (AA(LDA,N), UU(LDU,M), VT(LDVT,N), SS(N), WORK(LWMAX))
      !AA(1:LDA,1:N) = C(1:LDA,1:N)

      !LWORK = -1
      !CALL DGESVD('All','All',M,N,AA,LDA,SS,UU,LDU, VT, LDVT,WORK, LWORK, INFO )
      !LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
      !LWORK = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
      !CALL DGESVD('All','All',M,N,AA,LDA,SS,UU,LDU, VT, LDVT,WORK, LWORK, INFO )

      !IF( INFO.GT.0 ) THEN
      !   WRITE(*,*)'The algorithm computing SVD failed to converge.'
      !   STOP
      !END IF

      !--------------------------------------------------------------------
      !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:
      !Sumx=Sum(W(1:Ny)); 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:
      !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

      
      !--------------------------------------------------------
      ! Write to NC files:
      Fnu="U_Dr_Cli_"//trim(Dat_src)//"_"//trim(varn)//".nc"
      Fnv="V_Dr_Cli_"//trim(Dat_src)//"_"//trim(varn)//".nc"
      Fnab="AB_Dr_Cli_"//trim(Dat_src)//"_"//trim(varn)//".nc"
      Call Write_to_NC (Fnu,Fnv,Fnab)

      !--------------------------------------------------------
      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

         ! 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:
         Sumx=Sum(W(1:Ny)); 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:
         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 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 1962-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
       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 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)) 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_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)) 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)
        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"