! Maximum Covariance Analysis of Drought Inde:wqx

      Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Yr1=1960, Yr2=2019, Nidx=24
      Integer, Parameter :: Nmode=6, Nreg=6
      !Real, Parameter :: Slon1=70, Slon2=280, Slat1=-20, Slat2=40
      Real, Parameter :: Slon1=110, Slon2=290, Slat1=-30, 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
      !---------------------------------------------------------------- 

      ! 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 
      Real, Allocatable :: Homo(:,:,:),Hete(:,:,:)      ! Homo and Hetero maps (Nlon,Nlat,Nmode)
      Real, Allocatable :: p_Ho(:,:,:),p_He(:,:,:)      ! Pval

      Integer :: Mov_Ave_Step                    ! Buoc truot trong Moving average                
      Real :: Sumx,Corr,prob

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

      Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab,Dir_Out
      !Character (Len=15)  :: Dat_src="Ens" ! "APH" ! "Ens"
      Character (Len=15)  :: Dat_src="CRU" ! "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 = 13           ! Smooth with moving average

      Ntime=(Yr2-Yr1+1)*Nmon

      !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 ! Ntime=",Ntime
 
      Fname="/work/projects/G09/SEACLID/SST/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)
      print*,"Reading SST ... done !"

      Nt=Ntime; Nx=Landpoint; Ny=SeaPoint

      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 = .false.          ! Normalize data (true/false)
      !Calc_mtx=.false.              ! Calculate C-matrix (true/false)
      Calc_mtx=.true.              ! Calculate C-matrix (true/false)
      SVD_MatLab = .True.          ! Use SVD procedure from MatLab (true) or R (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 (SVD_MatLab) then         ! MatLab
            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         ! R
            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 MatLab hoac R

         if (SVD_MatLab) then         ! Use MatLab run
            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                          ! Use R run
            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

      endif

      Dir_Out="/home/tanpv/Progs/2021/NC/"
      !--------------------------------------------------------------------
      ! Write to NC files:
      Fnu=trim(Dir_Out)//"U_Dr_SST.nc"
      Fnv=trim(Dir_Out)//"V_Dr_SST.nc"
      Fnab=trim(Dir_Out)//"AB_Dr_SST.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_SST.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 Get_UVW (Fname,X,Nx,Ny,varn)
      Character (Len=*) :: Fname,varn
      Integer :: N1,N2,ncid, dimid, varid
      Real, Allocatable :: X(:,:)
      Logical :: status

      status = nf90_open(fname, nf90_NoWrite, ncid)
      status = nf90_inq_dimid(ncid,"x",dimid)
      status = nf90_inquire_dimension(ncid,dimid,len=Nx)
      status = nf90_inq_dimid(ncid,"y",dimid)
      status = nf90_inquire_dimension(ncid,dimid,len=Ny)

      Allocate (X(Nx,Ny))
      status = nf90_inq_varid(ncid,trim(varn),varid)
      status = nf90_get_var(ncid,varid,X)
      status =  nf90_close(ncid)
      End Subroutine Get_UVW

      !--------------------------------------------------------------------
      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
         Real :: Sumx
         Integer Nmin 
         Real, Allocatable :: X1d(:), Y1d(:)
         ! Sign inversion of U and V:
         !U = -U; V = -V

         print*,"Nx,Ny,Nmode: ",Nx,Ny,Nmode
         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:
         Allocate (X1d(Nt), Y1d(Nt))
         Do i=1,Nmode
            X1d(1:Nt)=Ak(i,1:Nt); Y1d(1:Nt)=Bk(i,1:Nt)
            Call Correlation (X1d,Y1d,Nt,Corr,Prob, Miss)
            Corre(i)=Corr; Probs(i)=prob
         Enddo
         Deallocate (X1d,Y1d)
         ! Accum (%) of singular values:
         Nmin=Min(Nx,Ny)
         !print*,"(30F6.1)",W(1:Nmin)
         Sumx=Sum(W(1:Nmin)); Accum(:)=0
         Do i=1,Nmode
            Percent(i)=W(i)/Sumx*100
            Accum(i)=Sum(Percent(1:i))
         Enddo
         print*,"Var: ",Percent(:)
         print*,"Cum: ",Accum(:)
         ! 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 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(:)
        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 (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_sst_data (Fname)
        Character (Len=*) :: Fname
        !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:
        !M = Mov_Ave_Step
        !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
                 Call Anomaly (X1d,Ntime)
                 Y(ipoint,1:Nt)=X1d(1:Ntime)
              endif
           Enddo
        Enddo
        Ntime = Nt

      End Subroutine Read_sst_data
      !------------------------------------------------------------------------------
      Subroutine Anomaly (X,Ntime)
      Real :: X(Ntime), Y(Ntime), TB(Nmon)
      Integer :: ic, it, imon, imon1
      Real :: tmp

      if (mod(Ntime,Nmon)/=0) print*,"Something wrong !"
      Do imon=1,Nmon
         tmp=0; ic=0
         Do it=1,Ntime
            imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (imon1==imon) then
               if (X(it)/=Miss) then
                  ic=ic+1; tmp=tmp+X(it)
               endif
            endif
         Enddo
         if (ic>0) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif
      Enddo
      Do it=1,Ntime
         imon=mod(it,Nmon); if (imon==0) imon=Nmon
         Y(it)=X(it)-TB(imon)
      Enddo
      X(:) = Y(:)
      End Subroutine Anomaly

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