! Doc KQ tinh U,V,W tu R va tinh lai Uk,Vk,A,B 

      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=80, Slon2=270, Slat1=-10, Slat2=30

      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) :: 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                 
      Real :: Sumx,Corr,prob

      Real, Allocatable :: xx1d(:), yy1d(:)
      !--------------------------------------------------------------------

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

      !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="/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)


      Nt=Ntime; Nx=Landpoint; Ny=SeaPoint

      Allocate (Xc(Nx,Nt),Yc(Ny,Nt),Xbar(Nx),Ybar(Ny))
      Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nx,Nt), B(Ny,Nt))
      !Allocate (C(Nx,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))
      
      !--------------------------------------------------------------------
      Call Centered_Data

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

      !print"(30F6.2)",C, "C"; pause
      !print*,"Calculating Covariance matrix... Done!"

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

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

      !--------------------------------------------------------------------
      !U = -U; V = -V

      Fnu="/work/users/Exchange/new_program/data_u.nc"
      Fnv="/work/users/Exchange/new_program/data_v.nc"
      Fnab="/work/users/Exchange/new_program/data_d.nc"
      !Fnab="/work/users/Exchange/new_test/old_d.nc"
      !Fnu="/work/users/Exchange/old_u.txt"
      !Fnv="/work/users/Exchange/old_v.txt"
      !Fnab="/work/users/Exchange/old_d.txt"
      Call Read_UVW (Fnu,Fnv,Fnab)

      print*,"Reading U,V,W... dome"
      !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
      print*,"Calculating Correlations... dome"
      
      ! 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_"//trim(Dat_src)//"_"//trim(varn)//".nc"
      Fnv="V_SST_"//trim(Dat_src)//"_"//trim(varn)//".nc"
      Fnab="AB_"//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 (Fu,Fv,Fw)
        Character (Len=*) :: Fu,Fv,Fw
        Real, Allocatable :: UU(:,:)
        Real, Allocatable :: xlo(:), xla(:)
        Integer :: nlon,nlat,i,j

        ! U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nx,Nt), B(Ny,Nt)

      Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"var") 
      !Open (1,file=Fu)
      Do i=1,Nx
         U(i,1:Ny) = UU(1:Ny,i)
         !Read (1,*) (U(i,j),j=1,Ny)
      Enddo
      !Close (1)
      Deallocate (xlo,xla,UU)
      Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"var")
      !Open (1,file=Fv)
      Do i=1,Ny
         V(i,1:Ny) = UU(1:Ny,i)
         !Read (1,*) (V(i,j),j=1,Ny)
      Enddo
      !Close(1)
      Deallocate (xlo,xla,UU)
      Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"var")
      W(1:Ny) = UU(1,1:Ny)
      !print*,W(1:Ny)
      Deallocate (xlo,xla,UU)
      !Open (1,file=Fw)
      !Do i=1,Ny
      !   Read (1,*) W(i)
      !Enddo

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

      Ak = MatMul(Transpose(Uk),Xc)
      Bk = MatMul(Transpose(Vk),Yc)


      End Subroutine Read_UVW
      !--------------------------------------------------------------------
      Subroutine Write_to_NC (Fnu,Fnv,Fnab)
        Character (Len=*) :: Fnu, Fnv, Fnab  ! Filename of U, V and A,B
        Real, Allocatable :: xmod(:), X2D(:,:), xtime(:)
        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))

        Do i=1,Nmode; xmod(i)=i; Enddo
        Do i=1,Nt; xtime(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)

        !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
        Integer :: i,j

        ! Centering:
        Do j=1,Nx
           Xbar(j)=Ave(X(j,1:Nt),Nt,Miss)
           Do i=1,Nt
              Xc(j,i)=X(j,i)-Xbar(j)
           Enddo
        Enddo
        Do j=1,Ny
           Ybar(j)=Ave(Y(j,1:Nt),Nt,Miss)
           Do i=1,Nt
              Yc(j,i)=Y(j,i)-Ybar(j)
           Enddo
        Enddo
      End Subroutine Centered_Data
      !--------------------------------------------------------------------
      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 :: 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
        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

        Allocate (Y(SeaPoint,Ntime)) 
        ipoint=0
        Do ilat=1, Nlat_sst
           Do ilon=1, Nlon_sst
              if (Mask_sst(ilon,ilat)/=Miss) then
                 ipoint=ipoint+1
                 Y(ipoint,1:Ntime)=YY(ilon,ilat,1:Ntime)
              endif
           Enddo
        Enddo


        !Fn="/home/tanpv/Data/SST_globe/SST_MCA_SEA.nc"
        !Call wr_ncdf_3d (Fn,Lon_sst,Lat_sst,X1d,Nlon_sst,Nlat_sst,Ntime,YY,"sst", &
        !"days since 1961-1-1", "SST over SEA", "C", "lon", "lat", "time", Miss)

      End Subroutine Read_sst_data
      !--------------------------------------------------------------------
      Subroutine Read_drou_data (Fname,varn)
        Character (Len=*) :: Fname,varn
        Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:)
        Integer :: 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
        Allocate (X(Landpoint,Ntime)) !,Reg(Landpoint))
        ipoint=0
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              !tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (Mask(ilon,ilat)/=Miss) then
              !if (.not.Ar1_Empty(tmp,Ntime)) then
                 ipoint=ipoint+1
                 !Reg(ipoint)=Reg2D(ilon,ilat)                       ! 1D Reg
                 X(ipoint,1:Ntime)=PDSI(ilon,ilat,1:Ntime)
                 Dr_IDX(ilon,ilat,1:Ntime)=PDSI(ilon,ilat,1:Ntime)  ! For further use
              endif
           Enddo
        Enddo
        !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 "Sub_MCA.f90"
      !--------------------------------------------------------------------
      END
      Include "/home/tanpv/Mylibs/Student_Prob.f90"