! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10)
      Use NetCDF
      Real, Parameter :: Miss=-99.

      Real*8, Allocatable :: Xbar(:), Std(:), Rx(:), D(:), V(:)
      Real*8, Allocatable :: X(:,:), Y(:,:), YT(:,:)
      ! X(nt,nx): SL ban dau, Y(nt,nx): SL quy tam hoac chuan hoa, YT-Chuyen vi cua Y 

      Real*8, Allocatable :: Fx(:,:), Ax(:,:)   ! X = Fx*A'x
      ! Fx, Ax la cac mang tinh trung gian

      Real*8, Allocatable :: EOF(:,:), PCs(:,:)
      ! EOF(nt,K), PCs(nx,K): Cac ham thoi gian va ham khong gian

      Real*8, Allocatable :: EVtor(:,:), Lamda(:,:)
      Real, Allocatable :: EValues(:), ECumulate(:)
      !---------------------------------------------------------------------------------
      !For PCA's calculations:
      Real*8, Allocatable :: RR(:,:), A(:), R(:)

      Integer :: Nt, Nx  ! Kich thuoc thuc cua mang X
      Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nth0, Nthre, K
      Real :: Threshold
      Character (Len=250) :: Fname, Fname1
      Real :: Factor

      Namelist /Model_Info/ Fname, switch, standard, PC_Stard, Nmode

      Integer :: i,j,L
      ! Mo ta thong tin SL dau vao:
      !---------------------------------------------------------------------------------
      Real, allocatable :: xlon(:), xlat(:), xtime(:)
      Integer :: Nlon, Nlat, Ntime
      Logical :: Status

      Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:)
      Real, allocatable ::  Mask(:,:), Mask_all(:,:)
      Real, allocatable :: Dr_IDX(:,:,:)
      Integer :: Nst, iidx
      Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12
      Character (Len=5), Parameter :: IDX_nm(NDr_idx)= &
                                      (/"pdsi","spi12","spe12"/)
                                   !(/"PDSI","SPI12","SPE12"/)
      !---------------------------------------------------------------------------------
      !For Writing EOFs, PCs:
      Real, allocatable :: PCA(:,:,:), xmode(:), Xtmp(:), EOF1(:,:)
      !---------------------------------------------------------------------------------

      Open (1,file="Namelist_PCA_New.txt",Status="old")
      Read (1, Model_Info)
      Close (1)
      !---------------------------------------------------------------------------------
      ! Reading data:
      Call Read_Drough_IDX (Fname)
      Nt = Ntime; Nx = Nst

      ! Danh cho qua trinh tinh EOFs, PCs:
      !---------------------------------------------------------------------------------
      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Factor = Sqrt(Real(Nt))


      Allocate (X(Nt,Nx), Y(Nt,Nx), YT(Nx,Nt))

      Allocate (Xbar(Mcount), Std(Mcount), Rx(Mcount*Mcount))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (EValues(Mcount), V(Mcount*Mcount))              ! Tri rieng va Vecto rieng (1D)
      Allocate (Fx(Ncount,Mcount), Ax(Mcount,Mcount))           ! EOF va PCA
      Allocate (D(Mcount), EVtor(Mcount,Mcount))                ! Working arr va matran Vecto rieng
      Allocate (Lamda(Mcount,Mcount))
      Allocate (EOF(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
      !---------------------------------------------------------------------------------
      !For PCA's calculations:
      Allocate (RR(Mcount,Mcount), A(Mcount*Mcount), R(Mcount*Mcount))
      !---------------------------------------------------------------------------------
      !For Writing EOFs, PCs:
      Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode), EOF1(Nmode,Nt))
      !---------------------------------------------------------------------------------
      ! Loop over all drought indices:
      Fname="pca_drou_idx.nc"; Fname1="eof_timeseries_drou_idx.nc"
      Do iidx=1,NDr_idx
         X(1:Nt,1:Nx) = Dr_IDX(1:Nt,1:Nx,iidx)
         ! Tinh cac EOFs, PCs
         Call Covar(Nt,Nx)
         Call jacobi
         Call eigsrt
         print*, SUM(EValues, Mask=EValues>=0)
         K = Nmode
         Call Calc_AX
         if (PC_Stard==1) then 
            Call Calc_EOF_PCA_Standard
         else
            Call Calc_EOF_PCA
         endif
         !------------------------------------------------------------------------------
         ! Write results to files:
         ! Fname="pca_drou_idx.nc"; Fname1="eof_timeseries_drou_idx.nc"
         Call Write_EOF_PCA (Fname,Fname1,iidx)
         print*, "===========> Complete for ",IDX_nm(iidx)
      Enddo
      !Deallocate(Xbar,Std,Rx,EValues,V,Fx,Ax,D,EVtor,Lamda,EOF,PCs,ECumulate)
      !Enddo
      print*,"Program is normal terminating ..."
      CONTAINS
      !---------------------------------------------------------------------------------
      Subroutine Read_Drough_IDX (Fname)
         Integer :: ist, ilon, ilat, it, Nt, imon, imon1, iyr
         Integer :: Nyear, Nmon, nlat1, nlat0
         Character (Len=*) :: Fname
         Real :: tmp, xmean(12), ic(12)
         Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:)

         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, PDSI, "ori_pdsi")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "SPEI_12")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "SPI_12")
         Deallocate (xlon, xlat)
         Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land")
         Nyear=30 ! 1976-2005
         Nmon =12
         Ntime = Nt - Nmon

         !Nlat1 = nlat
         !Nlat0 = 0 !15
         !Nlat  = nlat-Nlat0      ! Tru vung Australia (1-15 ~ 15S-11S)
         !Allocate (X0(nlat)); X0(1:nlat)=xlat(Nlat0+1:Nlat1)
         !Deallocate (xlat); Allocate (xlat(nlat))
         !xlat(:) = X0(:)
         !Deallocate (X0)

         Allocate (Mask_all(nlon,nlat))
         Mask_all=Mask

         Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime),X3(nlon, nlat, Ntime))
         X1(:,1:nlat,1:Ntime)=PDSI(:,Nlat0+1:Nlat1,Nmon+1:Nt)
         X2(:,1:nlat,1:Ntime)=SPEI12(:,Nlat0+1:Nlat1,Nmon+1:Nt)
         X3(:,1:nlat,1:Ntime)=SPI12(:,Nlat0+1:Nlat1,Nmon+1:Nt)

         Deallocate (PDSI,SPEI12,SPI12)
         Allocate (PDSI(nlon,nlat,Ntime),SPEI12(nlon,nlat,Ntime),SPI12(nlon,nlat,Ntime))
         PDSI=X1; SPEI12=X2; SPI12=X3

         Deallocate (X1, X2, X3)

         ist=0
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               if (Mask(ilon,ilat)>0) then
                  if (.not.Empty_Ar1(PDSI  (ilon,ilat,1:Ntime), Ntime).AND. &
                      .not.Empty_Ar1(SPI12 (ilon,ilat,1:Ntime), Ntime).AND. &
                      .not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then
                     ist=ist+1
                  else
                     Mask_all(ilon,ilat)=Miss
                  endif
               endif
            Enddo
         Enddo
         Nst=ist
         Allocate (Dr_IDX(Ntime,Nst,NDr_idx))
         Dr_IDX(:,:,:) = Miss
         ist=0
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               if (Mask(ilon,ilat)>0) then
                  if (.not.Empty_Ar1(PDSI  (ilon,ilat,1:Ntime), Ntime).AND. &
                      .not.Empty_Ar1(SPI12 (ilon,ilat,1:Ntime), Ntime).AND. &
                      .not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then
                     ist=ist+1
                     Dr_IDX(1:Ntime,ist,1)=PDSI  (ilon,ilat,1:Ntime)
                     Dr_IDX(1:Ntime,ist,2)=SPI12 (ilon,ilat,1:Ntime)
                     Dr_IDX(1:Ntime,ist,3)=SPEI12(ilon,ilat,1:Ntime)
                  endif
               endif
            Enddo
         Enddo
         Deallocate (PDSI,SPEI12,SPI12)
         print*,"Reading data ... Done !"
      End Subroutine Read_Drough_IDX
      !---------------------------------------------------------------------------------
      Logical Function Empty_Ar1(X,N)
         Integer :: N,i
         Real :: X(N)
         Logical :: Empty
         Empty=.False.
         Do i=1,N
            if (X(i)==Miss.OR.IsNaN(X(i))) then
               Empty=.True.; Exit
            endif
         Enddo
         Empty_Ar1=Empty
      End Function Empty_Ar1
      !---------------------------------------------------------------------------------
      Subroutine Write_EOF_PCA (F_pca, F_eof, iidx)
         Character (Len=*) :: F_eof, F_pca
         Character (Len=255) :: var_ln, t_unit
         Integer :: i,j,k,iidx
         Real :: tmp, Xtmp(Nmode)
         !Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:)

         !-------------------------------------------------------------------------         
         PCA(:,:,:)=Miss; k=0
         Do j=1, Nlat
            Do i=1, Nlon
               if (Mask_all(i,j)>0) then
                  k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode)
               endif
            Enddo
         Enddo
         Do i=1, Nmode; xmode(i)=i; Enddo         ! PCA's modes
         var_ln="Spatial distribution of PCAs of "//trim(IDX_nm(iidx))
         t_unit="months since 1976-01-01 "
         if (iidx==1) then
            Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, &
              trim(IDX_nm(iidx)), trim(t_unit),trim(var_ln), &
               "-", "lon", "lat","mode", Miss)
         else
            Call wr_ncdf_3d_apd (F_pca,PCA,trim(IDX_nm(iidx)),nlon,nlat, &
                 Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss)
         endif
         !-------------------------------------------------------------------------
         ECumulate(1:Nmode) = EValues(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, ECumulate, &
                 trim(IDX_nm(iidx))//"_eig", Nmode, &
                "Eigen values of "//trim(IDX_nm(iidx)), "-", "mode", Miss)
         !X(1:nmodes) = variance(1:nmodes)
         !Call wr_ncdf_1d_apd (F_pca, X, trim(IDX_nm(iidx))//"_var",nmodes, &
         !       "Variances of "//trim(IDX_nm(iidx)), "-", "mode", Miss )

         tmp=0.0
         Do i=1,Mcount
            tmp=tmp+EValues(i)
            if (I<=Nmode) then
               ECumulate(i)=0
               Do j=1,I
                  ECumulate(i)=ECumulate(i)+EValues(j)
               Enddo
            endif
         Enddo
         Do i=1,Nmode
            ECumulate(i)=ECumulate(i)/tmp
         Enddo

         Call wr_ncdf_1d_apd (F_pca, ECumulate, &
                trim(IDX_nm(iidx))//"_cum",Nmode, &
                "Cumuluted Eigenvalues of "//trim(IDX_nm(iidx)), &
                "-", "mode", Miss )
         !-------------------------------------------------------------------------
         Do i=1, Nmode; EOF1(i,1:Nt)=EOF(1:Nt,i); Enddo
         var_ln="Timeseries of eofs of "//trim(IDX_nm(iidx))
         if (iidx==1) then
            Call wr_ncdf_2d (F_eof,xmode,xtime, Nmode,Nt, EOF1, &
                     trim(IDX_nm(iidx)), &
                     trim(var_ln), "-", "mode", "mon", Miss)
         else
            Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nmode,Nt, &
                     trim(var_ln), "-", "mode", "mon", Miss)
         endif
         print*,"Write_EOF_PCA: Done !!!"
      End Subroutine Write_EOF_PCA

      !---------------------------------------------------------------------------------
      Subroutine Calc_AX
      Integer :: i,j
      Real :: SQ
      DO j=1,K
         SQ=Sqrt(EValues(j))
         DO i=1,Mcount
            Ax(I,J)=EVtor(i,j)*SQ
         ENDDO
      ENDDO
      End Subroutine Calc_AX
      !---------------------------------------------------------------------------------
      Subroutine Calc_EOF_PCA_Standard
      Integer :: i,j
      Real*8 :: Modul_Fx
        if (switch==0) then
           Ax(1:Nx,1:K)=EVtor(1:Nx,1:K)
           Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K))
           Do j=1,K
              !Fx(:,j)=Fx(:,j)/(EValues(j))/Factor ! Fx Chuan hoa
              Fx(:,j)=Fx(:,j)/Factor               ! Voi PCs chuan hoa (Ax=EVtor)
           Enddo
           EOF(1:Nt,1:K) = Fx(1:Nt,1:K)
           PCs(1:Nx,1:K) = Ax(1:Nx,1:K)
        else if (switch==1) then
           Do j=1, K
              !Ax(:,j)=Ax(:,j)/sqrt(EValues(j))
              Lamda(j,j)=Sqrt(EValues(j))
           Enddo
           YT = Transpose(Y)
           !Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,1:K))
           Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),EVtor(1:Nt,1:K))
           Do j=1,K
              Modul_Fx=0
              Do i=1, Ncount
                 Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j)
              Enddo
              Modul_Fx = Sqrt(Modul_Fx)
              Fx(:,j)=Fx(:,j)/Modul_Fx
           Enddo
           !Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K))
           EOF(1:Nt,1:K) = Ax(1:Nt,1:K)
           PCs(1:Nx,1:K) = Fx(1:Ncount,1:K)
        endif
      End Subroutine Calc_EOF_PCA_Standard
      !---------------------------------------------------------------------------------
      Subroutine Calc_EOF_PCA
      Integer :: i,j
      Real*8 :: Modul_Fx
        if (switch==0) then
           Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K))
           Do j=1,K
              Fx(:,j)=Fx(:,j)/(EValues(j))/Factor
           Enddo
           EOF(1:Nt,1:K) = Fx(1:Nt,1:K)
           PCs(1:Nx,1:K) = Ax(1:Nx,1:K)
        else if (switch==1) then
           Do j=1, K
              Ax(:,j)=Ax(:,j)/sqrt(EValues(j))
              Lamda(j,j)=Sqrt(EValues(j))
           Enddo
           YT = Transpose(Y)
           Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,1:K))
           Do j=1,K
              Modul_Fx=0
              Do i=1, Ncount
                 Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j)
              Enddo
              Modul_Fx = Sqrt(Modul_Fx)
              Fx(:,j)=Fx(:,j)/Modul_Fx
           Enddo
           Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K))
           EOF(1:Nt,1:K) = Ax(1:Nt,1:K)
           PCs(1:Nx,1:K) = Fx(1:Ncount,1:K)
        endif
        print*,"Calculating EOF & PCs ... Done !"
      End Subroutine Calc_EOF_PCA
      !---------------------------------------------------------------------------------
      SUBROUTINE Covar(Nt,Nx)
      Integer :: Nt, Nx, N1, L, J, K, IO, I
      Real*8 :: Sx, Xtb
      !---------------------------------------------------
      Where (IsNaN(X)); X=Miss; Endwhere
      Do j=1, Nx
         Sx=0; Xtb=0; N1=0
         Do i=1, Nt
            if (X(i,j)/=Miss) then
               n1=n1+1; Xtb=Xtb+X(i,j); Sx=Sx+X(i,j)*X(i,j)
            Endif
         Enddo
         if (n1>nt/2) then
            XBAR(j)=Xtb/real(n1); STD(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j))
            if (IsNaN(Xbar(j))) print*,"XBAR NaN",j
         endif
         Do i=1, Nt
            if (X(i,j)==Miss) then
               X(i,j)=XBAR(j) ! Fill-up Miss by mean values
            endif
            Y(i,j) = (X(i,j)-XBAR(j))
            if (standard==1) Y(i,j) = Y(i,j)/STD(j)
         Enddo
      Enddo
      Where (IsNaN(Y)); Y=Miss; Endwhere
      if (switch==1) then
         RR = MatMul(Y,Transpose(Y))
      else if (switch==0) then
         RR = MatMul(Transpose(Y),Y)
      endif
      RR(:,:)=RR(:,:)/real(Nt)
      ! Store R in 1D-Array of Rx to be returned
      Do I=1,Mcount
         DO J=1,Mcount
            if ((I-J) < 0) then
               L=I+(J*J-J)/2
            else
               L=J+(I*I-I)/2
            endif
            RX(L)=RR(i,j)
         ENDDO
      ENDDO
      print*,"Calculating Covariance ... Done !"
      End SUBROUTINE Covar
      !---------------------------------------------------------------------------------
      SUBROUTINE eigsrt 
      INTEGER :: n,np
      Real*8, Allocatable :: D(:), V(:,:)
      INTEGER :: i,j,k
      REAL*8 :: p

      Allocate (D(Mcount),V(Mcount,Mcount))
      N=Mcount
      D(:)=EValues(:); V(:,:)=EVtor(:,:)

      do  i=1,n-1
        k=i
        p=d(i)
        do  j=i+1,n
          if(d(j).ge.p)then
            k=j
            p=d(j)
          endif
        enddo
        if(k.ne.i)then
          d(k)=d(i)
          d(i)=p
          do  j=1,n
            p=v(j,i)
            v(j,i)=v(j,k)
            v(j,k)=p
          enddo
        endif
      enddo
      EValues(:)=D(:)
      EVtor(:,:)=V(:,:)
      Deallocate (D, V)
      return
      END SUBROUTINE eigsrt
      !---------------------------------------------------------------------------------
      SUBROUTINE jacobi 
      INTEGER :: n,np,nrot,NMAX
      INTEGER :: i,ip,iq,j
      REAL*8 :: c,g,h,s,sm,t,tau,theta,tresh !,b(NMAX),z(NMAX)
      Real*8, Allocatable :: A(:,:), D(:), V(:,:), B(:), Z(:)

      Allocate (A(Mcount,Mcount),D(Mcount),V(Mcount,Mcount),B(Mcount),Z(Mcount))
      A(:,:) = RR(:,:); N=Mcount
      do  ip=1,n
        do  iq=1,n
          v(ip,iq)=0.
        enddo
        v(ip,ip)=1.
      enddo
      do  ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.
      enddo
      nrot=0
      do  i=1,50
        sm=0.
        do  ip=1,n-1
          do iq=ip+1,n
            sm=sm+abs(a(ip,iq))
          enddo
        enddo
        if(sm.eq.0.) then
          EValues(:)=D(:)
          EVtor(:,:)=V(:,:)
          Deallocate (A, D, V, B, Z)
          return
        endif
        if(i.lt.4)then
          tresh=0.2*sm/n**2
        else
          tresh=0.
        endif
        do  ip=1,n-1
          do  iq=ip+1,n
            g=100.*abs(a(ip,iq))
            if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))).and. &
               (abs(d(iq))+g.eq.abs(d(iq))))then
              a(ip,iq)=0.
            else if(abs(a(ip,iq)).gt.tresh)then
              h=d(iq)-d(ip)
              if(abs(h)+g.eq.abs(h))then
                t=a(ip,iq)/h
              else
                theta=0.5*h/a(ip,iq)
                t=1./(abs(theta)+sqrt(1.+theta**2))
                if(theta.lt.0.)t=-t
              endif
              c=1./sqrt(1+t**2)
              s=t*c
              tau=s/(1.+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0.
              do  j=1,ip-1
                g=a(j,ip)
                h=a(j,iq)
                a(j,ip)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              enddo
              do  j=ip+1,iq-1
                g=a(ip,j)
                h=a(j,iq)
                a(ip,j)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              enddo
              do  j=iq+1,n
                g=a(ip,j)
                h=a(iq,j)
                a(ip,j)=g-s*(h+g*tau)
                a(iq,j)=h+s*(g-h*tau)
              enddo
              do  j=1,n
                g=v(j,ip)
                h=v(j,iq)
                v(j,ip)=g-s*(h+g*tau)
                v(j,iq)=h+s*(g-h*tau)
              enddo
              nrot=nrot+1
            endif
          enddo
        enddo
        do  ip=1,n
          b(ip)=b(ip)+z(ip)
          d(ip)=b(ip)
          z(ip)=0.
        enddo
      enddo
      Deallocate (A, D, V, B, Z)
      !pause 'too many iterations in jacobi'
      return
      END SUBROUTINE jacobi

      !----------------------------------------------------------------------------------!
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"

      !---------------------------------------------------------------------------------
      END Program